GIG

The Generalized Inverse Gaussian is the parent of normix’s positive distributions: an exponential family on \((0, \infty)\) whose log-partition is Bessel-valued. Its density with parameters \((p, a, b)\) is

\[ f(x \mid p, a, b) = \frac{(a/b)^{p/2}}{2\,K_p(\sqrt{ab})}\, x^{p-1}\exp\!\Big(-\tfrac12\big(a x + b/x\big)\Big), \qquad x > 0, \]

where \(K_p\) is the modified Bessel function of the second kind. It nests the Gamma (\(b \to 0\)), InverseGamma (\(a \to 0\)), and InverseGaussian (\(p = -\tfrac12\)).

Density gallery

For \(a > 0\) and \(b > 0\) the density tends to 0 at the origin and has one interior mode for every real \(p\). The parameters \(a\) and \(b\) set the scale and the two tails. The monotone Gamma shape is the limit \(b \to 0\) with \(p \le 1\).

Show code cell source

Hide code cell source

import jax
jax.config.update("jax_enable_x64", True)
import jax.numpy as jnp
import numpy as np
import matplotlib.pyplot as plt

from normix import GIG
from normix.utils.plotting import set_theme

set_theme()
np.set_printoptions(precision=5, suppress=True)

x = jnp.linspace(1e-3, 6.0, 400)
fig, ax = plt.subplots()
for p in [-1.0, 0.5, 1.5, 3.0]:
    g = GIG(p=jnp.array(p), a=jnp.array(2.0), b=jnp.array(1.5))
    ax.plot(np.asarray(x), np.asarray(jax.vmap(g.pdf)(x)), label=f"p = {p:g}")
ax.set_xlabel("x"); ax.set_ylabel("density")
ax.set_title("GIG densities (a = 2, b = 1.5)")
ax.legend()
plt.show()

Parametrizations

Stored attributes: p, a, b. As an exponential family with sufficient statistic \(t(x) = [\log x,\; 1/x,\; x]\):

Parametrization	Value
Classical	\(p \in \mathbb{R},\; a > 0,\; b > 0\)
Natural \(\theta\)	\([\,p - 1,\; -b/2,\; -a/2\,]\)
Expectation \(\eta = \nabla\psi\)	\([\,\mathbb{E}[\log X],\; \mathbb{E}[1/X],\; \mathbb{E}[X]\,]\)

Inverting \(\eta \mapsto \theta\) is a strictly convex but stiff Bessel problem, solved by an \(\eta\)-rescaled multi-start Newton iteration. See The Generalized Inverse Gaussian Distribution for the derivation.

Quick usage

gig = GIG(p=jnp.array(1.5), a=jnp.array(2.0), b=jnp.array(1.5))

print("mean/var/std:", float(gig.mean()), float(gig.var()), float(gig.std()))
print("pdf(1.0):    ", float(gig.pdf(jnp.array(1.0))))
print("5% quantile: ", float(gig.ppf(jnp.array(0.05))))

samples = gig.rvs(10_000, seed=0)          # Devroye exact sampler
fitted = GIG.fit_mle(samples)              # moment-match + multi-start Newton
print("fit (p,a,b): ", float(fitted.p), float(fitted.a), float(fitted.b))

mean/var/std: 2.0490381056766593 1.6740381056766491 1.2938462449907444
pdf(1.0):     0.4056682816793732
5% quantile:  0.5593202783302175
fit (p,a,b):  1.262973741884533 1.862191533234287 1.7036540588048583

See also

• API: GIG

• Theory: The Generalized Inverse Gaussian Distribution

• Tutorial: The Generalized Inverse Gaussian — Bessel log-partition and the multi-start solver in depth