Distributions

Each mixture family (Variance Gamma, Normal-Inverse Gamma, Normal-Inverse Gaussian, Generalized Hyperbolic) documents four classes from one module: the joint \((X,Y)\) law, the marginal-of-\(X\) law, a scalar (\(d=1\)) convenience wrapper, and the factor-analysis variant.

Univariate Distributions

Gamma

Gamma distribution as an exponential family.

\[p(x \mid \alpha, \beta) = \frac{\beta^\alpha}{\Gamma(\alpha)} x^{\alpha-1} e^{-\beta x}, \quad x > 0\]

Exponential family structure:

\[h(x) = 1, \quad t(x) = [\log x,\; x]\]

\[\theta = [\alpha-1,\; -\beta], \quad \theta_1 > -1,\; \theta_2 < 0\]

\[\psi(\theta) = \log\Gamma(\theta_1+1) - (\theta_1+1)\log(-\theta_2)\]

\[\eta = [\psi(\alpha) - \log\beta,\; \alpha/\beta] \quad \text{(digamma, mean)}\]

class normix.distributions.gamma.Gamma(alpha, beta)

Bases: ExponentialFamily

Gamma(\(\alpha\), \(\beta\)) distribution — shape \(\alpha > 0\), rate \(\beta > 0\).

Parameters:

• alpha (Array)

• beta (Array)

alpha: Array

beta: Array

natural_params()

\(\theta\) from stored classical parameters.

Return type:

Array

static sufficient_statistics(x)

\(t(x)\) for a single unbatched observation.

Parameters:

x (Array)

Return type:

Array

static log_base_measure(x)

\(\log h(x)\) for a single unbatched observation.

Parameters:

x (Array)

Return type:

Array

mean()

E[X]. Subclasses should override with analytical formulas.

Return type:

Array

var()

Var[X]. Subclasses should override with analytical formulas.

Return type:

Array

raw_moment(k)

Raw moment \(E[X^k] = \Gamma(\alpha+k)/(\Gamma(\alpha)\,\beta^k)\).

Finite for all \(k > -\alpha\).

Parameters:

k (Array)

Return type:

Array

raw_moments(ks)

Vectorised raw_moment() over orders ks.

Parameters:

ks (Array)

Return type:

Array

mode()

Mode \((\alpha - 1)/\beta\) for \(\alpha \ge 1\).

For \(\alpha < 1\) the density diverges at \(x=0\) and the formula returns a non-positive value; we clip to LOG_EPS so the return value is always a valid positive sample location.

Return type:

Array

cdf(x)

CDF F(x). Subclasses should override with analytical formulas.

Parameters:

x (Array)

Return type:

Array

ppf(q)

Quantile function \(F^{-1}(q) = \mathrm{gammaincinv}(\alpha, q) / \beta\).

Parameters:

q (Array)

Return type:

Array

rvs(n, seed=42)

Sample n observations from \(\mathrm{Gamma}(\alpha, \beta)\) via JAX PRNG.

Parameters:

• n (int)

• seed (int)

Return type:

Array

classmethod from_natural(theta)

Construct from natural parameters θ. Subclasses must override.

Parameters:

theta (Array)

Return type:

Gamma

to_gig(*, boundary_eps=0.0)

Exact embedding into the GIG family.

Gamma(\(\alpha\), \(\beta\)) is the \(b \to 0\) limit of GIG(\(p = \alpha,\; a = 2\beta,\; b\)). With boundary_eps = 0 the embedding stores b = 0 exactly; pass a small positive value to stay in the strict interior of GIG’s domain (matters only for downstream expectation_params calls on the lifted GIG).

Parameters:

boundary_eps (float)

classmethod from_expectation(eta, *, theta0=None, maxiter=100, tol=1e-12, backend='jax', alpha_min=None, **kwargs)

Closed-form \(\eta \to \theta\) via Newton on \(\psi(\alpha) - \log\alpha = \eta_1 - \log\eta_2\).

\(\eta = [E[\log X],\; E[X]]\) → \(\alpha\) from digamma inversion, \(\beta = \alpha / \eta_2\).

Parameters:

• backend (str) – 'jax' (default): lax.fori_loop Newton (JIT-compatible). 'cpu': scipy.special digamma/polygamma (no XLA tracing).

• alpha_min (float, optional) – Lower bound on the fitted shape \(\alpha\). When set, the digamma-inversion result is clamped to jnp.maximum(alpha, alpha_min) before \(\beta = \alpha/\eta_2\) is formed, so the projected estimate stays self-consistent. This is the opt-in VG estimand control (the ghyp “fix-\(\lambda\)” analogue): it restricts the estimator to a region where the VG marginal likelihood is bounded. None (default) leaves the estimate unconstrained. JIT-safe (no Python branch on traced values).

• eta (Array)

• maxiter (int)

• tol (float)

Return type:

Gamma

Inverse Gamma

InverseGamma distribution as an exponential family.

\[p(x \mid \alpha, \beta) = \frac{\beta^\alpha}{\Gamma(\alpha)} x^{-\alpha-1} e^{-\beta/x}, \quad x > 0\]

Exponential family structure:

\[h(x) = 1, \quad t(x) = [-1/x,\; \log x]\]

\[\theta = [\beta,\; -(\alpha+1)], \quad \theta_1 > 0,\; \theta_2 < -1\]

\[\psi(\theta) = \log\Gamma(-\theta_2-1) - (-\theta_2-1)\log\theta_1 = \log\Gamma(\alpha) - \alpha\log\beta\]

\[\eta = [-\alpha/\beta,\; \log\beta - \psi(\alpha)] \quad (E[-1/X],\; E[\log X])\]

class normix.distributions.inverse_gamma.InverseGamma(alpha, beta)

Bases: ExponentialFamily

InverseGamma(\(\alpha\), \(\beta\)) — shape \(\alpha > 0\), rate \(\beta > 0\).

Parameters:

• alpha (Array)

• beta (Array)

alpha: Array

beta: Array

natural_params()

\(\theta\) from stored classical parameters.

Return type:

Array

static sufficient_statistics(x)

\(t(x)\) for a single unbatched observation.

Parameters:

x (Array)

Return type:

Array

static log_base_measure(x)

\(\log h(x)\) for a single unbatched observation.

Parameters:

x (Array)

Return type:

Array

mean()

E[X]. Subclasses should override with analytical formulas.

Return type:

Array

var()

Var[X]. Subclasses should override with analytical formulas.

Return type:

Array

raw_moment(k)

Raw moment \(E[X^k] = \beta^k\,\Gamma(\alpha-k)/\Gamma(\alpha)\).

Finite for \(k < \alpha\).

Parameters:

k (Array)

Return type:

Array

raw_moments(ks)

Vectorised raw_moment() over orders ks.

Parameters:

ks (Array)

Return type:

Array

mode()

Mode \(\beta / (\alpha + 1)\) (closed form, valid for all \(\alpha > 0\)).

Return type:

Array

cdf(x)

