PDQutils (version 0.1.6)

dapx_gca: Approximate density and distribution via Gram-Charlier A expansion.

Description

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

Usage

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 \(\mbox{log}(f)\).

q

where to evaluate the approximate distribution.

log.p

logical; if TRUE, probabilities p are given as \(\mbox{log}(p)\).

lower.tail

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

Value

The approximate density at x, or the approximate CDF at

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: $$\int_{-\infty}^{\infty} p_n(x) p_m(x) w(x) \mathrm{d}x = 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 \ge 0} 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).

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

See Also

qapx_cf

Examples

# 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')
## ---------------------------------------------
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')
## ---------------------------------------------

# 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')
## ---------------------------------------------

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')
## ---------------------------------------------

# 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')
## ---------------------------------------------