dapx_gca: Approximate density and distribution via Gram-Charlier A... In PDQutils: PDQ Functions via Gram Charlier, Edgeworth, and Cornish Fisher Approximations

Description

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

Usage

 ```1 2 3 4 5 6 7``` ```dapx_gca(x, raw.moments, support=NULL, basis=c('normal','gamma','beta','arcsine','wigner'), basepar=NULL, log=FALSE) papx_gca(q, raw.moments, support=NULL, basis=c('normal','gamma','beta','arcsine','wigner'), basepar=NULL, lower.tail=TRUE, log.p=FALSE) ```

Arguments

 `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 `c(-Inf,Inf)` for the normal basis, `c(0,Inf)` for the gamma basis, and `c(0,1)` for the Beta, and `c(-1,1)` for the arcsine and wigner. `basis` the basis under which to perform the approximation. `'normal'` gives the classical 'A' series expansion around the PDF and CDF of the normal distribution via Hermite polynomials. `'gamma'` expands around a gamma distribution with parameters `basepar\$shape` and `basepar\$scale`. `'beta'` expands around a beta distribution with parameters `basepar\$shape1` and `basepar\$shape2`. `basepar` the parameters for the base distribution approximation. If `NULL`, the shape and rate are inferred from the first two moments and/or from the `support` as appropriate. `log` logical; if TRUE, densities f are given as log(f). `q` where to evaluate the approximate distribution. `log.p` logical; if TRUE, probabilities p are given as log(p). `lower.tail` whether to compute the lower tail. If false, we approximate the survival function.

Details

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).

Value

The approximate density at `x`, or the approximate CDF at

Note

Monotonicity of the CDF is not guaranteed.

Author(s)

Steven E. Pav [email protected]

References

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`

Examples

 ``` 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) ```

Example output

```Loading required package: polynom
[1] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
[38] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
[75] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
[112] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
[149] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
[186] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
[223] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
[260] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
[297] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
[334] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
[371] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
[408] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
[445] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
[482] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
[1] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
[38] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
[75] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
[112] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
[149] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
[186] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
[223] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
[260] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
[297] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
[334] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
[371] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
[408] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
[445] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
[482] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
```

PDQutils documentation built on May 30, 2017, 3:42 a.m.