CDF F(x). Subclasses should override with analytical formulas.

Parameters:

x (Array)

Return type:

Array

ppf(q)

Quantile function \(F^{-1}(q) = \beta / \mathrm{gammaincinv}(\alpha, 1-q)\).

Follows from \(F(x) = 1 - P(\alpha, \beta/x) = q\).

Parameters:

q (Array)

Return type:

Array

rvs(n, seed=42)

Sample n observations from \(\mathrm{InvGamma}(\alpha, \beta)\) via JAX PRNG.

Uses the relation: if \(X \sim \mathrm{Gamma}(\alpha, 1/\beta)\) then \(1/X \sim \mathrm{InvGamma}(\alpha, \beta)\).

Parameters:

• n (int)

• seed (int)

Return type:

Array

classmethod from_natural(theta)

Construct from natural parameters θ. Subclasses must override.

Parameters:

theta (Array)

Return type:

InverseGamma

to_gig(*, boundary_eps=0.0)

Exact embedding into the GIG family.

InverseGamma(\(\alpha\), \(\beta\)) is the \(a \to 0\) limit of GIG(\(p = -\alpha,\; a,\; b = 2\beta\)). With boundary_eps = 0 the embedding stores a = 0 exactly; pass a small positive value to stay in the strict interior of GIG’s domain.

Parameters:

boundary_eps (float)

classmethod from_expectation(eta, *, theta0=None, maxiter=100, tol=1e-12, backend='jax', **kwargs)

\(\eta = [-\alpha/\beta,\; \log\beta - \psi(\alpha)]\).

\(\beta = \alpha / (-\eta_1)\); solve \(\psi(\alpha) - \log\alpha = -\eta_2 - \log(-\eta_1)\) via Newton.

Parameters:

• backend (str) – 'jax' (default): lax.fori_loop Newton (JIT-compatible). 'cpu': scipy.special digamma/polygamma (no XLA tracing).

• eta (Array)

• maxiter (int)

• tol (float)

Return type:

InverseGamma

Inverse Gaussian

Inverse Gaussian (Wald) distribution as an exponential family.

\[f(x \mid \mu, \lambda) = \sqrt{\frac{\lambda}{2\pi}}\, x^{-3/2} \exp\!\left(-\frac{\lambda(x-\mu)^2}{2\mu^2 x}\right), \quad x > 0\]

Exponential family structure:

\[h(x) = (2\pi)^{-1/2} x^{-3/2}, \quad t(x) = [x,\; 1/x]\]

\[\theta = \Bigl[-\tfrac{\lambda}{2\mu^2},\; -\tfrac{\lambda}{2}\Bigr], \quad \theta_1 < 0,\; \theta_2 < 0\]

\[ \begin{align}\begin{aligned}\psi(\theta) = -\tfrac{1}{2}\log(-2\theta_2) - \sqrt{(-2\theta_1)(-2\theta_2)}\\\bigl(\tfrac{1}{2}\log(2\pi)\ \text{is absorbed into}\ \log h(x) = -\tfrac{1}{2}\log(2\pi) - \tfrac{3}{2}\log x\bigr)\end{aligned}\end{align} \]

\[\eta = [E[X],\; E[1/X]] = [\mu,\; 1/\mu + 1/\lambda]\]

class normix.distributions.inverse_gaussian.InverseGaussian(mu, lam)

Bases: ExponentialFamily

InverseGaussian(\(\mu\), \(\lambda\)) — mean \(\mu > 0\), shape \(\lambda > 0\).

Parameters:

• mu (Array)

• lam (Array)

mu: Array

lam: Array

natural_params()

\(\theta\) from stored classical parameters.

Return type:

Array

static sufficient_statistics(x)

\(t(x)\) for a single unbatched observation.

Parameters:

x (Array)

Return type:

Array

static log_base_measure(x)

\(\log h(x)\) for a single unbatched observation.

Parameters:

x (Array)

Return type:

Array

log_density_power(alpha)

Log density-power integral via the exact GIG embedding.

InverseGaussian has a non-constant base measure, \(\log h(x) = -\tfrac{1}{2}\log(2\pi) - \tfrac{3}{2}\log x\), so the constant-base-measure formula in ExponentialFamily does not apply. The distribution is exactly \(\mathrm{GIG}(p=-\tfrac12,\, a=\lambda/\mu^2,\, b=\lambda)\), where \(\log h \equiv 0\) and the \(\log x\) term is promoted to a sufficient statistic — so the density-power integral (and hence entropy(), varentropy(), and renyi(), which all route through this method) are evaluated on the GIG embedding.

Parameters:

alpha (Array)

Return type:

Array

mean()

E[X]. Subclasses should override with analytical formulas.

Return type:

Array

var()

Var[X]. Subclasses should override with analytical formulas.

Return type:

Array

raw_moment(k)

Raw moment \(E[X^k]\) via the exact GIG embedding.

Parameters:

k (Array)

Return type:

Array

raw_moments(ks)

Vectorised raw_moment() over orders ks.

Parameters:

ks (Array)

Return type:

Array

mode()

Mode \(\mu\bigl(\sqrt{1 + 9\mu^2/(4\lambda^2)} - 3\mu/(2\lambda)\bigr)\).

Closed-form maximiser of the IG density on \((0, \infty)\).

Return type:

Array

cdf(x)

CDF of the Inverse Gaussian distribution (log-space stable).

\[F(x) = \Phi(t_1) + \exp\!\bigl(2\lambda/\mu + \log\Phi(-t_2)\bigr)\]

where \(t_1 = \sqrt{\lambda/x}\,(x/\mu - 1)\) and \(t_2 = \sqrt{\lambda/x}\,(x/\mu + 1)\). The second term uses log_ndtr to avoid overflow when \(\lambda/\mu\) is large.

Parameters:

x (Array)

Return type:

Array

quantile_table()

Frozen PINV table for amortised ppf() / sampling.

The analytical cdf() does not use this table; hold the return value when evaluating many quantiles at fixed parameters.

Return type:

QuantileTable

ppf(q)

Quantile function via a PINV table built from log_prob().

Trapezoidal-CDF lookup on \(w = \log x\), seeded at \(\log\) mode().

Parameters:

q (Array)

Return type:

Array

rvs(n, seed=42)

Sample n observations from \(\mathrm{InvGaussian}(\mu, \lambda)\) via JAX PRNG.

Uses the algorithm from Michael, Schucany & Haas (1976):

1. \(\nu \sim \mathcal{N}(0,1)\), \(y = \nu^2\)

2. \(x = \mu + \frac{\mu^2 y}{2\lambda} - \frac{\mu}{2\lambda}\sqrt{4\mu\lambda y + \mu^2 y^2}\)

3. \(z \sim \mathrm{Uniform}(0,1)\); return \(x\) if \(z \le \mu/(\mu+x)\), else \(\mu^2/x\)

Uses jnp.where for vectorized branching over the full sample array.

Parameters:

• n (int)

• seed (int)

Return type:

Array

classmethod from_natural(theta)

