dgig function

Generalized Inverse Gaussian Distribution

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

dgig(x, Theta, KOmega = NULL)
pgig(q, Theta, small = 10^(-6), tiny = 10^(-10), deriv = 0.3,
     subdivisions = 100, accuracy = FALSE, ...)
qgig(p, Theta, small = 10^(-6), tiny = 10^(-10), deriv = 0.3,
     nInterpol = 100, subdivisions = 100, ...)
rgig(n, Theta)
rgig1(n, Theta)
ddgig(x, Theta, KOmega = NULL, ...)
gigBreaks (Theta, small = 10^(-6), tiny = 10^(-10), deriv = 0.3, ...)

Arguments

• x, q: Vector of quantiles.

• p: Vector of probabilities.

• n: Number of observations to be generated.

• Theta: Parameter vector taking the form c(lambda,chi,psi)

  for rgig, or c(chi,psi) for rgig1.

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

• small: Size of a small difference between the distribution function and zero or one. See Details .

• tiny: Size of a tiny difference between the distribution function and zero or one. See Details .

• deriv: Value between 0 and 1. Determines the point where the derivative becomes substantial, compared to its maximum value. See Details .

• accuracy: Uses accuracy calculated by integrate

  to try and determine the accuracy of the distribution function calculation.

• subdivisions: The maximum number of subdivisions used to integrate the density returning the distribution function.

• nInterpol: The number of points used in qhyperb for cubic spline interpolation (see splinefun) of the distribution function.

• ...: 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)}%f(x)=(psi/chi)^{lambda/2}/%(2K_lambda(sqrt(psi chi)))x^(lambda-1)%exp(-(1/2)(chi x^(-1)+psi x))

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 (1982).

Use gigChangePars to convert from the $(delta,gamma)$, $(alpha,beta)$, or $(omega,eta)$ parameterisations to the $(chi,psi)$, parameterisation used above.

pgig breaks the real line into eight regions in order to determine the integral of dgig. The break points determining the regions are found by gigBreaks, based on the values of small, tiny, and deriv. In the extreme tails of the distribution where the probability is tiny according to gigCalcRange, the probability is taken to be zero. For the generalized inverse Gaussian distribution the leftmost breakpoint is not affected by the value of tiny but is always taken as 0. In the inner part of the distribution, the range is divided in 6 regions, 3 above the mode, and 3 below. On each side of the mode, there are two break points giving the required three regions. The outer break point is where the probability in the tail has the value given by the variable small. The inner break point is where the derivative of the density function is deriv times the maximum value of the derivative on that side of the mode. In each of the 6 inner regions the numerical integration routine safeIntegrate (which is a wrapper for integrate) is used to integrate the density dgig.

qgig use the breakup of the real line into the same 8 regions as pgig. For quantiles which fall in the 2 extreme regions, the quantile is returned as -Inf or Inf as appropriate. In the 6 inner regions splinefun is used to fit values of the distribution function generated by pgig. The quantiles are then found using the uniroot function.

pgig and qgig may generally be expected to be accurate to 5 decimal places. Unfortunately, when lambda is less than -0.5, the accuracy may be as little as 3 decimal places.

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

Returns

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$. An estimate of the accuracy of the approximation to the distribution function may be found by setting accuracy = TRUE in the call to phyperb

which then returns a list with components value and error.

ddgig gives the derivative of dgig.

gigBreaks returns a list with components: - xTiny: Takes value 0 always.

• xSmall: Value such that probability to the left is less than small.

• lowBreak: Point to the left of the mode such that the derivative of the density is deriv times its maximum value on that side of the mode

• highBreak: Point to the right of the mode such that the derivative of the density is deriv times its maximum value on that side of the mode

• xLarge: Value such that probability to the right is less than small.

• xHuge: Value such that probability to the right is less than tiny.

• modeDist: The mode of the given generalized inverse Gaussian distribution.

References

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

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

Author(s)

David Scott d.scott@auckland.ac.nz , Richard Trendall, and Melanie Luen.

See Also

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

Examples

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