[R] R codes for g-and-h distribution

[Ben Bolker]

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