Construct from natural parameters θ. Subclasses must override.

Parameters:

theta (Array)

Return type:

InverseGaussian

to_gig(*, boundary_eps=0.0)

Exact embedding into the GIG family.

InverseGaussian(\(\mu\), \(\lambda\)) is GIG with \(p = -1/2,\; a = \lambda/\mu^2,\; b = \lambda\). No boundary approximation: the embedding lands strictly inside GIG’s domain. boundary_eps is accepted for API uniformity with Gamma/InvGamma and ignored.

Parameters:

boundary_eps (float)

classmethod from_expectation(eta, *, theta0=None, maxiter=100, tol=1e-12, **kwargs)

Closed-form from \(\eta = [E[X],\; E[1/X]] = [\mu,\; 1/\mu + 1/\lambda]\):

\(\mu = \eta_1\), \(\lambda = 1/(\eta_2 - 1/\eta_1)\).

Parameters:

• eta (Array)

• maxiter (int)

• tol (float)

Return type:

InverseGaussian

Generalized Inverse Gaussian

Generalized Inverse Gaussian (GIG) distribution as an exponential family.

\[f(x \mid p, a, b) = \frac{(a/b)^{p/2}}{2 K_p(\sqrt{ab})} \, x^{p-1} \exp\!\left(-\frac{ax + b/x}{2}\right), \quad x > 0\]

Exponential family structure:

\[h(x) = 1, \quad t(x) = [\log x,\; 1/x,\; x]\]

\[\theta = [p-1,\; -b/2,\; -a/2], \quad \theta_2 \le 0,\; \theta_3 \le 0\]

\[\psi(\theta) = \log 2 + \log K_p(\sqrt{ab}) + \tfrac{p}{2}\log(b/a), \quad p = \theta_1+1,\; a = -2\theta_3,\; b = -2\theta_2\]

\[\eta = [E[\log X],\; E[1/X],\; E[X]]\]

Special cases:

• \(b \to 0,\; p > 0\): GIG → \(\mathrm{Gamma}(p,\; a/2)\)

• \(a \to 0,\; p < 0\): GIG → \(\mathrm{InvGamma}(-p,\; b/2)\)

• \(p = -1/2\): GIG → InverseGaussian

η→θ rescaling (reduces Fisher condition number):

\[s = \sqrt{\eta_2/\eta_3}, \quad \tilde{\eta} = \bigl(\eta_1 + \tfrac{1}{2}\log(\eta_2/\eta_3),\; \sqrt{\eta_2\eta_3},\; \sqrt{\eta_2\eta_3}\bigr)\]

Solve \(\tilde{\eta} \to \tilde{\theta}\) with symmetric GIG (\(\tilde{a} = \tilde{b}\)), then unscale.

Log-Partition Triad Overrides:

• _log_partition_from_theta : JAX, uses log_kv(backend='jax')

• _grad_log_partition : analytical Bessel ratios (5 \(K_\nu\) calls)

• _hessian_log_partition : analytical 11-Bessel Hessian in \(\theta\)-space

• _log_partition_cpu : numpy + log_kv(backend='cpu')

• _grad_log_partition_cpu : analytical Bessel ratios via scipy.kve

• _hessian_log_partition_cpu: central FD on _log_partition_cpu

class normix.distributions.generalized_inverse_gaussian.GeneralizedInverseGaussian(p, a, b)

Bases: ExponentialFamily

Generalized Inverse Gaussian distribution.

Stored: \(p\) (shape, any real), \(a > 0\), \(b > 0\).

Parameters:

• p (Array)

• a (Array)

• b (Array)

p: Array

a: Array

b: Array

to_gig(*, boundary_eps=0.0)

Identity embedding into the GIG family.

GIG is already in GIG coordinates, so this returns self. It exists so that every subordinator family exposes a uniform to_gig() for the shared prior-to-posterior conjugacy map in the EM E-step (see _posterior_gig_params()). boundary_eps is accepted for API uniformity and ignored.

Parameters:

boundary_eps (float)

Return type:

GeneralizedInverseGaussian

natural_params()

\(\theta\) from stored classical parameters.

Return type:

Array

static sufficient_statistics(x)

\(t(x)\) for a single unbatched observation.

Parameters:

x (Array)

Return type:

Array

static log_base_measure(x)

\(\log h(x)\) for a single unbatched observation.

Parameters:

x (Array)

Return type:

Array

static expectation_params_batch(p, a, b, backend='jax')

Vectorized η for arrays of (p, a, b), each shape (N,). Returns (N, 3) array where columns are [E_log_X, E_inv_X, E_X].

backend=’jax’ : vmap over scalar JAX grad backend=’cpu’ : vectorized scipy.kve (6 C-level array calls)

Parameters:

backend (str)

Return type:

Array

mean()

\(E[X] = \eta_3\) from expectation parameters.

Return type:

Array

raw_moment(k)

Raw moment \(E[X^k] = (b/a)^{k/2}\,K_{p+k}(\sqrt{ab})/K_p(\sqrt{ab})\).

Requires \(a, b > 0\). For the Gamma / InverseGamma boundary embeddings use those classes’ own raw_moment() instead.

Parameters:

k (Array)

Return type:

Array

raw_moments(ks)

Vectorised raw_moment() over orders ks (shared \(K_p\)).

Parameters:

ks (Array)

Return type:

Array

var()

\(\mathrm{Var}[X] = E[X^2] - E[X]^2\) via three Bessel evaluations.

Uses the moment formula \(E[X^r] = (b/a)^{r/2}\,K_{p+r}(\sqrt{ab})/K_p(\sqrt{ab})\) rather than the full 11-Bessel Fisher Hessian entry [2, 2].

Return type:

Array

mode()

Interior mode \(\bigl((p-1) + \sqrt{(p-1)^2 + ab}\bigr) / a\).

Closed-form positive critical point of the log-density. For \(p \ge 1\) this is the unique global maximum on \((0,\infty)\); for \(p < 1\) the density diverges at 0 and this returns the interior local maximum.

Return type:

Array

cdf(x)

CDF \(F(x) = P(X \le x)\).

Trapezoidal CDF on a \(w = \log x\) grid built from log_prob(); seeded at \(\log\) mode(). In the degenerate regimes (\(\sqrt{ab} <\) GIG_DEGEN_THRESHOLD) delegates to the limiting Gamma / InverseGamma CDF for accuracy.

JIT-compatible: the degeneracy test uses jax.lax.cond() (no host float() casts).

Parameters:

x (Array)

Return type:

Array

ppf(q)

Quantile (inverse CDF) \(F^{-1}(q)\) via the PINV table.

JIT-compatible: see cdf().

Parameters:

q (Array)

Return type:

Array

quantile_table()

Frozen PINV table (non-degenerate regime only).

Degenerate Gamma / InverseGamma limits bypass the table in cdf() / ppf(); hold this object when evaluating many quantiles at fixed non-degenerate parameters.

Return type:

QuantileTable

rvs(n, seed=42, method='devroye')

