EM for generalized hyperbolic distributions

Expectation–maximization fits a latent-variable model by alternating posterior expectations and complete-data likelihood maximization. The general method is due to Dempster, Laird, and Rubin, 1977, Maximum likelihood from incomplete data via the EM algorithm. This GH specialization follows the upstream derivation, which credits Hu (2005), Calibration of multivariate generalized hyperbolic distributions using the EM algorithm.

Read ordinary EM for the likelihood inequality, and GH for the joint density and parameter conventions. We treat independent observations $x_1,\ldots,x_n\in\mathbb R^d$, with positive-definite $\Sigma$ and interior GIG parameters $a,b>0$, unless a boundary case is stated.

Complete and observed data¶

Write the classical parameter tuple as $\vartheta=(\mu,\gamma,\Sigma,p,a,b)$; reserve $\theta$ for the joint family’s natural coordinates. The model is

For $d=1$, the conditioning note uses $W=Y$, $\beta=\gamma$, and $\sigma^2=\Sigma$. The observed likelihood uses $f_X(x;\vartheta)=\int f(x,y;\vartheta)\,dy$. The joint GH law is an exponential family, which makes its expected log likelihood easier to optimize.

Conditional distribution of the mixing variable¶

Expanding the normal quadratic in the joint density and collecting powers of $y$ gives

Set $r=p-d/2$, $A=a+\gamma^\top\Sigma^{-1}\gamma$, and $B(x)=b+(x-\mu)^\top\Sigma^{-1}(x-\mu)$. All are evaluated at the current iterate when used in an E-step. Then

Interior assumptions guarantee $A,B(x)>0$ for every observation. The GIG normalizing integral therefore supplies every real posterior power moment:

Differentiating at $s=0$ yields

The three expectations of $Y$, $Y^{-1}$, and $\log Y$ must be evaluated separately. Replacing $Y$ by its posterior mean inside nonlinear statistics would give a different algorithm.

E-step in one statistic order¶

Use the joint sufficient-statistic blocks $t(x,y)=(\log y,y^{-1},y,x,x/y,xx^\top/y)$, as in the GIG and GH notes. At iterate $\vartheta_k$ let

The E-step averages all six blocks:

The source EM note uses the first three slots in a different order. Here $\widehat\eta_1$ always means the log moment, $\widehat\eta_2$ the inverse moment, and $\widehat\eta_3$ the first moment. This order also governs the mixture architecture.

M-step as an expectation-to-parameter map¶

Holding the old posterior fixed, maximize

In natural coordinates this equals $\langle\theta,\widehat\eta^{(k)}\rangle-\psi(\theta)$ plus a constant. For an attained interior maximum of the full joint family, the new model matches the posterior sufficient-statistic averages. Suppressing $(k)$ on $\widehat\eta$, the normal-block update is

The two vector updates solve $\widehat\eta_4=\mu+\gamma\widehat\eta_3$ and $\widehat\eta_5=\mu\widehat\eta_2+\gamma$. The covariance formula then follows from the sixth joint moment. The last line is the GIG likelihood problem; it usually requires numerical optimization. Only the normal block is closed form in the general GH case.

These formulas require a nonzero denominator, a positive-definite recovered covariance, and a feasible attained mixing-law optimum. Positive-definiteness can fail with insufficiently varied data. Convexity of the GIG subproblem does not remove its domain constraints or numerical conditioning problems.

With exact posteriors and an M-step that increases $Q$, the EM likelihood inequality ensures

This concerns the observed likelihood. It does not guarantee a global MLE, a unique stationary point, convergence of the parameter sequence, or a particular speed. Approximate numerical steps should be checked against their actual objective.

Scale normalization and equivariance¶

The GH scale action is

It represents the same observed law by replacing $Y$ with $Y'=cY$. The posterior parameters become $(r,A/c,cB(x))$. Consequently,

The E-step output transforms as

Substitution in (9), including the rescaled GIG optimum, proves $F(T_c\vartheta)=T_cF(\vartheta)$ for the exact unconstrained EM map $F$ when the maximizers are unique (or selected consistently under scaling).

To impose $|\Sigma|=1$ after an update, choose $c=|\Sigma|^{1/d}$ and transform all coupled parameters:

The observed density and its likelihood are unchanged. This is a choice of scale representative, not a penalty or a guarantee of numerical stability. Fixed bounds, shrinkage, approximate solvers, or inconsistent choices among maximizers can break the exact equivariance argument.

ECM and multiple cycles¶

