# [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

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