Sample n observations from \(\mathrm{GIG}(p, a, b)\).

Parameters:

• n (int) – Sample size.

• seed (int) – Integer seed for JAX PRNG (or scipy random_state for 'scipy').

• method (str) –

  Sampling algorithm:

  • 'devroye' (default) — Transformed density rejection (TDR) on \(\log(x)\), pure JAX, no Bessel functions.

  • 'pinv' — Numerical inverse CDF (CPU table build + JAX sampling), no Bessel. Best for large n with fixed parameters.

  • 'scipy' — scipy.stats.geninvgauss (CPU, original fallback).

Return type:

Array

to_gamma()

KL projection onto the Gamma family.

Minimises \(D_{\mathrm{KL}}(\mathrm{GIG}\,\|\,q)\) over \(q \in \mathrm{Gamma}\) by matching the Gamma sufficient statistics under the source GIG: \(E_{q^*}[\log X] = \eta_1,\; E_{q^*}[X] = \eta_3\). Solved via Gamma.from_expectation().

to_inverse_gamma()

KL projection onto the InverseGamma family.

Matches \(E[-1/X] = -\eta_2,\; E[\log X] = \eta_1\).

to_inverse_gaussian()

KL projection onto the InverseGaussian family.

Matches \(E[X] = \eta_3,\; E[1/X] = \eta_2\). The closed form \(\lambda = 1/(\eta_2 - 1/\eta_3)\) is well-defined whenever \(\eta_2 > 1/\eta_3\) (Jensen, true for every non-degenerate GIG); InverseGaussian.from_expectation() clamps near degenerate inputs.

classmethod from_natural(theta)

Construct from natural parameters θ. Subclasses must override.

Parameters:

theta (Array)

Return type:

GeneralizedInverseGaussian

classmethod from_expectation(eta, *, theta0=None, maxiter=500, tol=1e-10, backend='jax', method='newton', verbose=0)

\(\eta \to \theta\) via \(\eta\)-rescaling + optimization.

Rescaling makes the Fisher matrix symmetric (\(\tilde{a} = \tilde{b}\)), reducing condition number by up to \(10^{30}\) for extreme \(a/b\) ratios.

Parameters:

• theta0 (jax.Array, optional) – Warm-start point (required for JAX solvers; if None, uses multi-start CPU solver with Gamma/InvGamma/InvGauss seeds).

• backend (str) – 'jax' (default, JIT-able) or 'cpu' (scipy, more robust).

• method (str) – 'newton', 'lbfgs', or 'bfgs'.

• eta (Array)

• maxiter (int)

• tol (float)

• verbose (int)

Return type:

GeneralizedInverseGaussian

normix.distributions.generalized_inverse_gaussian.GIG

alias of GeneralizedInverseGaussian

Multivariate Distributions

Multivariate Normal distribution.

Stored: mu (d,), L_Sigma (d×d) lower-triangular Cholesky of \(\Sigma\). All linear algebra via L_Sigma — never form \(\Sigma^{-1}\) explicitly.

Exponential family structure

\[t(x) = [x,\; \operatorname{vec}(xx^\top)], \quad \theta = [\Sigma^{-1}\mu,\; -\tfrac{1}{2}\operatorname{vec}(\Sigma^{-1})], \quad \log h(x) = 0\]

\[\psi(\theta) = \tfrac{1}{2}\mu^\top\Sigma^{-1}\mu - \tfrac{1}{2}\log|\Sigma^{-1}| + \tfrac{d}{2}\log(2\pi)\]

where \(\operatorname{vec}\) uses row-major order (numpy.ndarray.ravel()).

All parametrization conversions are analytical (closed-form):

• classical \(\leftrightarrow\) natural: natural_params / from_natural

• natural \(\to\) expectation: _grad_log_partition (analytical override)

• expectation \(\to\) classical: from_expectation

No Bregman solver is ever invoked. fit_mle computes \(\hat\eta = n^{-1}\sum_i t(x_i)\) and calls from_expectation (closed-form). log_prob overrides the inherited EF formula with a direct Cholesky computation for efficiency.

class normix.distributions.normal.MultivariateNormal(mu, L_Sigma)

Bases: ExponentialFamily

Multivariate Normal distribution as an exponential family.

Parameters:

• mu (jax.Array) – (d,) mean vector.

• L_Sigma (jax.Array) – (d, d) lower-triangular Cholesky factor of \(\Sigma\).

Notes

Natural parameters: \(\theta = [\Sigma^{-1}\mu,\; -\tfrac{1}{2}\operatorname{vec}(\Sigma^{-1})]\). Sufficient statistics: \(t(x) = [x,\; \operatorname{vec}(xx^\top)]\). Log-partition: \(\psi(\theta) = \tfrac{1}{2}\mu^\top\Sigma^{-1}\mu - \tfrac{1}{2}\log|\Sigma^{-1}| + \tfrac{d}{2}\log(2\pi)\). Log base measure: \(\log h(x) = 0\).

mu: Array

L_Sigma: Array

natural_params()

\(\theta = [\Sigma^{-1}\mu,\; -\tfrac{1}{2}\operatorname{vec}(\Sigma^{-1})]\)

Return type:

Array

static sufficient_statistics(x)

\(t(x) = [x,\; \operatorname{vec}(xx^\top)]\).

Parameters:

x (Array)

Return type:

Array

static log_base_measure(x)

\(\log h(x) = 0\) (base measure is Lebesgue).

Parameters:

x (Array)

Return type:

Array

classmethod from_classical(mu, sigma)

Construct from mean μ and covariance matrix Σ.

Return type:

MultivariateNormal

classmethod from_expectation(eta, **_kwargs)

Closed-form inversion \(\eta \to \theta \to (\mu, L_\Sigma)\).

\(\eta = [E[X],\; \operatorname{vec}(E[XX^\top])]\), so

\[\mu = \eta_1, \qquad \Sigma = \operatorname{reshape}(\eta_2, d, d) - \mu\mu^\top\]

Parameters:

• eta (jax.Array) – Expectation parameters of shape (d + d²,).

• **_kwargs – Ignored (accepts backend, theta0, etc. for API compatibility).

Return type:

MultivariateNormal

classmethod from_natural(theta)

Construct from natural parameters \(\theta = [\theta_1, \theta_2]\).

Recovers \(\Lambda = -2\,\mathrm{reshape}(\theta_2, d, d)\), then \(\mu = \Lambda^{-1}\theta_1\) and \(L_\Sigma = \mathrm{chol}(\Lambda^{-1})\).

Parameters:

theta (Array)

Return type:

MultivariateNormal

log_prob(x)

\(\log f(x) = -\tfrac{d}{2}\log(2\pi) - \tfrac{1}{2}\log|\Sigma| - \tfrac{1}{2}\|L_\Sigma^{-1}(x-\mu)\|^2\).

Overrides the inherited EF formula for numerical efficiency (Cholesky-direct).

Parameters:

x (Array)

Return type:

Array

mean()