Meng and Rubin’s 1993 ECM framework replaces one M-step by conditional maximizations of parameter blocks. The upstream note discusses its multi-cycle version (MCECM) for GH, also citing McNeil, Frey, and Embrechts, Quantitative Risk Management (2010). A two-cycle construction in the unconstrained redundant coordinates is:

1. Compute the posterior moments at $\vartheta_k$, then maximize the normal block $(\mu,\gamma,\Sigma)$ with $(p,a,b)$ fixed.

2. Recompute the posterior moments at this intermediate parameter and maximize the GIG block with the new normal block fixed.

Each conditional maximization must increase the corresponding current $Q$. A full coupled scale normalization can then be applied without changing the observed law. If instead one imposes a determinant constraint inside a conditional maximization while holding the mixing parameters fixed, that is a constrained optimization problem: simply normalizing the unconstrained covariance is not a justification for its ascent property.

EM and valid MCECM updates share the likelihood-ascent argument. Neither has a general guarantee of reaching the global MLE, and their numerical costs and convergence rates depend on the problem.

Special mixing families¶

Write $l=\mathbb E[\log Y]$, $u=\mathbb E[Y^{-1}]$, and $v=\mathbb E[Y]$ for the posterior averages at the relevant E-step. Denote the digamma function by $\psi_0$ and its derivative by $\psi_1$, to distinguish them from the log-partition $\psi$.

Gamma mixing: variance gamma¶

For $Y\sim\operatorname{Gamma}(\alpha,\rho)$ with rate $\rho$, $\mathbb E[\log Y]=\psi_0(\alpha)-\log\rho$ and $\mathbb E[Y]=\alpha/\rho$. The mixing-law update solves

A Newton proposal is

Safeguards must keep $\alpha>0$ and verify improvement. The posterior is GIG with $r=\alpha-d/2$, $A=2\rho+\gamma^\top\Sigma^{-1}\gamma$, and $B=q(x)$. At $x=\mu$, $B=0$; posterior propriety and required inverse moments then need separate checking. For example, this gamma posterior has finite inverse moment only if $r>1$. The interior GIG formula does not cover that point without a boundary analysis.

Inverse-Gaussian mixing: NIG¶

For an inverse Gaussian with mean $m$ and shape $\lambda$, $\mathbb E[Y]=m$ and $\mathbb E[Y^{-1}]=1/m+1/\lambda$. The closed-form mixing update is

It requires $uv>1$; equality is a constant-mixing limit. Its GIG embedding has $p=-1/2$, $a=\lambda/m^2$, $b=\lambda$. The posterior is generally GIG with order $-1/2-d/2$, rather than another inverse Gaussian.

Inverse-gamma mixing¶

For $Y\sim\operatorname{InvGamma}(\alpha,\rho)$, $\mathbb E[\log Y]=\log\rho-\psi_0(\alpha)$ and $\mathbb E[Y^{-1}]=\alpha/\rho$. Hence

with Newton proposal

This is again a scalar solve, not a general closed-form shape estimate. For $\gamma=0$ the posterior has $A=0$ and is inverse gamma with shape $\alpha+d/2$; a finite posterior first moment requires $\alpha+d/2>1$ when that moment is needed. Symmetry-constrained updates can require fewer statistics than the full variance–mean model.

The normal block still uses inverse and first posterior moments when both location and mean coupling are free. Reducing the mixing-law statistics does not automatically eliminate these normal-block requirements.

Numerical considerations and implementation¶

Large $|p-d/2|$ or extreme $\sqrt{AB(x)}$ can overflow or underflow direct Bessel evaluations. Ratios, log-space evaluation, and stable order derivatives matter. The covariance update can suffer cancellation or rank deficiency; determinant normalization does not repair its condition number. Penalties change the optimization objective and need their own ascent analysis.

See why not gradient descent for the optimization comparison and exponential-family core for the expectation-to-natural map. Runnable fitters remain in the upstream EM guide and package API. Continue with online EM, penalized shrinkage, and the mathematical EM framework. The full EM design and Bessel/solver design remain upstream implementation references.

Source and adaptation¶

Adapted from xshi19/normix, docs/theory/em_algorithm.md, at revision 763bb3608920661a012cf089888d349fbf680aad (2026-09-13 import). Copyright (c) 2020 xshi19. Licensed under MIT. The pinned source records the original version. Notation, mathematical qualifications, and links were adapted for this site; API recipes, executable package cells, and notebook plots are omitted. No upstream benchmark execution or formal-proof verification is claimed for this adaptation.