Monte Carlo Error in Importance-Sampled Likelihoods

Why Monte Carlo Error Needs Its Own Check

An importance-sampled likelihood is a numerical estimate of an integral, not an exact evaluation of it. Even with the model and proposal fixed, different random draws can produce different estimates. In an NLME analysis, that simulation variability can affect an individual’s estimated marginal likelihood and, when likelihood contributions are combined, the objective used for estimation or comparison.

The key distinction is between the likelihood estimate and its logarithm. Under appropriate conditions, averaging importance weights gives an unbiased estimate of the likelihood. Taking the logarithm changes the situation: because the logarithm is concave, the expected log of a random likelihood estimate is generally below the log of the true likelihood. Monte Carlo variability and this logarithmic bias are related, but they are not the same quantity.

Definition: For individual $i$, let $L_i(\theta)$ be the marginal likelihood contribution at parameter value $\theta$: $$ L_i(\theta)=\int p(y_i\mid b_i;\theta)\,p(b_i;\theta)\,db_i. $$ If $b_i^{(m)}$ are independent draws from a proposal density $q_i(b_i)$, define the importance weight $$ W_i^{(m)}= \frac{p(y_i\mid b_i^{(m)};\theta)\,p(b_i^{(m)};\theta)} {q_i(b_i^{(m)})}. $$ The $M$-draw importance-sampling estimate is $\widehat L_i=M^{-1}\sum_{m=1}^M W_i^{(m)}$. Its Monte Carlo error is its random deviation from $L_i(\theta)$ caused by the finite simulation.

Conditions and Weight-Based Monte Carlo Error

The importance-sampling identity requires that the proposal cover the target integrand: wherever $p(y_i\mid b_i;\theta)p(b_i;\theta)$ is positive, $q_i(b_i)$ must also be positive. The draws must follow the stated proposal distribution. Under these conditions, each weight has expectation $L_i(\theta)$, so the sample mean is unbiased. To estimate its variance using the usual sample variance, we additionally require independent draws and a finite weight variance.

Write $\bar W_i=\widehat L_i$, and let $s_{W_i}$ be the sample standard deviation of the $M$ weights. Under the stated conditions, the standard error of the estimated mean is estimated by $s_{W_i}/\sqrt{M}$. This estimates uncertainty in the likelihood scale. It is not itself an uncertainty measure for the log-likelihood.

\widehat{\operatorname{SE}}(\widehat L_i)=\frac{s_{W_i}}{\sqrt{M}}, \qquad s_{W_i}^2=\frac{1}{M-1}\sum_{m=1}^{M}(W_i^{(m)}-\bar W_i)^2.

The relative Monte Carlo standard error, $\widehat{\operatorname{SE}}(\widehat L_i)/\widehat L_i$, expresses the estimated simulation error relative to the estimated likelihood. It can be useful when likelihood contributions have very different scales. The effective sample size (ESS), introduced in the earlier material on importance weights, is another weight-based diagnostic: it summarizes how concentrated the sampled weights are. Neither relative error nor ESS should be treated as an exact guarantee that the integral is accurate. In particular, a finite sample may not reveal rare, extremely large weights.

Worked Example: Standard Error From Five Importance Weights

Suppose five hypothetical draws for one individual produce weights $0.60,\ 0.80,\ 1.00,\ 1.20,\ 1.40$, in likelihood-contribution units. Assume these are independent draws from a proposal that covers the target integrand, and that the weights have finite variance.

The estimated likelihood is their mean: $$ \widehat L_i=\frac{0.60+0.80+1.00+1.20+1.40}{5} =\frac{5.00}{5}=1.00. $$

The squared deviations from the mean sum to $0.16+0.04+0+0.04+0.16=0.40$. Thus the sample variance and sample standard deviation are $$ s_{W_i}^2=\frac{0.40}{5-1}=0.1000, \qquad s_{W_i}=\sqrt{0.1000}=0.3162. $$

The estimated Monte Carlo standard error is $0.3162/\sqrt{5}=0.1414$ likelihood-contribution units. Relative to the estimate $1.00$, this is $14.14\%$. The estimate is unbiased under the assumptions, but five weights give a fairly imprecise result; the calculation does not establish that the true contribution is close to $1.00$.

Independent Replicates and Logarithmic Bias

A direct way to assess simulation variability is to repeat the calculation independently. For replicate $r$, use a fresh set of $M$ proposal draws to obtain $\widehat L_i^{(r)}$. If there are $R$ independent replicates, their sample standard deviation estimates the Monte Carlo standard deviation of an $M$-draw likelihood estimate. The standard error of the average across the $R$ replicates is that sample standard deviation divided by $\sqrt R$. These are distinct summaries: the first describes variability of one $M$-draw result, while the second describes the precision of the average of $R$ such results.