\(E[X] = \mu\).

Return type:

Array

cov()

\(\mathrm{Cov}[X] = \Sigma = L_\Sigma L_\Sigma^\top\).

Return type:

Array

rvs(n, seed=42)

Draw n i.i.d. samples via JAX PRNG.

Returns:

Shape (n, d).

Return type:

jax.Array

Parameters:

• n (int)

• seed (int)

property d: int

Dimensionality.

property dim: int

Dimensionality (alias for d).

sigma()

Covariance matrix \(\Sigma = L_\Sigma L_\Sigma^\top\) (alias for cov()).

Return type:

Array

Mixture Distributions

Variance Gamma

Variance Gamma (VG) distribution.

Special case of GH with GIG → Gamma subordinator (\(b \to 0\), \(p > 0\)). \(Y \sim \mathrm{Gamma}(\alpha, \beta)\), i.e. GIG(\(p = \alpha\), \(a = 2\beta\), \(b \to 0\)).

Stored: \(\mu\), \(\gamma\), \(L_\Sigma\) (Cholesky of \(\Sigma\)), \(\alpha\) (shape), \(\beta\) (rate) of Gamma.

class normix.distributions.variance_gamma.JointVarianceGamma(mu, gamma, L_Sigma, alpha, beta)

Bases: JointNormalMixture

Joint \(f(x,y)\): \(X\mid Y \sim \mathcal{N}(\mu+\gamma y, \Sigma y)\), \(Y \sim \mathrm{Gamma}(\alpha, \beta)\).

GIG limit: \(p = \alpha\), \(a = 2\beta\), \(b \to 0\).

Parameters:

• mu (Array)

• gamma (Array)

• L_Sigma (Array)

• alpha (Array)

• beta (Array)

alpha: Array

beta: Array

subordinator()

Return the fitted subordinator distribution.

Return type:

ExponentialFamily

natural_params()

\(\theta = [\alpha-1-d/2,\; -\tfrac{1}{2}\mu^\top\Lambda\mu,\; -(\beta+\tfrac{1}{2}\gamma^\top\Lambda\gamma),\; \Lambda\gamma,\; \Lambda\mu,\; -\tfrac{1}{2}\mathrm{vec}(\Lambda)]\)

(Gamma subordinator: \(p=\alpha\), \(a=2\beta\), \(b\to 0\)).

Return type:

Array

classmethod from_classical(*, mu, gamma, sigma, alpha, beta)

classmethod from_natural(theta)

Recover classical parameters from \(\theta\).

\(\alpha = \theta_1 + 1 + d/2\), \(\beta = -\theta_3 - \gamma_{\mathrm{quad}}\).

Parameters:

theta (Array)

Return type:

JointVarianceGamma

class normix.distributions.variance_gamma.VarianceGamma(joint)

Bases: NormalMixture

Marginal Variance Gamma distribution f(x).

Parameters:

joint (JointVarianceGamma)

classmethod from_classical(*, mu, gamma, sigma, alpha, beta)

Return type:

VarianceGamma

log_prob(x)

Marginal VG log-density (own formula, no GH delegation).

\[f(x) \propto \left(\frac{q}{2c}\right)^{\nu/2} K_\nu\!\left(\sqrt{2qc}\right) \exp(\gamma^\top\Sigma^{-1}(x-\mu))\]

where \(\nu = \alpha - d/2\), \(c = \beta + \tfrac{1}{2}\gamma^\top\Lambda\gamma\), \(q = (x-\mu)^\top\Lambda(x-\mu)\).

\(q,\,2c\) are floored at LOG_EPS before forming \(\sqrt{2qc}\), matching the floor used in the \((q/2c)^{\nu/2}\) ratio — otherwise \(x=\mu\) sends the true \(q=0\) into log_kv unfloored while the ratio uses the floored value, breaking the cancellation needed for a finite limit.

Parameters:

x (Array)

Return type:

Array

property alpha: Array

\(\alpha\) — Gamma shape (forwarded from the joint).

property beta: Array

\(\beta\) — Gamma rate (forwarded from the joint).

to_generalized_hyperbolic(*, boundary_eps=0.0)

Exact embedding into the GeneralizedHyperbolic family.

Parameters:

boundary_eps (float)

class normix.distributions.variance_gamma.UnivariateVarianceGamma(joint)

Bases: _UnivariateNormalMixtureMixin, VarianceGamma

Univariate (d=1) Variance Gamma distribution.

Sibling of VarianceGamma for 1-D problems: exposes scalar mean/var/std, (n,)-shaped rvs, and cdf/ppf backed by a PINV table over the marginal log-density. EM, fit, replace, and regularisation are inherited from VarianceGamma.

classmethod from_classical(*, mu, gamma, sigma, alpha, beta)

Build from scalar or 1-D classical parameters.

mu, gamma may be scalars or (1,); sigma may be a scalar variance, (1,), or (1, 1).

Return type:

UnivariateVarianceGamma

class normix.distributions.variance_gamma.FactorVarianceGamma(mu, gamma, F, D, *, alpha, beta)

Bases: FactorNormalMixture

Factor-analysis Variance Gamma: \(Y \sim \mathrm{Gamma}(\alpha, \beta)\), \(\Sigma = F F^\top + \mathrm{diag}(D)\).

GIG limit of the subordinator: \(p = \alpha\), \(a = 2\beta\), \(b \to 0\).

Parameters:

• mu (Array)

• gamma (Array)

• F (Array)

• D (Array)

classmethod from_classical(*, mu, gamma, F, D, alpha, beta)

Return type:

FactorVarianceGamma

property alpha: Array

property beta: Array

log_prob(x)

Marginal VG log-density evaluated with Woodbury Σ-solve.

Parameters:

x (Array)

Return type:

Array

Normal Inverse Gamma

Normal-Inverse Gamma (NInvG) distribution.

Special case of GH with GIG → InverseGamma subordinator (\(a \to 0\), \(p < 0\)). \(Y \sim \mathrm{InvGamma}(\alpha, \beta)\), i.e. GIG(\(p = -\alpha\), \(a \to 0\), \(b = 2\beta\)).

Stored: \(\mu\), \(\gamma\), \(L_\Sigma\) (Cholesky of \(\Sigma\)), \(\alpha\) (shape), \(\beta\) (rate) of InverseGamma.

class normix.distributions.normal_inverse_gamma.JointNormalInverseGamma(mu, gamma, L_Sigma, alpha, beta)

Bases: JointNormalMixture

Joint \(f(x,y)\): \(X\mid Y \sim \mathcal{N}(\mu+\gamma y, \Sigma y)\), \(Y \sim \mathrm{InvGamma}(\alpha, \beta)\).

GIG limit: \(p = -\alpha\), \(a \to 0\), \(b = 2\beta\).

The \(\gamma y\) mean term is the skewness that generalizes this beyond the textbook Normal-inverse-gamma conjugate prior — see NormalInverseGamma for its role in the marginal density.

Parameters:

• mu (Array)

• gamma (Array)

