[R] R codes for g-and-h distribution
Ben Bolker
bolker at ufl.edu
Fri Jul 27 21:28:30 CEST 2007
filame uyaco <filams0704 <at> yahoo.com> writes:
>
>
> hi!
>
> I would like to ask help how to generate numbers from g-and-h distribution.
This distribution is like
> normal distribution but span more of the kurtosis and skewness plane. Has R
any package on how to generate
> them?
>
Someone else asked for this in 2005, but I didn't see any
answers. If the wiki weren't down I would put this up there ...
## http://www.itl.nist.gov/div898/software/dataplot/refman2/auxillar/ghpdf.htm
## G(p,g,h) = exp((g*Zp)-1)*exp((h*Zp^2)/2)/g
[n.b. formatting on this page is very confusing, I inserted
parentheses]
## with Zp denoting the normal percent point function of p. When g = 0
## and h = 0, the g-and-h distribution reduces to a standard normal
## distribution.
## The g-and-h probability density function is computed by taking the
## numerical derivative of the cumulative distribution function (which is
## turn computed by numerically inverting the percent point function
## using a bisection method).
## [from
## http://www.itl.nist.gov/div898/handbook/eda/section3/eda362.htm:
## The percent point function (ppf) is the inverse of the cumulative
## distribution function. For this reason, the percent point function
## is also commonly referred to as the inverse distribution function.
## thus R's "q" functions are percent point functions
qgh <- function(q,g,h) {
Zp <- qnorm(q)
## not vectorized!
if (g==0) Zp else (exp(g*Zp)-1)*exp((h*Zp^2/2))/g
}
## since the quantile function is defined, it makes generating
## random values easy!
rgh <- function(n,g,h) {
qgh(runif(n),g,h)
}
## eps sets distance from 1 for searching
## strategy could probably be improved (first approx. based on
## qnorm??)
pgh <- function(p,g,h,eps=1e-7) {
uniroot(function(z) { qgh(z,g,h) - p}, interval=c(eps,1-eps))$root
}
dgh <- function(x,g,h,log=FALSE,ndep=1e-3,...) {
## crude vectorization in x (not g or h)
if (length(x)>1) return(sapply(x,dgh,g=g,h=h,log=log,ndep=ndep,...))
r <- (pgh(x+ndep,g,h)-pgh(x,g,h))/ndep
if (log) log(r) else r
}
## examples ...
set.seed(1001)
## should be approx. normal
r1 = rgh(10000,1e-5,0)
r2 = rgh(10000,1e-5,0)
r3 = rgh(10000,1e-5,0)
## plot 3 different samples
plot(density(r1))
lines(density(r2))
lines(density(r3))
curve(dnorm(x),col=2,add=TRUE)
curve(dgh(x,1e-5,0),col=4,add=TRUE)
## some slight numerical glitches -- e.g. around
## x=-2.6
r4 = rgh(50000,0.4,0.4)
plot(density(r4,n=1024),xlim=c(-10,10))
curve(dnorm(x),add=TRUE,col=2)
curve(dgh(x,0.4,0.4),add=TRUE,col=4)
legend("topleft",c("density","g-h","normal"),
lty=rep(1,3),col=c(1,4,2))
## note glitch again, this time around
## y=-3.89
More information about the R-help
mailing list