NormalInverseGamma

The Normal-Inverse Gamma distribution is the normal variance-mean mixture

\[ X \mid Y \sim \mathcal{N}(\mu + \gamma Y,\; \Sigma Y), \qquad Y \sim \mathrm{InvGamma}(\alpha, \beta), \]

with an InverseGamma subordinator — equivalently the GIG limit \(a \to 0,\; p = -\alpha,\; b = 2\beta\). The polynomial right tail of the mixing law gives the marginal genuinely heavy (power-law) tails, heavier than the VarianceGamma.

Density gallery

As with every mixture, \(\gamma\) controls skewness about \(\mu\).

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 NormalInverseGamma, UnivariateNormalInverseGamma
from normix.utils.plotting import set_theme

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

x = jnp.linspace(-6.0, 6.0, 400)
fig, ax = plt.subplots()
for g in [-1.5, 0.0, 1.5]:
    ninvg = UnivariateNormalInverseGamma.from_classical(
        mu=0.0, gamma=g, sigma=1.0, alpha=3.0, beta=2.0)
    ax.plot(np.asarray(x), np.asarray(jax.vmap(ninvg.pdf)(x)), label=rf"$\gamma$ = {g:g}")
ax.set_xlabel("x"); ax.set_ylabel("density")
ax.set_title(r"NormalInverseGamma marginals ($\mu$=0, $\sigma$=1, $\alpha$=3, $\beta$=2)")
ax.legend()
plt.show()

Parametrizations

Built from the shared location/shape block and the InverseGamma subordinator:

Symbol	Attribute	Meaning
\(\mu\)	mu	location \((d,)\)
\(\gamma\)	gamma	skewness \((d,)\)
\(\Sigma = L_\Sigma L_\Sigma^\top\)	L_Sigma	dispersion Cholesky \((d, d)\)
\(\alpha\)	alpha	InverseGamma shape
\(\beta\)	beta	InverseGamma rate

The marginal is not an exponential family, but the joint \((X, Y)\) is — its natural parametrization and the EM machinery are in The Generalized Hyperbolic Distribution.

Quick usage

Raw \((\gamma, \Sigma, \beta)\) are identified only up to the scale gauge \(Y \mapsto cY\) (see EM Algorithm for Generalized Hyperbolic Distributions); compare the invariants \(\mu\), \(\gamma E[Y]\), and \(E[Y]\,\Sigma\) instead.

mu = jnp.array([0.0, 0.0])
gamma = jnp.array([0.3, -0.4])
Sigma = jnp.array([[1.0, 0.3], [0.3, 1.0]])

ninvg = NormalInverseGamma.from_classical(mu=mu, gamma=gamma, sigma=Sigma, alpha=3.0, beta=2.0)
print("mean:", np.asarray(ninvg.mean()))
print("cov:\n", np.asarray(ninvg.cov()))

X = ninvg.rvs(2_000, seed=0)
result = NormalInverseGamma.default_init(X).fit(X, max_iter=50, tol=1e-3)
fit = result.model
ey = float(fit.joint.subordinator().mean())
print("converged:", bool(result.converged), "n_iter:", int(result.n_iter))
print("mu:", np.asarray(fit.mu))
print("gamma * E[Y]:", np.asarray(fit.gamma) * ey)   # true [0.3, -0.4] (E[Y]=1)
print("E[Y] * Sigma:\n", ey * np.asarray(fit.sigma()))
print("mean:", np.asarray(fit.mean()))
print("cov:\n", np.asarray(fit.cov()))

mean: [ 0.3 -0.4]
cov:
 [[1.09 0.18]
 [0.18 1.16]]
converged: True n_iter: 33
mu: [0.0104  0.03065]
gamma * E[Y]: [ 0.27368 -0.39732]
E[Y] * Sigma:
 [[0.94443 0.27664]
 [0.27664 0.9602 ]]
mean: [ 0.28408 -0.36667]
cov:
 [[1.01626 0.17236]
 [0.17236 1.11159]]

The \(d = 1\) sibling UnivariateNormalInverseGamma adds cdf / ppf; the FactorNormalInverseGamma scales to high dimensions.

See also

• API: NormalInverseGamma

• Theory: The Generalized Hyperbolic Distribution

• Tutorial: Normal variance-mean mixtures