Generalized Inverse Gaussian {GeneralizedHyperbolic}

Generalized Inverse Gaussian Distribution

Description

Density function, cumulative distribution function, quantile function and random number generation for the generalized inverse Gaussian distribution with parameter vector param. Utility routines are included for the derivative of the density function and to find suitable break points for use in determining the distribution function.

Usage

dgig(x, chi = 1, psi = 1, lambda = 1,
     param = c(chi, psi, lambda), KOmega = NULL)
pgig(q, chi = 1, psi = 1, lambda = 1,
     param = c(chi,psi,lambda), lower.tail = TRUE,
     ibfTol = .Machine$double.eps^(0.85), nmax = 200)
qgig(p, chi = 1, psi = 1, lambda = 1,
     param = c(chi, psi, lambda), lower.tail = TRUE,
     method = c("spline", "integrate"),
     nInterpol = 501, uniTol = 10^(-7),
     ibfTol = .Machine$double.eps^(0.85), nmax =200, ...)
rgig(n, chi = 1, psi = 1, lambda = 1,
     param = c(chi, psi, lambda))
rgig1(n, chi = 1, psi = 1, param = c(chi, psi))
ddgig(x, chi = 1, psi = 1, lambda = 1,
      param = c(chi, psi, lambda), KOmega = NULL)

Arguments

x, q

Vector of quantiles.

p

Vector of probabilities.

n

Number of observations to be generated.

chi

A shape parameter that by default holds a value of 1.

psi

Another shape parameter that is set to 1 by default.

lambda

Shape parameter of the GIG distribution. Common to all forms of parameterization. By default this is set to 1.

param

Parameter vector taking the form c(chi, psi, lambda) for rgig, or c(chi, psi) for rgig1.

method

Character. If "spline" quantiles are found from a spline approximation to the distribution function. If "integrate", the distribution function used is always obtained by integration.

lower.tail

Logical. If TRUE, probabilities are P(X\leq x), otherwise as P(X>x).

KOmega

Sets the value of the Bessel function in the density or derivative of the density. See Details.

ibfTol

Value of tolerance to be passed to incompleteBesselK by pgig.

nmax

Value of maximum order of the approximating series to be passed to incompleteBesselK by pgig.

nInterpol

The number of points used in qgig for cubic spline interpolation (see splinefun) of the distribution function.

uniTol

Value of tol in calls to uniroot. See uniroot.

...

Passes arguments to uniroot. See Details.

Details

The generalized inverse Gaussian distribution has density

f(x)=\frac{(\psi/\chi)^{\frac{\lambda}{2}}}% {2K_\lambda(\sqrt{\psi\chi})}x^{\lambda-1}% e^{-\frac{1}{2}\left(\chi x^{-1}+\psi x\right)}

for x>0, where K_\lambda() is the modified Bessel function of the third kind with order \lambda.

The generalized inverse Gaussian distribution is investigated in detail in Jörgensen (1982).

Use gigChangePars to convert from the (\delta,\gamma), (\alpha,\beta), or (\omega,\eta) parameterizations to the (\chi,\psi), parameterization used above.

pgig calls the function incompleteBesselK from the package DistributionUtils to integrate the density function dgig. This can be expected to be accurate to about 13 decimal places on a 32-bit computer, often more accurate. The algorithm used is due to Slavinsky and Safouhi (2010).

Calculation of quantiles using qgig permits the use of two different methods. Both methods use uniroot to find the value of x for which a given q is equal F(x) where F denotes the cumulative distribution function. The difference is in how the numerical approximation to F is obtained. The obvious and more accurate method is to calculate the value of F(x) whenever it is required using a call to pghyp. This is what is done if the method is specified as "integrate". It is clear that the time required for this approach is roughly linear in the number of quantiles being calculated. A Q-Q plot of a large data set will clearly take some time. The alternative (and default) method is that for the major part of the distribution a spline approximation to F(x) is calculated and quantiles found using uniroot with this approximation. For extreme values (for which the tail probability is less than 10^{-7}), the integration method is still used even when the method specifed is "spline".

If accurate probabilities or quantiles are required, tolerances (intTol and uniTol) should be set to small values, say 10^{-10} or 10^{-12} with method = "integrate". Generally then accuracy might be expected to be at least 10^{-9}. If the default values of the functions are used, accuracy can only be expected to be around 10^{-4}. Note that on 32-bit systems .Machine$double.eps^0.25 = 0.0001220703 is a typical value.

Generalized inverse Gaussian observations are obtained via the algorithm of Dagpunar (1989).

Value

dgig gives the density, pgig gives the distribution function, qgig gives the quantile function, and rgig generates random variates. rgig1 generates random variates in the special case where \lambda = 1.

ddgig gives the derivative of dgig.

Author(s)

David Scott [e-mail], Richard Trendall, and Melanie Luen.

References

Dagpunar, J.S. (1989). An easily implemented generalised inverse Gaussian generator. Commun. Statist. -Simula., 18, 703–710.

Jörgensen, B. (1982). Statistical Properties of the Generalized Inverse Gaussian Distribution. Lecture Notes in Statistics, Vol. 9, Springer-Verlag, New York.

Slevinsky, Richard M., and Safouhi, Hassan (2010) A recursive algorithm for the G transformation and accurate computation of incomplete Bessel functions. Appl. Numer. Math., In press.

See Also

safeIntegrate, integrate for its shortfalls, splinefun, uniroot and gigChangePars for changing parameters to the (\chi,\psi) parameterization, dghyp for the generalized hyperbolic distribution.

Examples

param <- c(2, 3, 1)
gigRange <- gigCalcRange(param = param, tol = 10^(-3))
par(mfrow = c(1, 2))
curve(dgig(x, param = param), from = gigRange[1], to = gigRange[2],
      n = 1000)
title("Density of the\n Generalized Inverse Gaussian")
curve(pgig(x, param = param), from = gigRange[1], to = gigRange[2],
      n = 1000)
title("Distribution Function of the\n Generalized Inverse Gaussian")
dataVector <- rgig(500, param = param)
curve(dgig(x, param = param), range(dataVector)[1], range(dataVector)[2],
      n = 500)
hist(dataVector, freq = FALSE, add = TRUE)
title("Density and Histogram\n of the Generalized Inverse Gaussian")
DistributionUtils::logHist(dataVector, main =
        "Log-Density and Log-Histogram\n of the Generalized Inverse Gaussian")
curve(log(dgig(x, param = param)), add = TRUE,
      range(dataVector)[1], range(dataVector)[2], n = 500)
par(mfrow = c(2, 1))
curve(dgig(x, param = param), from = gigRange[1], to = gigRange[2],
      n = 1000)
title("Density of the\n Generalized Inverse Gaussian")
curve(ddgig(x, param = param), from = gigRange[1], to = gigRange[2],
      n = 1000)
title("Derivative of the Density\n of the Generalized Inverse Gaussian")