Gauss-Hermite Quadrature for Gaussian Random Effects

Why use Gauss–Hermite quadrature?

A marginal likelihood in a nonlinear mixed-effects model averages each subject’s conditional likelihood over that subject’s random-effects distribution. Even when the distribution is Gaussian, the conditional likelihood is generally a nonlinear function of the random effects, so the integral is not usually available in closed form. Gauss–Hermite quadrature replaces that integral with a finite weighted sum. The Gaussian density is built into the quadrature rule, concentrating the calculation where the random-effects distribution places its mass.

The central practical point is that Gauss–Hermite nodes are not automatically values of a standard normal random variable. The nodes and weights depend on the convention used for the Hermite polynomials. This tutorial uses the physicists’ convention, based on the weight $e^{-x^2}$, and shows explicitly how it maps to a Gaussian random-effects integral.

From a Gaussian integral to a weighted sum

Definition: The $m$-point physicists’ Gauss–Hermite rule uses nodes $x_1,\ldots,x_m$ and weights $w_1,\ldots,w_m$ to approximate $$ \int_{-\infty}^{\infty} e^{-x^2} f(x)\,dx \approx \sum_{j=1}^{m} w_j f(x_j). $$ The nodes are the roots of the physicists’ Hermite polynomial $H_m(x)$, and the weights are chosen for this convention. In particular, the rule integrates polynomials of degree at most $2m-1$ exactly.

For a scalar random effect $b\sim N(\mu,\sigma^2)$, write $b=\mu+\sigma z$, where $z\sim N(0,1)$. Then set $z=\sqrt{2}x$. Since $dz=\sqrt{2}\,dx$, the standard-normal density and the change of variable give

\mathbb{E}[f(b)] =\frac{1}{\sqrt{2\pi}}\int_{-\infty}^{\infty} f(\mu+\sigma z)e^{-z^2/2}\,dz =\frac{1}{\sqrt{\pi}}\int_{-\infty}^{\infty} e^{-x^2}f(\mu+\sqrt{2}\sigma x)\,dx.

Applying the quadrature rule yields

\mathbb{E}[f(b)] \approx \frac{1}{\sqrt{\pi}}\sum_{j=1}^{m} w_j f(\mu+\sqrt{2}\sigma x_j).

Equivalently, define the standard-normal nodes $z_j=\sqrt{2}x_j$ and normalized weights $\alpha_j=w_j/\sqrt{\pi}$. Then $\sum_j\alpha_j=1$, and the approximation is $\sum_j\alpha_j f(\mu+\sigma z_j)$. These normalized weights make clear that the procedure approximates an expectation, not an unnormalized integral.

Why the rule is exact for low-degree polynomials

The polynomial exactness property explains an important feature of the method: increasing the number of nodes raises the degree of polynomials integrated exactly. For the physicists’ Hermite polynomial $H_m$, the nodes $x_j$ are its roots, and $H_m$ is orthogonal to every polynomial of degree less than $m$ under the weight $e^{-x^2}$.

To see the exactness, take any polynomial $p$ of degree at most $2m-1$. Polynomial division writes $p=qH_m+r$, where both $q$ and $r$ have degree at most $m-1$. Orthogonality gives $\int e^{-x^2}q(x)H_m(x)\,dx=0$. At every quadrature node, $H_m(x_j)=0$, so the quadrature sum for $qH_m$ is also zero. The remaining polynomial $r$ has degree at most $m-1$; the weights are interpolatory and integrate it exactly. Thus the rule integrates $p=qH_m+r$ exactly. For a general nonlinear likelihood, exactness does not follow, and the weighted sum is an approximation whose adequacy must be checked.

Quadrature for a marginal likelihood

Let $y_i$ denote the observations for subject $i$, let $b_i$ be that subject’s random-effects vector, and let $\theta$ denote the population and residual parameters in the conditional model. Write the conditional likelihood as $L_i(b_i;\theta)=p(y_i\mid b_i;\theta)$. If $b_i\sim N(\mu,\Omega)$, the subject-specific marginal likelihood contribution is

L_i^{\mathrm{marg}}(\theta) =\int_{\mathbb{R}^q} L_i(b_i;\theta)\, \phi_q(b_i;\mu,\Omega)\,db_i.

Here $\phi_q$ is the $q$-dimensional Gaussian density. When $\Omega$ is positive definite, choose a matrix $C$ satisfying $CC^{\mathsf T}=\Omega$, such as its Cholesky factor, and set $b_i=\mu+Cz$, where the components of $z$ are independent standard normal variables. The Gaussian density and Jacobian then reduce the integral to an expectation under $z\sim N(0,I_q)$. A product of one-dimensional Gauss–Hermite rules gives

