The Importance Sampling Identity and Unbiased Estimators

Why Reweight Samples?

An NLME marginal likelihood averages the conditional data likelihood over the random-effects distribution. When the conditional likelihood is concentrated in a small region of random-effects space, drawing random effects directly from their population distribution may rarely sample that region. Importance sampling rewrites the same integral as an expectation under a chosen proposal density, with a weight that corrects for the difference between the proposal and the target density. The identity is exact; the Monte Carlo estimate of the identity is random.

The distinction between those statements matters. Importance sampling does not, by itself, approximate the model’s marginal likelihood with a different target. It provides a way to estimate the original integral. Whether that estimate is unbiased depends on the sampling and support conditions, and its variability depends on how well the proposal covers the parts of the integral that contribute most.

Marginal Likelihood and the Importance Sampling Identity

For subject \(i\), let \(y_i\) denote the observed data, \(b_i\) the random effects, and \(\theta\) the model parameters. Write \(f(y_i\mid b_i;\theta)\) for the conditional data density and \(p(b_i\mid\theta)\) for the random-effects density. The subject-specific marginal likelihood is

L_i(\theta)=\int f(y_i\mid b_i;\theta)\,p(b_i\mid\theta)\,db_i.

Using standard-normal coordinates, write \(b_i=T_\theta(z)\), where \(z\) has standard-normal density \(\phi_d(z)\) in \(d\) dimensions. The transformation \(T_\theta\) maps those coordinates to the random-effects scale; for a Gaussian random-effects distribution, it can be an affine transformation. Define

h_i(z;\theta)=f(y_i\mid T_\theta(z);\theta).

Then the same marginal likelihood is

L_i(\theta)=\int h_i(z;\theta)\,\phi_d(z)\,dz.

Let \(q(z)\) be a proposal density. It must be a normalized density on the same space as \(z\), and it must be positive wherever \(h_i(z;\theta)\phi_d(z)\) is positive, apart from sets of measure zero. Under this support condition, multiply and divide the integrand by \(q(z)\):

\begin{aligned} L_i(\theta) &=\int h_i(z;\theta)\,\phi_d(z)\, \frac{q(z)}{q(z)}\,dz\\ &=\int \frac{h_i(z;\theta)\phi_d(z)}{q(z)}\,q(z)\,dz\\ &=E_q\!\left[ \frac{h_i(Z;\theta)\phi_d(Z)}{q(Z)} \right]. \end{aligned}

Here \(E_q\) means expectation when \(Z\) has density \(q\). The factor \(\phi_d(Z)/q(Z)\) corrects for drawing from \(q\) rather than from the standard-normal density. Equivalently, the weighted integrand is the conditional likelihood multiplied by the ratio of the target random-effects density to the proposal density.

Definition: The importance-sampling identity expresses the marginal likelihood as \(L_i(\theta)=E_q[W_i(Z;\theta)]\), where \(W_i(z;\theta)=h_i(z;\theta)\phi_d(z)/q(z)\). The identity holds when the proposal covers the nonzero target integrand and the integral is finite.

Unbiased Monte Carlo Estimation

Suppose \(Z_1,\ldots,Z_M\) are independent draws from \(q\). Define one importance weight per draw and average them:

\widehat L_i(\theta) =\frac{1}{M}\sum_{m=1}^{M} \frac{h_i(Z_m;\theta)\phi_d(Z_m)}{q(Z_m)}.

If the marginal likelihood is finite, the support condition holds, and the weights have a finite expectation, then this estimator is unbiased for the marginal likelihood at fixed \(\theta\). The proof uses linearity of expectation and the identity just derived:

\begin{aligned} E_q[\widehat L_i(\theta)] &=\frac{1}{M}\sum_{m=1}^{M} E_q[W_i(Z_m;\theta)]\\ &=\frac{1}{M}\sum_{m=1}^{M}L_i(\theta) =L_i(\theta). \end{aligned}

The independence of the draws is not needed for this expectation calculation, but it is the usual setup for the standard Monte Carlo variance formula. If the draws are independent and the weight has finite second moment, then

\operatorname{Var}_q\!\left(\widehat L_i(\theta)\right) =\frac{1}{M}\left\{ \int \frac{h_i(z;\theta)^2\phi_d(z)^2}{q(z)}\,dz -L_i(\theta)^2 \right\}.