• L_Sigma (Array)

• alpha (Array)

• beta (Array)

alpha: Array

beta: Array

subordinator()

Return the fitted subordinator distribution.

Return type:

ExponentialFamily

natural_params()

\(\theta = [-(\alpha+1)-d/2,\; -(\beta+\tfrac{1}{2}\mu^\top\Lambda\mu),\; -\tfrac{1}{2}\gamma^\top\Lambda\gamma,\; \Lambda\gamma,\; \Lambda\mu,\; -\tfrac{1}{2}\mathrm{vec}(\Lambda)]\)

(InverseGamma subordinator: \(p=-\alpha\), \(a\to 0\), \(b=2\beta\)).

Return type:

Array

classmethod from_classical(*, mu, gamma, sigma, alpha, beta)

classmethod from_natural(theta)

Recover classical parameters from \(\theta\).

\(\alpha = -(\theta_1 + d/2) - 1\), \(\beta = -\theta_2 - \mu_{\mathrm{quad}}\).

Parameters:

theta (Array)

Return type:

JointNormalInverseGamma

class normix.distributions.normal_inverse_gamma.NormalInverseGamma(joint)

Bases: NormalMixture

Marginal Normal-Inverse Gamma distribution \(f(x)\).

Not the textbook Normal-inverse-gamma conjugate prior — that distribution is a pure variance mixture, \(X\mid Y\sim\mathcal N(\mu,\Sigma Y)\), whose marginal is a plain Student-t with no special functions. This class is the GH-family variance-mean mixture (see JointNormalInverseGamma): \(X\mid Y\sim\mathcal N(\mu+\gamma Y,\Sigma Y)\). The extra \(\gamma Y\) term is why log_prob() needs log_kv(): expanding the quadratic form adds a term linear in \(Y\), turning the \(\int_0^\infty dY\) integral from a Gamma integral into a Generalized Inverse Gaussian one, \(2(b/a)^{p/2}K_p(\sqrt{ab})\) with \(a=\gamma^\top\Sigma^{-1}\gamma\). At \(\gamma=0\) this reduces exactly to the textbook result \(f(x)=t\big(x\mid\nu=2\alpha,\ \hat\mu=\mu,\ \hat\Sigma=\tfrac{\beta}{\alpha}\Sigma\big)\) — see The Generalized Hyperbolic Distribution “Special Cases”.

Parameters:

joint (JointNormalInverseGamma)

classmethod from_classical(*, mu, gamma, sigma, alpha, beta)

log_prob(x)

Marginal NInvG log-density (own formula, no GH delegation).

GIG params: \(p=-\alpha\), \(a=\gamma^\top\Lambda\gamma\), \(b=2\beta+Q(x)\). The normalising integral is \(2(b/a)^{p/2} K_p(\sqrt{ab})\).

\(a,b\) are floored at LOG_EPS before forming \(\sqrt{ab}\), matching the floor used in the \((b/a)^{p/2}\) ratio. Without this, \(\gamma=0\) (the common no-skew case) sends the true \(a=0\) into log_kv unfloored while the ratio uses the floored value, breaking the cancellation between the two that is needed for a finite limit.

Parameters:

x (Array)

Return type:

Array

property alpha: Array

\(\alpha\) — InverseGamma shape (forwarded from the joint).

property beta: Array

\(\beta\) — InverseGamma rate (forwarded from the joint).

to_generalized_hyperbolic(*, boundary_eps=0.0)

Exact embedding into the GeneralizedHyperbolic family.

Parameters:

boundary_eps (float)

class normix.distributions.normal_inverse_gamma.UnivariateNormalInverseGamma(joint)

Bases: _UnivariateNormalMixtureMixin, NormalInverseGamma

Univariate (d=1) Normal-Inverse-Gamma distribution.

Sibling of NormalInverseGamma for 1-D problems; see UnivariateVarianceGamma for the contract.

classmethod from_classical(*, mu, gamma, sigma, alpha, beta)

Return type:

UnivariateNormalInverseGamma

class normix.distributions.normal_inverse_gamma.FactorNormalInverseGamma(mu, gamma, F, D, *, alpha, beta)

Bases: FactorNormalMixture

Factor-analysis Normal-Inverse-Gamma: \(Y \sim \mathrm{InvGamma}(\alpha, \beta)\), \(\Sigma = F F^\top + \mathrm{diag}(D)\).

GIG limit: \(p = -\alpha\), \(a \to 0\), \(b = 2\beta\).

Parameters:

• mu (Array)

• gamma (Array)

• F (Array)

• D (Array)

classmethod from_classical(*, mu, gamma, F, D, alpha, beta)

Return type:

FactorNormalInverseGamma

property alpha: Array

property beta: Array

log_prob(x)

Marginal \(\log f(x)\) for a single observation.

Parameters:

x (Array)

Return type:

Array

Normal Inverse Gaussian

Normal-Inverse Gaussian (NIG) distribution.

Special case of GH with GIG → InverseGaussian subordinator (\(p = -1/2\)). \(Y \sim \mathrm{InvGaussian}(\mu_{IG}, \lambda)\), i.e. GIG(\(p = -1/2\), \(a = \lambda/\mu_{IG}^2\), \(b = \lambda\)).

Stored: \(\mu\), \(\gamma\), \(L_\Sigma\) (Cholesky of \(\Sigma\)), \(\mu_{IG}\) (IG mean), \(\lambda\) (IG shape).

class normix.distributions.normal_inverse_gaussian.JointNormalInverseGaussian(mu, gamma, L_Sigma, mu_ig, lam)

Bases: JointNormalMixture

Joint \(f(x,y)\): \(X\mid Y \sim \mathcal{N}(\mu+\gamma y, \Sigma y)\), \(Y \sim \mathrm{InvGaussian}(\mu_{IG}, \lambda)\).

Stored: \(\mu_{IG}\) (IG mean) and \(\lambda\) (IG shape) directly. GIG params: \(p = -1/2\), \(a = \lambda/\mu_{IG}^2\), \(b = \lambda\).

Parameters:

• mu (Array)

• gamma (Array)

• L_Sigma (Array)

• mu_ig (Array)

• lam (Array)

mu_ig: Array

lam: Array

subordinator()

Return the fitted subordinator distribution.

Return type:

ExponentialFamily

natural_params()

\(\theta = [-3/2-d/2,\; -(\lambda/2+\tfrac{1}{2}\mu^\top\Lambda\mu),\; -(\lambda/(2\mu_{IG}^2)+\tfrac{1}{2}\gamma^\top\Lambda\gamma),\; \Lambda\gamma,\; \Lambda\mu,\; -\tfrac{1}{2}\mathrm{vec}(\Lambda)]\)

where \(p=-1/2\), \(a=\lambda/\mu_{IG}^2\), \(b=\lambda\), aligned with GIG natural parameters on \([\log y,\,1/y,\,y]\).

Return type:

Array

classmethod from_classical(*, mu, gamma, sigma, mu_ig, lam)