The replicate likelihood estimates can also be logged. The mean of those log estimates describes the average log result across replicates; it is not an unbiased estimate of $\log L_i(\theta)$. The following second-order argument explains why.

Justification: Let $\widehat L_i=L_i+\varepsilon$, where $L_i>0$, $E(\varepsilon)=0$, and $\widehat L_i>0$. When $\varepsilon/L_i$ is small enough for a Taylor approximation around $L_i$, $$ \log(L_i+\varepsilon) \approx \log L_i+\frac{\varepsilon}{L_i}-\frac{\varepsilon^2}{2L_i^2}. $$ Taking expectations and using $E(\varepsilon)=0$ gives $$ E\{\log(\widehat L_i)\} \approx \log L_i-\frac{\operatorname{Var}(\widehat L_i)}{2L_i^2}. $$ The second-order term is nonpositive, so the log estimate has an approximate downward bias. This is an approximation, not an exact bias formula; it is most useful when the likelihood estimate is reasonably concentrated and its variance is finite.

Replacing the unknown likelihood and variance in this expression with estimates gives the plug-in approximation $-\widehat{\operatorname{Var}}(\widehat L_i)/(2\widehat L_i^2)$. The log-scale standard error can also be approximated by the delta method as $\widehat{\operatorname{SE}}(\widehat L_i)/\widehat L_i$. Both are local approximations. When estimates are highly variable, replicate logs are especially useful for showing the actual spread and skewness, but a small number of replicates cannot precisely characterize the tails.

Worked Example: Replicate Logs and the Plug-In Bias Approximation

Consider five hypothetical, independent $M$-draw likelihood estimates for the same individual: $0.80,\ 0.92,\ 1.01,\ 1.08,\ 1.190476$. Assume each replicate uses independent proposal draws and satisfies the support and finite-variance conditions. The slightly more precise final value is retained so that the log calculation can be checked directly.

The replicate mean is $$ \bar L=\frac{0.80+0.92+1.01+1.08+1.190476}{5} =\frac{5.000476}{5}=1.0000952. $$ The logs sum to approximately $-0.0452606$, so their mean is $$ \frac{\log(0.80)+\log(0.92)+\log(1.01)+\log(1.08)+\log(1.190476)}{5} \approx \frac{-0.0452606}{5}=-0.0090521. $$

The sample variance of the five likelihood estimates is approximately $0.022295$. Using the replicate mean as a plug-in estimate of $L_i$, the second-order bias approximation is $$ -\frac{0.022295}{2(1.0000952)^2}\approx -0.01115. $$ The observed mean log estimate is below $\log(\bar L)$, which is approximately $0.0000952$, consistent with the direction predicted by concavity. Its observed difference, about $-0.009147$, is not expected to equal the plug-in approximation exactly: there are only five replicates, and the approximation uses estimated rather than known quantities.

The replicate standard deviation is approximately $\sqrt{0.022295}=0.1493$ likelihood units, describing the variability of an $M$-draw estimate. The standard error of the five-replicate average is $0.1493/\sqrt{5}\approx0.0668$. These values answer different questions and should not be confused with the mean of the replicate logs.

Assessing Error in a Combined Log-Likelihood

For conditionally independent individuals, the marginal likelihood is the product of the individual contributions, and the log-likelihood is the sum of their logarithms. If each individual’s likelihood is estimated by Monte Carlo, the uncertainty in that combined log-likelihood depends on the simulation variability of each contribution and on any dependence between their simulation errors.

With independent simulation draws across individuals, a first-order delta-method approximation gives the variance of each log contribution as $\operatorname{Var}(\widehat L_i)/L_i^2$. Independent contributions then have variances that add. If simulation errors are correlated, covariance terms must also be included; simply summing separate variances would then be incomplete. Replicates that recompute the whole combined log-likelihood can capture its total variation directly, provided each replicate is generated independently using the intended simulation scheme.

Worked Example: Combining Two Individuals’ Log-Scale Error

Suppose a hypothetical analysis uses $R=3$ independent replicates for each of two individuals. Their replicate likelihood estimates are $0.90,\ 1.00,\ 1.10$ for individual 1 and $1.80,\ 2.00,\ 2.20$ for individual 2. Assume the simulation sets for the two individuals are independent and that the likelihood estimates are sufficiently concentrated for a first-order delta approximation.

