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 ______________________________________________ [email protected] mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