Thus finite expectation is enough for unbiasedness, but a finite second moment is needed for finite variance. A proposal can satisfy the unbiasedness conditions and still produce a very unstable estimate if it assigns too little density to regions that contribute substantially to the integral. Increasing \(M\) reduces the variance at rate \(1/M\) when that variance is finite; it does not repair a proposal that misses target support.

Write the target integral.
Express the marginal likelihood in the coordinates being integrated over, including the random-effects density or its standard-normal equivalent.

Check proposal support.
Confirm that \(q(z)>0\) wherever the target integrand is nonzero, except possibly on sets of measure zero.

Form the corrected weight.
Multiply the conditional likelihood by the target density divided by the proposal density.

Average proposal draws.
The sample average estimates the marginal likelihood unbiasedly when the draws follow \(q\) and the required expectation is finite.

Worked Example: Sampling Directly From the Random-Effects Distribution

Consider a one-dimensional hypothetical model with \(Z\sim N(0,1)\) and \(Y\mid Z=z\sim N(z,1)\). For an observation \(y=0\), the conditional likelihood is \(h(z)=\exp(-z^2/2)/\sqrt{2\pi}\). If the proposal is the random-effects density itself, \(q(z)=\phi(z)\), the density ratio cancels and each importance weight is simply \(h(z)\).

The exact marginal likelihood is the density of \(Y\sim N(0,2)\) evaluated at zero:

L(0)=\frac{1}{\sqrt{4\pi}}=0.2821.

Suppose three hypothetical proposal draws, rounded for display, are \(z_1=0\), \(z_2=1\), and \(z_3=-1\). Their weights are

w_1=\frac{1}{\sqrt{2\pi}}=0.3989,\qquad w_2=w_3=\frac{\exp(-1/2)}{\sqrt{2\pi}}=0.2420.

The estimate is

\widehat L(0)=\frac{0.3989+0.2420+0.2420}{3}=0.2943.

This realized estimate is not equal to \(0.2821\), which is expected for a finite random sample. Unbiasedness describes the average over repeated sets of proposal draws, not a guarantee that any one estimate equals the integral.

Proposal Choice and the Meaning of Unbiasedness

The proposal can change the distribution of the weights without changing their expectation. The following example uses a wider Gaussian proposal to show the density correction explicitly. It also illustrates why unbiasedness alone does not ensure that one run is accurate.

Worked Example: Correcting for a Wider Gaussian Proposal

Use the same hypothetical model and observation \(y=0\), but now take \(q=N(0,2)\). The target random-effects density is \(\phi(z)=N(z;0,1)\), so

\frac{\phi(z)}{q(z)} =\sqrt{2}\exp(-z^2/4).

Multiplying by the conditional likelihood gives the weight function

w(z)=\frac{h(z)\phi(z)}{q(z)} =\frac{1}{\sqrt{\pi}}\exp(-3z^2/4).

For hypothetical rounded draws \(z_1=0\), \(z_2=1\), and \(z_3=-1\), the weights are \(0.5642\), \(0.2665\), and \(0.2665\). Therefore,

\widehat L(0)=\frac{0.5642+0.2665+0.2665}{3}=0.3657.

The exact marginal likelihood remains \(0.2821\): changing \(q\) does not change the target. The estimate \(0.3657\) is a possible finite-sample result, not evidence that the identity changed. The example’s weights vary more than in the idealized proposal proportional to the conditional distribution of \(Z\) given \(Y=0\), for which the weights are constant. In practice, that conditional distribution is generally not known exactly.

In an NLME likelihood with conditionally independent subjects and independent random effects, the full likelihood is the product of subject-specific marginal likelihoods, \(L(\theta)=\prod_i L_i(\theta)\). If each subject’s marginal likelihood is estimated unbiasedly using independent simulation streams across subjects, the product of those estimates is unbiased for the full likelihood: independence makes the expectation of the product equal the product of expectations. If the estimates share random draws, that argument does not apply automatically; their dependence must be accounted for.

Worked Example: Why Taking Logs Changes the Bias Question

