Stable Log-Domain Computation of Marginal Likelihoods

Why Compute Marginal Likelihoods on the Log Scale?

In a nonlinear mixed-effects model, a subject’s marginal likelihood integrates the conditional data likelihood over the distribution of that subject’s random effects. Importance sampling replaces this integral with an average of weighted contributions. Those contributions can be extremely small, extremely large, or spread across many orders of magnitude. Computing them first on the original scale can therefore produce underflow, overflow, or loss of precision, even when the final log-likelihood is a reasonable finite number.

The remedy is to keep contributions as log weights for as long as possible. The central calculation is the log-sum-exp identity: it evaluates the logarithm of a sum of exponentials by subtracting the largest exponent before exponentiating. This tutorial develops the calculation for marginal likelihood estimates, shows how to normalize log weights, and explains how to identify numerical failures.

Definition: For subject $i$, let $y_i$ denote the observed data, $b_i$ the random effects, and $\theta$ the population and residual-model parameters. The marginal likelihood contribution is $$ L_i(\theta)=\int p(y_i\mid b_i,\theta)\,p(b_i\mid\theta)\,db_i. $$ If random effects are sampled from a proposal density $q_i(b_i)$, define the importance weight at draw $b_i^{(m)}$ by $$ w_i^{(m)}=\frac{p(y_i\mid b_i^{(m)},\theta)\,p(b_i^{(m)}\mid\theta)} {q_i(b_i^{(m)})}, \qquad \ell_i^{(m)}=\log w_i^{(m)}. $$

The numerator and proposal density must be evaluated with respect to the same coordinates and measure. If a transformation of random effects is used, its density transformation must be included so that the ratio remains correct. The proposal must assign positive density almost everywhere with respect to the measure defined by the integrand $p(y_i\mid b_i,\theta)p(b_i\mid\theta)$ (equivalently, the integrand measure must be absolutely continuous with respect to the proposal measure). Otherwise, a set of positive integrand measure is omitted and the importance-sampling identity does not hold.

For independent draws from $q_i$, the usual Monte Carlo estimate is $\widehat L_i=M^{-1}\sum_{m=1}^M w_i^{(m)}$. The weights are nonnegative, and under the support condition and a finite integral, this estimate is unbiased for $L_i$. Its logarithm is the quantity used in a log-likelihood objective, but the logarithm of an unbiased estimate is not generally an unbiased estimate of $\log L_i$.

The Log-Sum-Exp Calculation

Suppose the computed log weights for one subject are $\ell^{(1)},\ldots,\ell^{(M)}$. Directly evaluating each $\exp(\ell^{(m)})$ may overflow or underflow. Define $a=\max_m \ell^{(m)}$. Factoring $\exp(a)$ out of the sum gives the stable identity below.

\log\left(\sum_{m=1}^M \exp(\ell^{(m)})\right) = a+\log\left(\sum_{m=1}^M \exp(\ell^{(m)}-a)\right).

Each shifted value $\ell^{(m)}-a$ is nonpositive, and at least one equals zero. Thus no exponential exceeds one, and the shifted sum is at least one. The log of the Monte Carlo estimate is

\log\widehat L_i = a+\log\left(\sum_{m=1}^M \exp(\ell_i^{(m)}-a)\right)-\log M.

The subtraction of $\log M$ matters: the sum estimates $M L_i$, while the Monte Carlo estimate is an average. Omitting that term inflates the estimate by a factor of $M$. The identity is exact in real arithmetic; finite-precision arithmetic still has limits, but shifting prevents the most common overflow failure and greatly reduces underflow risk.

Worked Example: Log Weights Near Zero on the Original Scale

Suppose a hypothetical subject has three importance-sampling log weights $-20$, $-21$, and $-25$. Compute the estimated marginal likelihood and its logarithm.

Here $M=3$ and $a=-20$. The shifted sum is $$ \exp(0)+\exp(-1)+\exp(-5) =1+0.367879+0.006738 =1.374617. $$

Therefore, $$ \log\widehat L_i =-20+\log(1.374617)-\log(3) \approx -20.7804. $$ On the original scale, $$ \widehat L_i =\frac{\exp(-20)+\exp(-21)+\exp(-25)}{3} =\exp(-20)\frac{1+\exp(-1)+\exp(-5)}{3} \approx 9.4443\times10^{-10}. $$

