Common Scaling and Jacobian Mistakes in Quadrature

Why Scaling Errors Matter

Quadrature approximates an integral by evaluating its integrand at selected nodes and combining those values with weights. The approximation can be numerically stable and use many nodes, yet still be wrong if the nodes, weights, and change of variables do not describe the same integral. In GLMM likelihood calculations, this kind of error can change a likelihood systematically; increasing the quadrature order will not repair an incorrectly transformed integral.

The central distinction is between a Gaussian-weighted expectation and an ordinary integral. A Gaussian expectation already includes a density that changes under rescaling. An ordinary integral does not: when its integration variable is transformed, the Jacobian must be included. Adaptive quadrature adds another scale, determined by local curvature, which must also agree with the transformed nodes.

Two Integral Forms, Two Different Bookkeeping Rules

Let \(Z\sim N(0,1)\), with density \(\phi(z)=(2\pi)^{-1/2}\exp(-z^2/2)\). If \(B\sim N(0,\sigma^2)\), then \(B=\sigma Z\), where \(\sigma\) is the standard deviation, not the variance. For a function \(f\),

The density transformation and the differential \(db=\sigma\,dz\) cancel in this expectation identity. Thus, if \(x_k,w_k\) are standard-normal quadrature nodes and weights, with weights summing to one, the approximation is \(\sum_k w_k f(\sigma x_k)\). There is no extra factor of \(\sigma\) in that sum.

By contrast, for an integral with respect to ordinary length, \(I=\int q(b)\,db\), setting \(b=a+sz\), for \(s>0\), gives \(db=s\,dz\). Therefore,

Here the factor \(s\) is required. Omitting it, or adding it to an expectation that already uses normalized Gaussian weights, changes the target quantity.

Jacobians for Ordinary and Adaptive Integrals

A useful check for any change of variables is to differentiate the transformation explicitly. For \(b=a+sz\), \(db/dz=s\), so \(db=s\,dz\). For an adaptive transformation, the scale is often chosen from the curvature of the log integrand at its mode. Let \(g(b)\) be twice differentiable, have an interior mode \(\hat b\), and have negative curvature there. Define \(H=-g''(\hat b)>0\) and \(s=H^{-1/2}\). The local change of variables \(b=\hat b+sz\) gives the following exact identity, not an approximation:

The Laplace approximation uses the local quadratic expansion of \(g\) to approximate this integral. Adaptive quadrature instead evaluates the transformed integrand at quadrature nodes; whichever method is used, the transformation’s Jacobian remains \(s\). A curvature-based scale is not interchangeable with its reciprocal.

Transformed Nodes Must Match the Curvature Scale

In adaptive quadrature, a common convention begins with standard-normal nodes \(x_k\) and weights \(w_k\), where \(\sum_k w_k=1\). For the ordinary integral above, or more generally for \(\int \exp\{g(b)\}\,db\), the transformation \(b=\hat b+s z\) yields

The identity follows by multiplying the expectation’s integrand by the standard-normal density: \(\exp(Z^2/2)\phi(Z)=1/\sqrt{2\pi}\). The corresponding quadrature approximation is \(s\sqrt{2\pi}\sum_k w_k\exp\{g(\hat b+sx_k)+x_k^2/2\}\). The same scale \(s\) must appear in the transformed argument and in the Jacobian outside the expectation.

Quadrature conventions also differ. Some Gauss–Hermite tables use nodes \(u_k\) and weights \(W_k\) for the weight function \(\exp(-u^2)\), with \(\int \exp(-u^2)\,du=\sqrt{\pi}\). The standard-normal conversion is \(z=\sqrt{2}\,u\) and normalized weights \(W_k/\sqrt{\pi}\). Treating these \(W_k\) as standard-normal weights fails even for a constant integrand: their sum is \(\sqrt{\pi}\), not 1. Always check the weight convention before combining a table with a formula.

Diagnostics and Common Failure Modes

Scaling mistakes are structural. Refining quadrature order may make the computed value converge more precisely to the wrong integral. The order-refinement and mode checks used in integration diagnostics are therefore most informative after the transformation has been verified algebraically.

• Check a constant Gaussian expectation. With \(f(b)=1\), a normalized Gaussian expectation must equal 1 for every variance. A result such as \(\sigma\), \(\sqrt{\pi}\), or another scale-dependent value signals a density or weight-normalization mismatch.
• Check a Gaussian integral with a known answer. For \(q(b)=\exp\{-(b-a)^2/(2v)\}\), the ordinary integral is \(\sqrt{2\pi v}\). Reproduce it under the same transformation and node convention used in the implementation.
• Check dimensions. An ordinary integral over a variable measured in units of \(b\) has units of \(b\) when its integrand is dimensionless. Its Jacobian supplies those units. A normalized expectation is dimensionless when its integrand is dimensionless.
• Check parameterization invariance. Two correct substitutions of the same integral must agree, provided each uses its own transformed nodes and Jacobian. A difference indicates inconsistent bookkeeping or numerical error.
• Separate numerical error from formula error. If a known Gaussian benchmark remains wrong as the quadrature order rises, investigate the scale, nodes, weights, and Jacobian before increasing the order further.

The main alternative is to use a fixed, nonadaptive rule or a direct Gaussian-expectation formulation. These can simplify bookkeeping when the target distribution is already represented by normalized Gaussian weights. Adaptive rules can concentrate nodes near a conditional mode and account for local curvature, but they introduce additional opportunities to confuse scale conventions. Neither approach is automatically more accurate: correctness of the transformation comes first, and numerical adequacy must then be assessed for the chosen rule.