classmethod from_natural(theta)

Recover classical parameters from \(\theta\).

From \(b = -2\theta_2 - 2\mu_{\mathrm{quad}} = \lambda\) and \(a = -2\theta_3 - 2\gamma_{\mathrm{quad}} = \lambda/\mu_{IG}^2\), so \(\mu_{IG} = \sqrt{b/a}\).

Parameters:

theta (Array)

Return type:

JointNormalInverseGaussian

log_density_power(alpha)

Log density-power integral via the exact JointGeneralizedHyperbolic embedding.

Unlike the VG / NInvG joints (whose shape \(p\) is a free coordinate read back from \(\theta_1\)), this class fixes \(p=-1/2\) and drops \(\theta_1\) from _log_partition_from_theta() (using \(K_{-1/2}\)). The density-power scaling \(\theta\mapsto\alpha\theta\) moves \(\theta_1\) (i.e. the GIG order) off \(-1/2\), so the constant-base-measure formula in ExponentialFamily must be evaluated on the GH embedding, where \(p\) is carried faithfully through log_kv. This drives entropy(), varentropy(), and renyi().

Parameters:

alpha (Array)

Return type:

Array

class normix.distributions.normal_inverse_gaussian.NormalInverseGaussian(joint)

Bases: NormalMixture

Marginal Normal-Inverse Gaussian distribution f(x).

Parameters:

joint (JointNormalInverseGaussian)

classmethod from_classical(*, mu, gamma, sigma, mu_ig, lam)

log_prob(x)

Marginal NIG log-density.

Uses \(K_{-1/2}(z) = \sqrt{\pi/(2z)}\,e^{-z}\) for the normalisation, leaving only one \(\log K_\nu\) call at order \(\nu = -1/2 - d/2\).

Parameters:

x (Array)

Return type:

Array

property mu_ig: Array

\(\mu_{IG}\) — InverseGaussian mean (forwarded from the joint).

property lam: Array

\(\lambda\) — InverseGaussian shape (forwarded from the joint).

regularize_a_eq_b()

Rescale so \(a = b = \sqrt{ab}\).

