gk R package

July 12, 2011

I’ve started tidying up some code from my thesis to create an R package called “gk”, containing functions for the g-and-k distribution.  This defined by its quantile function (i.e. its inverse cdf) Q(\Phi^{-1}(u)) where \Phi^{-1} is the standard normal cdf and
Q(z) = A + B \left( 1 + c \frac{1-\exp(-gz)}{1+\exp(-gz)} \right) (1+z^2)^k z \qquad \text(1)
To explain, if u is a random variable with a uniform distribution on [0,1], then z=\Phi^{-1}(u) follows a standard normal distribution and Q(z) follows the g-and-k distribution. Indeed this is a straightforward method (the inversion method) for simulation.

However, the cdf and density functions are not available in closed form making likelihood based inference hard (although possible through numerical inversion of (1)).  Inference by ABC methods is possible, and this has been much studied, both from interest in using the g-and-k distribution and as a testing ground for research in ABC methods (References include Allingham et al, Drovandi and Pettitt, and McVinish.)

The package contains distributional functions in the standard R format, plus some helper functions for, and applications of these.  These functions should be useful for exploratory implementations of ABC methods.  No such code is included (for now!) as my thesis ABC g-and-k code is messy and uses a lot of C routines.  One application function that is included is gkMLE which estimates MLEs under a model of iid g-and-k observations.

A few technical notes.  The helper functions z2gk and gk2z are the main workhorse functions.  z2gk converts standard normal quantiles to g-and-k quantiles by equation (1).  gk2z performs the inversion operation numerically (A technical point is that this required providing upper and lower bounds on Q^{-1}(x).  Some rough work detailing the bounds I used is here.)  The pgk dgk and gkMLE functions rely on gk2z so they are subject to numerical error and only limited confidence should be placed on their accuracy.