Approximate the probability density or cumulative distribution function of a distribution via its raw moments.

1 2 3 4 5 6 7 |

`x` |
where to evaluate the approximate density. |

`raw.moments` |
an atomic array of the 1st through kth raw moments of the probability distribution. |

`support` |
the support of the density function. It is assumed
that the density is zero on the complement of this open interval.
This defaults to |

`basis` |
the basis under which to perform the approximation. |

`basepar` |
the parameters for the base distribution approximation.
If |

`log` |
logical; if TRUE, densities |

`q` |
where to evaluate the approximate distribution. |

`log.p` |
logical; if TRUE, probabilities p are given
as |

`lower.tail` |
whether to compute the lower tail. If false, we approximate the survival function. |

Given the raw moments of a probability distribution, we can approximate the probability
density function, or the cumulative distribution function, via a Gram-Charlier
expansion on the standardized distribution. This expansion uses
some weighting function, *w*, typically the density of some 'parent'
probability distribution, and polynomials, *p_n* which are
orthogonal with respect to that weighting:

*integral_-inf^inf p_n(x) p_m(x) w(x) dx = h_n delta_mn.*

Let *f(x)* be the probability density of some random variable,
with cumulative distribution function *F(x)*. We express

*f(x) = sum_n c_n p_n(x) w(x)*

The constants *c_n* can be computed from the known moments
of the distribution.

For the Gram Charlier 'A' series, the weighting function is the PDF of the normal distribution, and the polynomials are the (probabilist's) Hermite polynomials. As a weighting function, one can also use the PDF of the gamma distribution (resulting in generalized Laguerre polynomials), or the PDF of the Beta distribution (resulting in Jacobi polynomials).

The approximate density at `x`

, or the approximate CDF at

Monotonicity of the CDF is not guaranteed.

Steven E. Pav shabbychef@gmail.com

Jaschke, Stefan R. "The Cornish-Fisher-expansion in the context of Delta-Gamma-normal approximations." No. 2001, 54. Discussion Papers, Interdisciplinary Research Project 373: Quantification and Simulation of Economic Processes, 2001. http://www.jaschke-net.de/papers/CoFi.pdf

S. Blinnikov and R. Moessner. "Expansions for nearly Gaussian distributions." Astronomy and Astrophysics Supplement 130 (1998): 193-205. http://arxiv.org/abs/astro-ph/9711239

`qapx_cf`

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 | ```
# normal distribution:
xvals <- seq(-2,2,length.out=501)
d1 <- dapx_gca(xvals, c(0,1,0,3,0), basis='normal')
d2 <- dnorm(xvals)
# they should match:
d1 - d2
qvals <- seq(-2,2,length.out=501)
p1 <- papx_gca(qvals, c(0,1,0,3,0))
p2 <- pnorm(qvals)
p1 - p2
xvals <- seq(-6,6,length.out=501)
mu <- 2
sigma <- 3
raw.moments <- c(2,13,62,475,3182)
d1 <- dapx_gca(xvals, raw.moments, basis='normal')
d2 <- dnorm(xvals,mean=mu,sd=sigma)
## Not run:
plot(xvals,d1)
lines(xvals,d2,col='red')
## End(Not run)
p1 <- papx_gca(xvals, raw.moments, basis='normal')
p2 <- pnorm(xvals,mean=mu,sd=sigma)
## Not run:
plot(xvals,p1)
lines(xvals,p2,col='red')
## End(Not run)
# for a one-sided distribution, like the chi-square
chidf <- 30
ords <- seq(1,9)
raw.moments <- exp(ords * log(2) + lgamma((chidf/2) + ords) - lgamma(chidf/2))
xvals <- seq(0.3,10,length.out=501)
d1g <- dapx_gca(xvals, raw.moments, support=c(0,Inf), basis='gamma')
d2 <- dchisq(xvals,df=chidf)
## Not run:
plot(xvals,d1g)
lines(xvals,d2,col='red')
## End(Not run)
p1g <- papx_gca(xvals, raw.moments, support=c(0,Inf), basis='gamma')
p2 <- pchisq(xvals,df=chidf)
## Not run:
plot(xvals,p1g)
lines(xvals,p2,col='red')
## End(Not run)
# for a one-sided distribution, like the log-normal
mu <- 2
sigma <- 1
ords <- seq(1,8)
raw.moments <- exp(ords * mu + 0.5 * (sigma*ords)^2)
xvals <- seq(0.5,10,length.out=501)
d1g <- dapx_gca(xvals, raw.moments, support=c(0,Inf), basis='gamma')
d2 <- dnorm(log(xvals),mean=mu,sd=sigma) / xvals
## Not run:
plot(xvals,d1g)
lines(xvals,d2,col='red')
## End(Not run)
``` |

Questions? Problems? Suggestions? Tweet to @rdrrHQ or email at ian@mutexlabs.com.

All documentation is copyright its authors; we didn't write any of that.