For NIG, \(a = \lambda/\mu_{IG}^2,\;b = \lambda\), so \(s = \sqrt{a/b} = 1/\mu_{IG}\). After rescaling \(\mu_{IG}' = 1\), i.e. the InverseGaussian has unit mean.

Return type:

NormalInverseGaussian

to_generalized_hyperbolic()

Exact embedding into the GeneralizedHyperbolic family.

No boundary approximation: NIG sits in the strict interior of GH’s parameter space (\(p = -1/2,\; a = \lambda/\mu_{IG}^2,\; b = \lambda\)).

class normix.distributions.normal_inverse_gaussian.UnivariateNormalInverseGaussian(joint)

Bases: _UnivariateNormalMixtureMixin, NormalInverseGaussian

Univariate (d=1) Normal-Inverse-Gaussian distribution.

Sibling of NormalInverseGaussian for 1-D problems; see UnivariateVarianceGamma for the contract.

classmethod from_classical(*, mu, gamma, sigma, mu_ig, lam)

Return type:

UnivariateNormalInverseGaussian

class normix.distributions.normal_inverse_gaussian.FactorNormalInverseGaussian(mu, gamma, F, D, *, mu_ig, lam)

Bases: FactorNormalMixture

Factor-analysis Normal-Inverse-Gaussian: \(Y \sim \mathrm{InvGaussian}(\mu_{IG}, \lambda)\), \(\Sigma = F F^\top + \mathrm{diag}(D)\).

GIG params: \(p = -1/2\), \(a = \lambda/\mu_{IG}^2\), \(b = \lambda\).

Parameters:

• mu (Array)

• gamma (Array)

• F (Array)

• D (Array)

classmethod from_classical(*, mu, gamma, F, D, mu_ig, lam)

Return type:

FactorNormalInverseGaussian

property mu_ig: Array

property lam: Array

log_prob(x)

Marginal \(\log f(x)\) for a single observation.

Parameters:

x (Array)

Return type:

Array

regularize_a_eq_b()

Rescale so \(a = b\). For NIG (\(a = \lambda/\mu_{IG}^2\), \(b = \lambda\)), this means \(\mu_{IG} \to 1\).

Return type:

FactorNormalInverseGaussian

Generalized Hyperbolic

Generalized Hyperbolic (GH) distribution.

Joint: \(X \mid Y \sim \mathcal{N}(\mu + \gamma y, \Sigma y)\), \(Y \sim \mathrm{GIG}(p, a, b)\).

Marginal: \(\mathrm{GH}(\mu, \gamma, \Sigma, p, a, b)\).

Marginal log-density (closed form via Bessel functions):

Let \(Q(x) = (x-\mu)^\top \Sigma^{-1}(x-\mu)\), \(A = a + \gamma^\top \Sigma^{-1} \gamma\).

\[\log f(x) = -\tfrac{d}{2}\log(2\pi) - \tfrac{1}{2}\log|\Sigma| + \tfrac{p}{2}(\log a - \log b) - \log K_p(\sqrt{ab}) + \tfrac{d/2-p}{2}\log\frac{A}{Q(x)+b} + \log K_{p-d/2}\!\left(\sqrt{A(Q(x)+b)}\right) + \gamma^\top \Sigma^{-1}(x - \mu)\]

Posterior: \(Y \mid X = x \sim \mathrm{GIG}(p - d/2,\; a + \gamma^\top\Sigma^{-1}\gamma,\; b + (x-\mu)^\top\Sigma^{-1}(x-\mu))\).

class normix.distributions.generalized_hyperbolic.JointGeneralizedHyperbolic(mu, gamma, L_Sigma, p, a, b)

Bases: JointNormalMixture

Joint \(f(x,y)\): \(X\mid Y \sim \mathcal{N}(\mu+\gamma y, \Sigma y)\), \(Y \sim \mathrm{GIG}(p, a, b)\).

Stored: mu, gamma, L_Sigma (from JointNormalMixture) + p, a, b (GIG parameters).

Parameters:

• mu (Array)

• gamma (Array)

• L_Sigma (Array)

• p (Array)

• a (Array)

• b (Array)

p: Array

a: Array

b: Array

subordinator()

Return the fitted subordinator distribution.

Return type:

ExponentialFamily

to_joint_generalized_hyperbolic(*, boundary_eps=0.0)

Identity: JointGH is already the divergence gauge.

Parameters:

boundary_eps (float)

natural_params()

\(\theta = [p-1-d/2,\; -(b/2+\tfrac{1}{2}\mu^\top\Sigma^{-1}\mu),\; -(a/2+\tfrac{1}{2}\gamma^\top\Sigma^{-1}\gamma),\; \Sigma^{-1}\gamma,\; \Sigma^{-1}\mu,\; -\tfrac{1}{2}\mathrm{vec}(\Sigma^{-1})]\)

The scalar coefficients on sufficient statistics \(1/y\) and \(y\) match the GIG convention \(\theta_{\mathrm{GIG}} = [p-1,\,-b/2,\,-a/2]\) on \(t_Y = [\log y,\,1/y,\,y]\).

Return type:

Array

classmethod from_classical(*, mu, gamma, sigma, p, a, b)

Construct from classical parameters.

classmethod from_natural(theta)

Recover classical parameters from \(\theta\).

\(p = \theta_1 + 1 + d/2\), \(b = -2\theta_2 - 2\mu_{\mathrm{quad}}\), \(a = -2\theta_3 - 2\gamma_{\mathrm{quad}}\).

Parameters:

theta (Array)

Return type:

JointGeneralizedHyperbolic

to_joint_variance_gamma()

KL projection onto JointVarianceGamma (Gamma subordinator).

Return type:

JointVarianceGamma

to_joint_normal_inverse_gamma()

KL projection onto JointNormalInverseGamma.

Return type:

JointNormalInverseGamma

to_joint_normal_inverse_gaussian()

KL projection onto JointNormalInverseGaussian.

Return type:

JointNormalInverseGaussian

class normix.distributions.generalized_hyperbolic.GeneralizedHyperbolic(joint)

Bases: NormalMixture

Marginal Generalized Hyperbolic distribution \(f(x)\).

Stores a JointGeneralizedHyperbolic. Provides:

• log_prob(x) — closed-form Bessel expression

• e_step, m_step — for EM fitting

• fit(X, ...) — convenience fitting method

Parameters:

joint (JointGeneralizedHyperbolic)

classmethod from_classical(*, mu, gamma, sigma, p, a, b)

Return type:

GeneralizedHyperbolic

log_prob(x)

Marginal \(\log f(x)\).

\[f(x) \propto \left(\frac{A}{Q(x)+b}\right)^{(d/2-p)/2} K_{p-d/2}\!\left(\sqrt{A(Q(x)+b)}\right) \exp\!\left(\gamma^\top\Sigma^{-1}(x-\mu)\right)\]

where \(Q(x) = (x-\mu)^\top\Sigma^{-1}(x-\mu)\), \(A = a + \gamma^\top\Sigma^{-1}\gamma\).

Parameters:

x (Array)

Return type:

Array

property p: Array

\(p\) — GIG order (forwarded from the joint).

property a: Array

\(a\) — GIG concentration (forwarded from the joint).

property b: Array

\(b\) — GIG concentration (forwarded from the joint).

m_step(eta, **kwargs)

Full M-step with warm-started + sanity-checked subordinator solve.

Overrides the base cls.from_expectation(eta) because the GIG solver benefits substantially from warm-starting at the current \(\theta\) and from a fall-back when the solver wanders into unsane regions; both are managed by m_step_subordinator().

Return type:

GeneralizedHyperbolic

m_step_subordinator(eta, **kwargs)

M-step for the subordinator only (MCECM Cycle 2).

Reads the subordinator-relevant fields of eta; normal parameters are read from self._joint and copied unchanged. Subclasses with iterative solvers may override to add warm-start or sanity-check fallbacks.

Return type:

GeneralizedHyperbolic

regularize_a_eq_b()

Rescale so \(a = b = \sqrt{ab}\) (orbit invariant).

Picks \(s = \sqrt{a/b}\). Idempotent: applying twice leaves the model unchanged.

Return type:

GeneralizedHyperbolic

to_variance_gamma()

KL projection onto the VarianceGamma family.

Return type:

VarianceGamma

to_normal_inverse_gamma()

KL projection onto the NormalInverseGamma family.

Return type:

NormalInverseGamma

to_normal_inverse_gaussian()

KL projection onto the NormalInverseGaussian family.

Return type:

NormalInverseGaussian

classmethod default_init(X)

Warm-start from the best of NIG / VG / NInvG sub-model fits.

Runs 5 EM iterations (JAX backend, Newton method) for each special case, converts each to GH parametrisation, and selects the candidate with the highest marginal log-likelihood. A moment-based fallback (\(p=1, a=1, b=1\)) is included as a fourth candidate.

Fully JAX-native: no try/except, no Python branching on data values.

Parameters:

X (Array)

Return type:

GeneralizedHyperbolic

class normix.distributions.generalized_hyperbolic.UnivariateGeneralizedHyperbolic(joint)

Bases: _UnivariateNormalMixtureMixin, GeneralizedHyperbolic

Univariate (d=1) Generalized Hyperbolic distribution.

Sibling of GeneralizedHyperbolic for 1-D problems; see UnivariateVarianceGamma for the contract.

classmethod from_classical(*, mu, gamma, sigma, p, a, b)

Return type:

UnivariateGeneralizedHyperbolic

class normix.distributions.generalized_hyperbolic.FactorGeneralizedHyperbolic(mu, gamma, F, D, *, p, a, b)

Bases: FactorNormalMixture

Factor-analysis Generalized Hyperbolic: \(Y \sim \mathrm{GIG}(p, a, b)\), \(\Sigma = F F^\top + \mathrm{diag}(D)\).

The GIG subordinator’s M-step uses warm-started numerical optimisation with the same sanity-check fall-back as the standard GeneralizedHyperbolic.

Parameters:

• mu (Array)

• gamma (Array)

• F (Array)

• D (Array)

classmethod from_classical(*, mu, gamma, F, D, p, a, b)

Return type:

FactorGeneralizedHyperbolic

property p: Array

property a: Array

property b: Array

log_prob(x)

Marginal \(\log f(x)\) for a single observation.

Parameters:

x (Array)

Return type:

Array

m_step(eta, **kwargs)

Full M-step with warm-started, sanity-checked GIG solve.

Mirrors GeneralizedHyperbolic.m_step(): do the closed-form factor M-step first, then run the subordinator update with warm-start and fall-back.

Return type:

FactorGeneralizedHyperbolic

m_step_subordinator(eta, **kwargs)

M-step for the subordinator only (MCECM cycle 2).

Reads the subordinator-relevant fields of eta; \((\mu, \gamma, F, D)\) are read from self and copied unchanged. Subclasses with iterative solvers may override to add a warm-start or sanity-check fallback.

Return type:

FactorGeneralizedHyperbolic

regularize_a_eq_b()

Rescale so \(a = b = \sqrt{ab}\) (orbit invariant).

Mirrors GeneralizedHyperbolic.regularize_a_eq_b(). Picks \(s = \sqrt{a/b}\).

Return type:

FactorGeneralizedHyperbolic

classmethod default_init(X, *, r=1)

Warm-start from the best of FactorNIG / FactorVG / FactorNInvG fits.

Mirrors GeneralizedHyperbolic.default_init(): runs a few EM iterations of each special-case factor family, converts each subordinator to the GIG embedding, and selects the candidate with the highest marginal log-likelihood. A moment-based fallback (\(p = 1, a = 1, b = 1\)) is included as a fourth candidate.

Fully JAX-native: no try/except, no Python branching on data values.

Parameters:

• X (Array)

• r (int)

Return type:

FactorGeneralizedHyperbolic