The estimate is positive and small, as expected from the log weights. The rounded original-scale result is approximately $9.44\times10^{-10}$, not zero. The log-scale result preserves a usable value without requiring the small exponentials to be summed directly.

Normalizing Log Weights

Some calculations require weights that sum to one, for example when computing a weighted average over sampled random effects. The normalized weight is $\widetilde w^{(m)}=w^{(m)}/\sum_r w^{(r)}$. In log form, its logarithm is the log weight minus the log-sum-exp of all log weights. Equivalently, subtract the largest log weight, exponentiate the shifted values, and divide by their shifted sum.

\log \widetilde w^{(m)} = \ell^{(m)}-\operatorname{LSE}(\ell^{(1)},\ldots,\ell^{(M)}), \qquad \widetilde w^{(m)} = \frac{\exp(\ell^{(m)}-a)} {\sum_{r=1}^M\exp(\ell^{(r)}-a)}.

Here $\operatorname{LSE}$ denotes the log-sum-exp. The normalized weights are nonnegative and sum to one, up to rounding. Their scale is relative: adding the same constant to every log weight changes neither the normalized weights nor their ordering.

Worked Example: Normalizing Very Large Log Weights

For three hypothetical draws, suppose the log weights are $1000$, $999$, and $998$. Direct exponentiation would overflow in ordinary floating-point arithmetic. Normalize them using $a=1000$.

The shifted exponentials are $1$, $\exp(-1)\approx0.367879$, and $\exp(-2)\approx0.135335$. Their sum is approximately $1.503215$. The normalized weights are therefore $$ \widetilde w^{(1)}=\frac{1}{1.503215}\approx0.6652,\quad \widetilde w^{(2)}=\frac{0.367879}{1.503215}\approx0.2447,\quad \widetilde w^{(3)}=\frac{0.135335}{1.503215}\approx0.0900. $$

They sum to approximately $1.000$, with small discrepancies possible from rounding. Their values depend only on the differences among the log weights, not on the common offset of $1000$. This calculation is valid for normalizing weights; it does not by itself give the marginal likelihood estimate, which also requires retaining the log-sum scale and subtracting $\log M$.

Extreme Values, Underflow, and Failure Checks

Log-domain computation is especially important when all contributions are tiny. If each $\ell^{(m)}$ is around $-1000$, direct exponentiation may return zero in standard floating-point arithmetic. A program may then report a zero likelihood and an undefined log likelihood, even though the correct log estimate is finite. Shifting by the largest log weight avoids this failure.

Worked Example: Contributions That Underflow Directly

Consider hypothetical log weights $-1000$, $-1001$, and $-1010$. Directly exponentiating them in ordinary double precision may produce three zeros. Instead, set $a=-1000$. The shifted sum is $$ 1+\exp(-1)+\exp(-10) \approx1+0.367879+0.0000454 =1.367925. $$

The stable log estimate is $$ \log\widehat L_i =-1000+\log(1.367925)-\log(3) \approx -1000.7853. $$

This is a finite log estimate even though representing $\widehat L_i$ itself may underflow to zero. If downstream calculations need the likelihood, retain its logarithm rather than exponentiating it merely to take a logarithm again.

After shifting, an exponential whose difference from the maximum is extremely negative can still underflow to zero. In a finite-precision calculation its relative contribution to the sum is then negligible at the available precision; the largest term remains represented as one. Do not confuse this harmless loss of a numerically negligible term with the serious failure of every original-scale contribution underflowing before summation.

Several edge cases need explicit handling. A log weight of $-\infty$ represents a zero contribution and can be included; if all log weights are $-\infty$, the estimated likelihood is zero and its log is $-\infty$. A NaN indicates an invalid computation and should not be silently treated as a weight. If a log weight is $+\infty$, subtracting the maximum can produce the indeterminate expression $+\infty-(+\infty)$; investigate the density calculations or use explicit special-case logic rather than applying the ordinary formula blindly.

Derivation and Practical Diagnostics