For individual 1, the mean is $1.00$, and the replicate sample variance is $(0.10^2+0^2+0.10^2)/(3-1)=0.0100$. The approximate variance of its log estimate is therefore $0.0100/(1.00)^2=0.0100$.

For individual 2, the mean is $2.00$, and the replicate sample variance is $(0.20^2+0^2+0.20^2)/(3-1)=0.0400$. Its approximate log-estimate variance is $0.0400/(2.00)^2=0.0100$.

Under independence, the approximate variance of the sum of the two log contributions is $0.0100+0.0100=0.0200$, giving an approximate Monte Carlo standard error of $\sqrt{0.0200}=0.1414$ on the log-likelihood scale. The corresponding second-order downward bias approximation for the sum is $-0.0200/2=-0.0100$. With only three replicates, these are illustrative, imprecise estimates; more independent replicates would be needed for a stable assessment.

Diagnostics, Alternatives, and Common Pitfalls

A useful diagnostic combines more than one view of the simulation. Examine the weight distribution, estimated standard error, relative standard error, and ESS; then check whether independent replicates tell a consistent story. Large changes in the likelihood or log-likelihood across replicates, a few weights dominating the estimate, or a large estimated error relative to the likelihood scale all suggest that the chosen simulation budget or proposal is not yielding a stable estimate. Increasing $M$ can reduce Monte Carlo error when the assumptions hold, but it does not correct a proposal that fails to cover important regions of the integrand.

Independent replicates are an empirical check, not a replacement for weight diagnostics. Replicate variation can be deceptively small if all runs miss a rare region that matters to the integral. Conversely, a very small ESS or a large weight-based standard error flags instability even if a few replicate estimates happen to be similar. When the likelihood estimate is positive but highly skewed, the log estimates may also be asymmetric, making a normal approximation to their uncertainty unreliable.

The main practical alternatives are to spend more draws with the same proposal, improve the proposal, or use a deterministic integration method where feasible. More draws target finite-simulation variability but may be inefficient if weights remain highly variable. A better proposal can reduce weight dispersion, while deterministic integration avoids random-draw error but has its own numerical limitations. The method choice does not remove the need to assess accuracy. In particular, an approximate likelihood from one method should not be treated as exact merely because its computation is repeatable.

Pitfall: treating an unbiased likelihood estimate as an unbiased log-likelihood. Unbiasedness means $E(\widehat L_i)=L_i$; it does not imply $E\{\log(\widehat L_i)\}=\log L_i$. Concavity generally makes the expected log smaller.
Pitfall: reporting $s_W/\sqrt M$ as a log-likelihood standard error. This quantity is on the likelihood scale. A first-order log-scale approximation divides it by $\widehat L_i$.
Pitfall: treating ESS as a guarantee. ESS summarizes the observed weight concentration, but cannot certify that an unobserved, influential region was sampled.
Pitfall: confusing replicate spread with the precision of the replicate average. The replicate standard deviation estimates the variability of one run; dividing it by $\sqrt R$ estimates the standard error of the average across $R$ runs.

Key takeaway: Use weight-based standard errors and independent replicates to assess simulation variability, and keep the likelihood and log-likelihood scales distinct. The importance-sampling likelihood estimate can be unbiased while its logarithm has a downward bias whose size depends, approximately, on relative variance.

Check Your Understanding

Answer each question, then compare with the solutions.

What two conditions make the sample mean of importance weights an unbiased estimate of a marginal likelihood contribution?
If the estimated likelihood standard error is $0.12$ and the estimated likelihood is $2.00$, what is the first-order approximate log-scale standard error?
Why does an unbiased likelihood estimate not imply an unbiased log-likelihood estimate?
Across four independent replicates, the standard deviation of the likelihood estimates is $0.08$. What does $0.08$ estimate, and what is the estimated standard error of their average?
Why do similar independent replicate estimates not, by themselves, prove that the importance-sampling calculation is accurate?

Solutions. (1) The proposal must have positive density wherever the target integrand is positive, and the draws must follow that proposal; then the importance weights have expectation equal to the integral. (2) The delta-method approximation is $0.12/2.00=0.060$ on the log scale. (3) The logarithm is concave, so the second-order expansion includes the nonpositive term $-\operatorname{Var}(\widehat L)/(2L^2)$; unbiasedness on the likelihood scale does not carry through the nonlinear transformation. (4) The value $0.08$ estimates the standard deviation of one replicate’s likelihood estimate. The standard error of the four-replicate average is $0.08/\sqrt{4}=0.040$. (5) All replicates can miss the same important region of the integrand; weight diagnostics and the proposal’s coverage assumptions must also be considered.