Suppose, for illustration, an unbiased estimator of a positive marginal likelihood \(L=0.30\) takes the value \(0.10\) with probability \(1/2\) and \(0.50\) with probability \(1/2\). Its expectation is

E[\widehat L]=\frac{0.10+0.50}{2}=0.30.

But the expectation of its logarithm is

E[\log\widehat L] =\frac{\log(0.10)+\log(0.50)}{2} =-1.498, \qquad \log L=\log(0.30)=-1.204.

The log estimate is lower on average in this example, even though the likelihood estimate is unbiased. In general, for a positive unbiased estimator with an integrable logarithm, Jensen’s inequality gives \(E[\log\widehat L]\leq\log E[\widehat L]=\log L\). Equality holds only when the estimator is constant with probability one. Thus one should not describe the logarithm of an unbiased Monte Carlo likelihood estimate as an unbiased log-likelihood estimate.

Assumptions, Diagnostics, and Alternatives

The identity and estimator rely on a clear set of conditions. The proposal must be a valid density on the integration space; samples must actually be drawn according to that proposal; and its support must cover the nonzero target integrand. The marginal likelihood must be finite for the unbiasedness statement. To use the usual finite-variance calculation, the squared weight must also have finite expectation. These conditions concern a fixed parameter value \(\theta\); if the proposal depends on \(\theta\), the identity still holds at each fixed value when the conditions hold there.

Some warning signs are immediate. A proposal that assigns zero density where the target integrand is positive makes the ratio undefined and omits part of the integral, so the estimator no longer estimates the desired marginal likelihood. A proposal with very small density in an important region can generate occasional very large weights. Then estimates may change substantially across independent runs or as more draws are added. Those patterns call for checking the proposal and its coverage, not simply asserting that the estimator is unbiased.

Direct Monte Carlo under the random-effects distribution is the special case \(q(z)=\phi_d(z)\), in which the density ratio is one. It is simple, but can be inefficient when most of the contribution to the likelihood comes from a region rarely visited under the population distribution. Importance sampling shifts the simulation distribution and corrects it by weighting. Deterministic quadrature instead evaluates the integrand at selected nodes, while Laplace approximations replace the integral by a local approximation. Those are alternatives with different numerical properties: importance sampling preserves the exact integral in expectation under its conditions, but has Monte Carlo variability; quadrature and Laplace approaches do not rely on random proposal draws but can incur approximation or numerical error.

The key limitation is that unbiasedness is a statement about repeated simulation, not about the quality of an individual run. It does not guarantee small variance, good numerical stability, or an unbiased log likelihood. Proposal support, the distribution of weights, and the distinction between estimating \(L_i(\theta)\) and \(\log L_i(\theta)\) must remain explicit when interpreting results.

Common pitfall: “The likelihood estimate is unbiased, so its logarithm is unbiased too.” This is incorrect because the logarithm is nonlinear. For a positive estimator, Jensen’s inequality gives \(E[\log\widehat L]\leq\log L\), with strict inequality when the estimator varies and the expectations exist.

Check Your Understanding

Answer each question, then compare with the solutions below.

Write the importance-sampling identity for \(L_i(\theta)=\int h_i(z;\theta)\phi_d(z)\,dz\) using a proposal \(q(z)\).
For \(M\) independent draws from \(q\), what is the importance-sampling estimator, and what is its expectation?
Why is the condition \(q(z)>0\) wherever \(h_i(z;\theta)\phi_d(z)>0\) important?
If \(E[\widehat L]=L>0\), what inequality relates \(E[\log\widehat L]\) and \(\log L\), when the logarithm is integrable?

Solutions. (1) \(L_i(\theta)=E_q[h_i(Z;\theta)\phi_d(Z)/q(Z)]\). (2) \(\widehat L_i(\theta)=M^{-1}\sum_{m=1}^M h_i(Z_m;\theta)\phi_d(Z_m)/q(Z_m)\); linearity of expectation gives \(E_q[\widehat L_i(\theta)]=L_i(\theta)\). (3) If the proposal is zero in a region with positive target integrand, that region cannot be sampled and its contribution is omitted, so the estimator does not represent the full integral. (4) Jensen’s inequality gives \(E[\log\widehat L]\leq\log E[\widehat L]=\log L\), with equality only if \(\widehat L\) is constant with probability one.