The derivation requires only factoring the largest exponential from a positive sum. Since $\ell^{(m)}=a+(\ell^{(m)}-a)$, exponentiation gives $\exp(\ell^{(m)})=\exp(a)\exp(\ell^{(m)}-a)$. Summing and taking logarithms yields the identity. For the Monte Carlo average, taking the logarithm of $M^{-1}\sum_m\exp(\ell^{(m)})$ adds the term $-\log M$. This also explains why a common shift $c$ in all log weights adds $c$ to the log estimate but leaves normalized weights unchanged.

Implementation checks should be built into the calculation:

Confirm that the log likelihood, random-effects log density, and proposal log density use the same draw and compatible coordinates. The log weight is their combination $\log p(y_i\mid b_i,\theta)+\log p(b_i\mid\theta)-\log q_i(b_i)$.
Check that the computed log estimate is finite when the data and proposal imply a positive finite estimate. A zero original-scale sum with finite log weights is a warning that direct exponentiation was used too early.
Verify that normalized weights sum to one within a numerical tolerance and that no weight is negative. Inspect the largest normalized weight; a value near one signals that one draw dominates the normalized sample.
Compare estimates across independent runs or larger simulation sizes. Numerical stability prevents arithmetic failure, but it does not guarantee small Monte Carlo error.
For a numerical check on concentration, the effective sample size can be computed stably from log weights as $\exp\{2\operatorname{LSE}(\ell)-\operatorname{LSE}(2\ell)\}$. This does not repair a poor proposal; it flags when few draws carry most of the total weight.

For independent subjects, the model’s joint log likelihood is the sum of subject-specific log likelihood contributions. If each contribution is estimated by simulation, summing their computed log estimates gives the simulated objective. However, because logarithms of random estimates are generally biased, that objective is not automatically an unbiased estimate of the exact joint log likelihood. Keep this distinction separate from the unbiasedness of the likelihood-scale importance estimate.

Alternatives and Limitations

Approach	Advantage	Limitation
Directly exponentiate and sum weights	Simple for values safely within the numerical range	Can underflow to zero, overflow, or lose small contributions
Log-sum-exp with maximum shifting	Stable for very small or large finite log weights; supports stable normalization	Cannot fix incorrect density ratios, inadequate sampling, or invalid values such as NaN
Deterministic numerical integration	A different way to approximate the marginal integral without random importance draws	Its accuracy and computational cost depend on the integral and the chosen integration method

Log-sum-exp is a numerical safeguard, not a statistical correction. It cannot compensate for a proposal that misses important regions, a coding error in the proposal density, or high Monte Carlo variability. Nor should every tiny normalized weight be interpreted as a failure: a broad range of weights may be mathematically correct, although severe concentration can indicate that the simulation is inefficient. The useful practice is to preserve log weights, inspect their spread and normalized concentration, and assess simulation error separately from floating-point stability.

Key takeaway: Keep importance contributions as log weights, subtract their maximum before exponentiating, and retain the log-sum scale when calculating the likelihood estimate. For normalized weights, use the same shifted values but divide by their shifted sum. Stable arithmetic prevents avoidable underflow and overflow; it does not guarantee a good proposal or precise Monte Carlo estimate.

Check Your Understanding

Answer each question, then compare with the solution immediately below it.

For log weights $-3$ and $-5$, write the stable expression for the log Monte Carlo estimate with $M=2$.
Solution: The maximum is $-3$, so the shifted sum is $1+\exp(-2)$. Thus $\log\widehat L=-3+\log(1+\exp(-2))-\log 2$, approximately $-3.566$.
Why is it incorrect to omit $-\log M$ from the log estimate?
Solution: Log-sum-exp without that term gives the logarithm of $\sum_m w^{(m)}$, not the logarithm of the average. Since $\widehat L=M^{-1}\sum_m w^{(m)}$, its logarithm includes $-\log M$.
All computed log weights for a subject equal $-\infty$. What does the calculation imply, and what should be checked?
Solution: Every represented contribution is zero, so the estimated likelihood is zero and its log is $-\infty$. Check whether the integrand is truly zero at every sampled draw or whether log densities underflowed, invalid values were introduced, or the proposal failed to cover relevant regions.
If the same finite constant is added to every log weight, what happens to the normalized weights and the log estimate?
Solution: The normalized weights do not change, because the common factor cancels between numerator and denominator. The log estimate increases by that constant, since the likelihood-scale weights have all been multiplied by its exponential.