L_i^{\mathrm{marg}}(\theta)\approx \sum_{j_1=1}^{m}\cdots\sum_{j_q=1}^{m} \left(\prod_{k=1}^{q}\alpha_{j_k}\right) L_i\!\left(\mu+C \begin{bmatrix} z_{j_1}\\ \vdots\\ z_{j_q} \end{bmatrix};\theta\right).

This formula uses the same one-dimensional normalized weights $\alpha_j$ and standard-normal nodes $z_j$ in each coordinate. The subject’s marginal log-likelihood contribution is the logarithm of the resulting weighted sum. It is not the weighted average of the conditional log likelihoods: in general, $\log\{\sum_j\alpha_j L_i(b_j;\theta)\}\ne\sum_j\alpha_j\log L_i(b_j;\theta)$. The latter averages log likelihoods rather than marginalizing likelihoods.

Worked Example: A Two-Node Expectation With a Shifted Gaussian

Suppose $b\sim N(0.2,0.25)$, so its mean is $0.2$ and its standard deviation is $0.5$. Use two-point Gauss–Hermite quadrature to approximate $\mathbb{E}[e^b]$. For this rule, the standard-normal nodes are $-1,1$, each with normalized weight $1/2$.

The corresponding random-effect values are $0.2+0.5(-1)=-0.3$ and $0.2+0.5(1)=0.7$. Therefore,

\mathbb{E}[e^b]\approx \frac{1}{2}e^{-0.3}+\frac{1}{2}e^{0.7} =\frac{1}{2}(0.7408182207)+\frac{1}{2}(2.013752707) =1.377285464\approx 1.37729.

The exact Gaussian expectation is $\exp(\mu+\sigma^2/2)=\exp(0.325)\approx1.38403$. Thus, in this example, the two-node rule underestimates the expectation by about $0.00675$. The exponential is not a polynomial of degree at most three, so the two-node exactness result does not guarantee an exact answer.

Worked Example: A Three-Node Marginal Likelihood Integral

Consider a scalar random effect $\eta\sim N(0,1)$ and a hypothetical conditional likelihood contribution $L(\eta)=e^{-\eta^2/2}$. The marginal likelihood contribution is $\mathbb{E}[L(\eta)]$. For three-point Gauss–Hermite quadrature, the standard-normal nodes are $-\sqrt{3},0,\sqrt{3}$, with normalized weights $1/6,2/3,1/6$.

Substituting these nodes and weights gives

L^{\mathrm{marg}}\approx \frac{1}{6}e^{-3/2}+\frac{2}{3}e^0+\frac{1}{6}e^{-3/2} =\frac{2}{3}+\frac{1}{3}e^{-3/2} =0.741043387\approx0.7410.

The exact integral is $1/\sqrt{2}\approx0.7071$, because multiplying the standard-normal density by $e^{-\eta^2/2}$ gives an integrand proportional to $e^{-\eta^2}$. The three-node value is therefore an approximation, not an exact marginal likelihood. In a real calculation, the relevant question is whether this error is acceptably small for the inference being made.

Worked Example: Correlated Random Effects and a Conditional Density

For a hypothetical subject, suppose $b=(b_1,b_2)^{\mathsf T}\sim N(0,\Omega)$, where $\Omega=\begin{bmatrix}1&0.5\\0.5&1\end{bmatrix}$. A hypothetical observation is $y=0$, with $Y\mid b\sim N(b_1+b_2,1)$. Thus $L(b)=\frac{1}{\sqrt{2\pi}}\exp\{-(b_1+b_2)^2/2\}$. Use two nodes per standard-normal coordinate.

A Cholesky factor is $C=\begin{bmatrix}1&0\\0.5&\sqrt{0.75}\end{bmatrix}$, so $b_1=z_1$ and $b_2=0.5z_1+\sqrt{0.75}z_2$. The two-point nodes are $z_1,z_2\in\{-1,1\}$, with weight $1/4$ for each pair. Consequently, $b_1+b_2=1.5z_1+\sqrt{0.75}z_2$. Across the four pairs, this sum takes values $\pm2.366025$ and $\pm0.6339746$.

The quadrature estimate is therefore

L^{\mathrm{marg}}\approx \frac{1}{4\sqrt{2\pi}} \left(2e^{-(2.366025)^2/2}+2e^{-(0.6339746)^2/2}\right) \approx 0.17530.

