Why Importance Sampling Is Used
In a nonlinear mixed-effects model, a subject's likelihood contribution averages the conditional data likelihood over the distribution of that subject's random effects. The resulting integral is often not available in closed form. Importance sampling estimates it by drawing random effects from a chosen proposal distribution and reweighting each draw so that the average still targets the original integral.
The proposal distribution matters. If it rarely draws random effects from regions that contribute substantially to the likelihood, a few unusually large weights may dominate the estimate. Increasing the number of draws can reduce simulation error, but it does not make a poorly chosen proposal efficient automatically.
The proposal must be a normalized density, and it must be positive wherever the target integrand \(f(y_i\mid\eta_i;\theta)g(\eta_i;\theta)\) is positive. This support condition is essential: if a region contributes to the integral but the proposal assigns it zero probability, no importance weight can recover its contribution. A proposal may depend on the observed data and on the current parameter values, provided the corresponding density is used in the denominator of every weight.
Why the Weighted Average Targets the Integral
The estimator follows from changing the density used to express the integral. Write the target integrand as \(h_i(\eta_i;\theta)=f(y_i\mid\eta_i;\theta)g(\eta_i;\theta)\). Wherever \(q_i\) is positive, multiply and divide by \(q_i\), then take an expectation under the proposal:
Thus each weight has expectation \(L_i(\theta)\), and the average of independent weights is an unbiased estimate of the likelihood contribution. If the weights have finite variance, the variance of that average is
This result explains the role of both proposal and sample count. A proposal that produces less variable weights reduces the variance for a given \(M\). Under a fixed proposal with finite weight variance, multiplying \(M\) by four divides the estimator's variance by four and its standard error by two. The result does not promise good precision when the weight variance is extremely large or infinite.
For independent subject-level draws, the estimated total log likelihood is usually formed as \(\sum_i \log \widehat L_i\). Although \(\widehat L_i\) is unbiased for \(L_i\), its logarithm is generally not unbiased: because the logarithm is concave, Jensen's inequality gives \(E[\log \widehat L_i]\leq \log L_i\). This matters when simulation noise is substantial; a stable-looking likelihood estimate should not be treated as exact merely because it is expressed on the log scale.
Worked Example: Estimating One Subject's Likelihood Contribution
Suppose a hypothetical one-dimensional random-effects likelihood integral is evaluated using four independent draws from a proposal. At each draw, the target integrand \(h_i=f(y_i\mid\eta_i)g(\eta_i)\) and proposal density \(q_i\) have the following values. The values are illustrative density evaluations, not observations from a real study.
| Draw | \(h_i(\eta_i^{(m)})\) | \(q_i(\eta_i^{(m)})\) | Weight \(h_i/q_i\) |
|---|---|---|---|
| 1 | 0.008 | 0.20 | 0.040 |
| 2 | 0.018 | 0.30 | 0.060 |
| 3 | 0.006 | 0.15 | 0.040 |
| 4 | 0.005 | 0.25 | 0.020 |
The estimate is the arithmetic mean of the four weights:
The estimated likelihood contribution is 0.0400 in the units of the conditional data density. The proposal density has been accounted for in each weight; averaging the target integrand values alone would not implement importance sampling. This estimate is based on only four draws, so it should not be interpreted as highly precise without further assessment.
Proposal Choice and Weight Stability
An effective proposal places draws in regions that matter to the integral and has tails broad enough to cover the target integrand. In a random-effects likelihood, the relevant region is often where the product of the conditional data likelihood and the random-effects density is appreciable. Sampling only from the random-effects distribution can be inefficient when the subject's data make the important region much narrower or shift it away from the population center.
A useful theoretical benchmark is a proposal proportional to the target integrand: \(q_i(\eta_i)\propto h_i(\eta_i)\). If that normalized density were available, the weights would be constant and the estimate would have zero sampling variance. In practice, its normalizing constant is the unknown integral itself, so this ideal is usually unavailable. It nevertheless explains why proposals that approximate the likelihood-relevant distribution can improve efficiency.
Worked Example: Comparing a Prior Proposal With a Data-Focused Proposal
Consider a hypothetical random effect with two possible values, common and rare. The random-effects probabilities are 0.99 and 0.01, and the conditional likelihoods for one subject's data are 0.01 and 0.80, respectively. The exact marginal contribution is
If draws come from the random-effects distribution itself, the proposal probabilities are 0.99 and 0.01. The importance weights equal the conditional likelihoods, 0.01 for the common value and 0.80 for the rare value. With 100 draws, the probability of drawing no rare value is
When no rare value is drawn, the estimate is 0.0100, well below the exact contribution 0.0179. The missing contribution is not corrected by the weights because no draw has reached the region responsible for it.
Now use a data-focused proposal that assigns probability 0.50 to each value. The weights are \((0.99)(0.01)/0.50=0.0198\) for the common value and \((0.01)(0.80)/0.50=0.0160\) for the rare value. Suppose ten independent draws happen to include six common and four rare values. Then
This estimate is close to 0.0179, and its weights are similar in size. The example illustrates the principle rather than a universal recipe: a proposal should allocate draws to important regions, but the target-to-proposal ratio must still be applied. A proposal that samples from the target integrand without correcting the weights would generally estimate a different quantity.
Assessing Monte Carlo Precision
For a set of \(M\) independent weights, let \(s_w^2\) be their sample variance. An estimated Monte Carlo standard error for the likelihood estimate is \(s_w/\sqrt{M}\). Dividing this by the estimated likelihood gives an estimated relative standard error. These quantities describe simulation variability under the proposal and model used; they do not measure model misspecification or uncertainty from the data.
The effective sample size (ESS) is a weight-concentration diagnostic:
For equal positive weights, ESS equals \(M\). When a few weights dominate, ESS is smaller. It is a diagnostic summary, not a guarantee that the estimate is accurate: draws may have missed an important region entirely, and a small ESS signals that the estimate may be unstable rather than providing a corrected estimate.
Worked Example: How Sample Count Affects Estimated Precision
Suppose a hypothetical pilot calculation produces four independent weights: 0.8, 1.0, 1.2, and 1.0. Their mean is
The deviations from the mean are \(-0.2, 0, 0.2, 0\). The sample variance and estimated standard error are therefore
The estimated relative standard error is \(0.08165/1.0=8.17\%\). The pilot weight coefficient of variation is \(\sqrt{0.02667}/1.0=0.1633\). If this variability remained representative as the number of independent draws increased, the approximate relative standard error would be \(0.1633/\sqrt{M}\): about 1.63% at \(M=100\) and 0.817% at \(M=400\). These are scaling estimates, not guaranteed errors for a future run; the actual weights can vary, especially when rare, large weights occur.
Numerical Practice and Failure Modes
When weights are very unequal, calculations on the log scale help avoid numerical underflow. If \(\ell_m=\log w_i^{(m)}\), the log estimate can be computed as
Subtracting the largest log weight keeps the exponentials at or below one; the adjustment by \(a\) restores the original sum. This stabilizes arithmetic but does not fix an inadequate proposal or reduce the sampling variability of the weights.
Warning signs include a very large maximum weight relative to the others, low ESS, materially different estimates across independent runs, and log-likelihood contributions that change substantially when the draw count is increased. Repeating an estimate with a different random seed can help reveal instability, but agreement across a few runs does not prove that all important regions have been sampled.
Importance sampling is one approach to likelihood evaluation, not the only one. Laplace approximation replaces the integral by a local approximation around a mode, while numerical quadrature evaluates the integrand at structured points. Importance sampling is flexible in dimension and can represent a non-Gaussian target through weighted draws, but its precision depends strongly on the proposal and the tails of the weights. These approaches have different approximation or numerical errors; none should be assumed accurate solely because an optimizer returns a result.
A common mistake is to argue that a large draw count guarantees a reliable estimate. A large \(M\) reduces the variance at the rate \(1/M\) only when the weight variance is finite and the draws genuinely cover the relevant support. If the proposal almost never visits an important region, the finite-run estimate can remain misleading. The correction is to check support and weight behavior, examine estimated Monte Carlo error and ESS, and determine whether the likelihood estimate is stable to additional draws or a more suitable proposal.
Check Your Understanding
Answer each question, then compare your reasoning with the solutions below.
- State the support condition required of a proposal density for importance sampling to estimate a random-effects likelihood integral.
- For target integrand values 0.006 and 0.012 and proposal density values 0.20 and 0.30 at two draws, calculate the two weights and the estimated integral.
- If independent weights have finite variance, what happens to the standard error when the draw count is multiplied by nine?
- Why can a low ESS be concerning, and why does it not prove that the estimate is biased?
Solutions. (1) The proposal must be positive everywhere the target integrand is positive. (2) The weights are \(0.006/0.20=0.030\) and \(0.012/0.30=0.040\); their mean is \((0.030+0.040)/2=0.035\). (3) Since standard error scales as \(1/\sqrt{M}\) under the stated conditions, multiplying \(M\) by nine divides the standard error by three. (4) Low ESS indicates that a small number of weights dominate, so the estimate may have high simulation variability. ESS is a diagnostic, not a test of bias; a low value alone does not establish that the estimate is biased.