As a check, $b_1+b_2$ is Gaussian with variance $1+1+2(0.5)=3$. Hence $Y$ has marginal distribution $N(0,4)$, and its exact density at zero is $1/\sqrt{2\pi(4)}\approx0.19947$. The two-node product rule gives a low estimate here. This calculation also shows why correlated effects must be transformed with a covariance factor: using independent random effects directly would represent a different model.

Assumptions, diagnostics, and alternatives

The Gaussian transformation and formulas above assume that the random-effects distribution in the integral is Gaussian with the stated mean and covariance. The covariance must be positive definite for the full-dimensional Cholesky transformation shown; a singular covariance requires reducing the integration to its nondegenerate subspace. The conditional likelihood must be integrable under that Gaussian distribution. Gauss–Hermite quadrature is most effective when the integrand is reasonably smooth over the region carrying substantial Gaussian probability. A sharp peak, strong skewness, or substantial likelihood mass far from the Gaussian center can make a fixed, population-centered rule inaccurate.

A basic diagnostic is to recompute the integral using more nodes and compare the marginal likelihood, or its log, at the precision needed for the analysis. Large changes indicate inadequate numerical resolution. Check the node transformation, weight normalization, covariance factor, and whether the calculation averages likelihoods rather than log likelihoods. In log-domain implementations, evaluate the weighted sum stably, especially if likelihood values vary greatly in magnitude. A stable sum prevents avoidable floating-point errors; it cannot repair a quadrature rule that misses important parts of the integrand.

Approach	How it integrates	Useful consideration
Fixed Gauss–Hermite quadrature	Evaluates the integrand at deterministic nodes weighted for a Gaussian distribution.	Repeatable and efficient in low dimensions when the Gaussian-centered nodes resolve the integrand.
Monte Carlo or importance sampling	Uses random draws from a target or proposal distribution and weights their contributions.	Can handle more dimensions flexibly, but estimates have simulation error and require monitoring of weights.
Laplace approximation	Approximates the integrand near a mode using a local Gaussian shape.	Can be economical when the local quadratic approximation is adequate; accuracy depends on the integrand’s shape near the mode.

These methods address the same integration problem differently. Fixed Gauss–Hermite quadrature uses deterministic nodes centered according to the Gaussian random-effects distribution; it does not automatically move those nodes toward an individual subject’s likelihood peak. The introductory material on marginal likelihood integrals establishes why conditional likelihoods must be integrated over random effects. Here, the key numerical choice is whether the fixed Gaussian-centered nodes adequately represent each such integral. The most consequential limitations are the dependence on node placement and the difficulty of resolving complicated integrands as the random-effects dimension grows.

Common pitfalls

Mixing Hermite conventions. Physicists’ nodes for the weight $e^{-x^2}$ are not themselves standard-normal nodes. Use $z_j=\sqrt{2}x_j$ and normalized weights $w_j/\sqrt{\pi}$.
Dropping the Gaussian normalization. The weights for an expectation must sum to one. An unnormalized weighted integral is not the marginal likelihood under a probability distribution.
Using independent coordinates for correlated effects. Apply a factor $C$ satisfying $CC^{\mathsf T}=\Omega$; otherwise the nodes represent the wrong covariance structure.
Assuming more nodes prove accuracy. A rule can still miss a poorly placed or sharply concentrated contribution. Compare orders and check whether estimates stabilize.
Averaging log likelihoods. Marginalize likelihood values first, then take the logarithm of the weighted sum.

Key takeaway: Transform Gaussian random effects to standard-normal coordinates, evaluate the conditional likelihood at the corresponding Gauss–Hermite nodes, and sum using normalized weights. The result is a numerical approximation to the marginal likelihood; its reliability depends on the integrand, the node order, and the correct covariance transformation.

Check Your Understanding

Try each question, then compare with the brief solutions below.

For $b\sim N(1,4)$, what random-effect values correspond to standard-normal nodes $-1$ and $1$?
Why do two-point Gauss–Hermite quadrature weights for an expectation sum to one?
For correlated random effects with covariance $\Omega$, what condition must the transformation matrix $C$ satisfy?
In a marginal likelihood calculation, should one average conditional likelihoods or conditional log likelihoods before taking a logarithm?
What does a substantial change in an integral after increasing the node order suggest?

Solutions. (1) The standard deviation is $2$, so $b=1+2z$ gives values $-1$ and $3$. (2) The normalized weights are $w_j/\sqrt{\pi}$; they sum to one because the rule integrates the constant function under the standard-normal density. (3) It must satisfy $CC^{\mathsf T}=\Omega$. (4) Average the conditional likelihood values with the quadrature weights, then take the logarithm of that weighted sum. (5) It suggests that the integral may not yet be numerically stable with respect to quadrature order; investigate the transformation, integrand, and node resolution before relying on it.