Tutorials › Nonlinear Mixed-Effects Models (NLME) › Adaptive Importance Sampling by Weighted Moment Updates

Importance Sampling and Adaptive Quadrature · Tutorial 184 of 797

Adaptive Importance Sampling by Weighted Moment Updates

You will learn to update a Gaussian importance proposal from weighted samples while controlling instability and separating adaptation from reliable likelihood estimation.

Advanced 11 min read

What You'll Learn

  • Define normalized importance weights and use them to estimate target means and covariances.
  • Derive a damped Gaussian proposal update and calculate its weighted moments.
  • Apply safeguards for small effective sample sizes, singular covariance estimates, and missed modes.
  • Distinguish adaptation samples from independent samples used for final likelihood estimation.
  • Diagnose whether adaptation improved overlap without treating a larger effective sample size as proof of accuracy.

Why Adapt an Importance Proposal?

For an individual in a nonlinear mixed-effects model, a marginal likelihood contribution is an integral over that individual's random effects. In coordinates \(z\), write its unnormalized integrand as \(g(z)\), so the quantity of interest is

$$ Z=\int g(z)\,dz. $$

An importance proposal \(q(z)\) lets us estimate this integral by drawing \(z_j\) from \(q\) and averaging the importance weights \(w_j=g(z_j)/q(z_j)\). If the proposal places little probability where \(g\) is substantial, a few draws can dominate the estimate. Earlier material on importance weights and effective sample size, and on heavy-tailed weights and defensive importance proposals, provides tools for recognizing and reducing such problems.

Adaptive importance sampling uses pilot draws and their weights to move a proposal toward regions that matter under the normalized target \(g(z)/Z\). The aim here is moment-based adaptation: estimate the target's mean and covariance by weighted averages, then use those estimates to update a Gaussian proposal. Adaptation changes the sampling distribution, not the target integral.

Definition: Given draws \(z_1,\ldots,z_n\) from a proposal \(q\), define \(w_j=g(z_j)/q(z_j)\) and \(\tilde w_j=w_j/\sum_{k=1}^n w_k\). The weighted estimates of the target mean and covariance are $$ \hat m=\sum_{j=1}^n\tilde w_jz_j,\qquad \hat C=\sum_{j=1}^n\tilde w_j(z_j-\hat m)(z_j-\hat m)^\mathsf{T}. $$

These are self-normalized estimates: the unknown normalizing constant \(Z\) cancels when weights are normalized. The covariance uses denominator \(\sum_j\tilde w_j=1\); it is not an ordinary sample covariance with an \(n-1\) correction. With suitable support and moment conditions, these weighted estimates converge to the mean and covariance of the normalized target as the number of draws grows. At finite sample sizes they can be noisy, especially when only a few draws carry most of the weight.

The Moment Update and Its Justification

Suppose the current proposal at iteration \(t\) is Gaussian, \(q_t(z)=N(m_t,C_t)\). Draw a pilot batch from it, calculate importance weights using \(q_t\), and estimate \(\hat m_t\) and \(\hat C_t\). The simplest update sets the next proposal's mean and covariance to those estimates. In practice, taking the entire step can be unstable, so a damped update is often used:

$$ m_{t+1}=(1-\lambda)m_t+\lambda\hat m_t,\qquad C_{t+1}=(1-\lambda)C_t+\lambda\hat C_t,\qquad 0<\lambda\leq 1. $$

The parameter \(\lambda\) controls how quickly the proposal changes. A smaller value retains more of the current proposal; \(\lambda=1\) makes the direct moment-matching update. Both covariance matrices in the blend should be symmetric positive definite. A convex combination of positive-definite matrices is positive definite, but the estimated \(\hat C_t\) may fail to be positive definite when the weighted points do not span the integration space. That situation calls for a safeguard, not an unexamined update.

Justification: For any integrable function \(h\), importance sampling gives \(E_q[w(Z)h(Z)]=\int h(z)g(z)\,dz\), while \(E_q[w(Z)]=Z\). Their ratio is the target expectation \(E_g[h(Z)/Z]\). Replacing these expectations by sample averages gives the self-normalized estimate \(\sum_j\tilde w_j h(z_j)\). Taking \(h(z)=z\) estimates the target mean; using the estimated mean in the weighted outer products estimates its covariance. Thus the update targets moments of the integrand's normalized distribution, rather than moments of the proposal draws alone.

This justification also identifies important conditions. The proposal must cover the target's support: wherever \(g(z)>0\), the proposal must allow sampling there. The weighted moments must exist, and in practice enough draws must reach important target regions for the estimates to be informative. Moment matching is not a guarantee of an accurate marginal likelihood. A Gaussian with the target's first two moments can still represent a skewed or multimodal target poorly.

Worked Example: A Damped One-Dimensional Update

Suppose a pilot batch produces locations \(z=(-1,0,2)\) with importance weights \(w=(1,2,1)\). The current proposal has mean \(m_t=0\) and variance \(C_t=1\), and the update uses \(\lambda=0.5\). These are hypothetical pilot values.

The total weight is \(1+2+1=4\), so the normalized weights are \((0.25,0.50,0.25)\). The weighted mean is

$$ \hat m_t=0.25(-1)+0.50(0)+0.25(2)=0.25. $$

The weighted variance around that mean is

$$ \hat C_t=0.25(-1-0.25)^2+0.50(0-0.25)^2+0.25(2-0.25)^2 =0.390625+0.03125+0.765625=1.1875. $$

Damping gives \(m_{t+1}=0.5(0)+0.5(0.25)=0.125\) and \(C_{t+1}=0.5(1)+0.5(1.1875)=1.09375\). The updated proposal is \(N(0.125,1.09375)\), with variance in squared \(z\)-units. The effective sample size for this batch is \((\sum w_j)^2/\sum w_j^2=16/6=2.667\) out of three draws. That small batch is adequate for illustrating the arithmetic, but it gives little basis for trusting the update in a real calculation.

Multivariate Updates and Numerical Safeguards

For a vector of random effects, the weighted covariance contains both marginal variances and cross-covariances. Preserving those cross-covariances can matter: an updated Gaussian proposal with an appropriate covariance can follow an elongated target region better than one that only adjusts coordinate-wise variances. The mechanics are the same as in one dimension, but a covariance estimate based on too few influential points is more likely to be singular or unstable.

Worked Example: Updating a Two-Dimensional Proposal

Consider hypothetical pilot points \(z_1=(0,0)\), \(z_2=(2,0)\), and \(z_3=(0,1)\), with weights \(1,2,1\). Let the current proposal have \(m_t=(0,0)\) and \(C_t=I\), the two-dimensional identity matrix, and use \(\lambda=0.5\).

The normalized weights are \((0.25,0.50,0.25)\), giving

$$ \hat m_t=0.25(0,0)+0.50(2,0)+0.25(0,1)=(1,0.25). $$

For the first coordinate, the weighted variance is \(0.25(0-1)^2+0.50(2-1)^2+0.25(0-1)^2=1\). For the second, it is \(0.25(0-0.25)^2+0.50(0-0.25)^2+0.25(1-0.25)^2=0.1875\). The weighted cross-covariance is \(0.25(0-1)(0-0.25)+0.50(2-1)(0-0.25)+0.25(0-1)(1-0.25)=-0.25\). Therefore,

$$ \hat C_t= \begin{pmatrix} 1&-0.25\\ -0.25&0.1875 \end{pmatrix}. $$

The damped update is \(m_{t+1}=(0.5,0.125)\) and

$$ C_{t+1}= \begin{pmatrix} 1&-0.125\\ -0.125&0.59375 \end{pmatrix}. $$

This matrix is positive definite: its first diagonal entry is positive and its determinant is \(1(0.59375)-(-0.125)^2=0.578125>0\). The updated proposal therefore has a valid covariance matrix. The example also has effective sample size \(16/(1+4+1)=2.667\), again a warning that these illustrative moments rest on very limited weighted information.

Useful safeguards address different failure modes; none can recover target regions that the pilot never visits.

  • Damp the update. Choose \(\lambda<1\) when moment estimates vary noticeably across pilot batches. Avoid a large move based on a low effective sample size.
  • Regularize the covariance. If \(\hat C_t\) is singular or nearly singular, add a small positive diagonal term or impose a lower bound on its eigenvalues before using it. The scale of the bound should be meaningful in the coordinates being sampled.
  • Retain adequate spread. Do not allow a proposal covariance to collapse merely because the weighted pilot points are tightly clustered. A covariance floor or a retained fraction of the previous covariance can help keep exploration possible.
  • Use a defensive component when appropriate. The earlier subject of heavy-tailed weights and defensive importance proposals covers mixtures that retain broad support. Moment adaptation should not remove that protection if it is needed.

When Weighted Moments Mislead

A weighted mean and covariance summarize the points that were sampled and received appreciable weights. If the proposal misses an important mode, the pilot cannot reveal that mode. In that case, the estimated moments can be numerically well behaved and still describe only part of the target. This is a fundamental limitation of adaptation based on one pilot cloud, not a defect that can be corrected by algebra.

Worked Example: A Pilot That Misses a Mode

Suppose the normalized target is a hypothetical equal mixture of two narrow components centered at \(z=-3\) and \(z=3\). Its symmetry gives target mean \(0\), and if each component has variance \(0.25\), its variance is \(9.25\): the within-component variance \(0.25\) plus the between-component variance \(9\).

Imagine instead that a pilot proposal samples only near the left component and returns \(z=(-3.5,-3,-2.5)\) with equal weights. Its weighted mean is \((-3.5-3-2.5)/3=-3\), and its weighted variance is \((0.25+0+0.25)/3=1/6\). A direct moment update would center the next Gaussian at \(-3\) with variance \(1/6\), despite the target's mean of \(0\) and variance \(9.25\). The arithmetic is correct for the sampled, weighted points; the pilot has failed to represent the full target. A broader or defensive proposal, multiple starting regions, and independent checks for missed modes are more relevant than further tightening this update.

After adaptation, inspect more than the updated moments. Check the effective sample size and the largest normalized weight in fresh draws; compare estimates across independent batches; examine whether important regions are covered; and verify that the covariance remains well-conditioned and the update has not repeatedly hit its safeguards. A rising effective sample size is encouraging, but it is not proof that all modes were found or that the marginal likelihood estimate is accurate.

For a final marginal likelihood estimate, the clearest practice is to freeze the adapted proposal and draw an independent production sample. Conditional on the pilot and hence on the fixed proposal, the ordinary importance-sampling average from this new sample estimates \(Z\) without using those same draws to choose the proposal. Reusing adaptation draws in an estimate can introduce dependence between the selected proposal and the resulting estimate. More elaborate adaptive schemes can account for this dependence, but simple pilot-then-production sampling keeps the interpretation direct.

How This Approach Fits With Alternatives

A fixed Gaussian proposal centered using conditional modes can be effective when the mode and local curvature describe the important target region. Weighted moment adaptation instead uses sampled target information to revise both center and spread; it can respond to skewed placement or correlations that an initial approximation misses. Neither approach ensures good coverage of separated modes. A defensive mixture can protect against tails or regions outside a narrow Gaussian, but it may spend draws broadly rather than concentrating them. These methods can be combined, provided the proposal density used in each importance weight is the actual density from which the draw was generated.

Moment adaptation also has limits even when it converges: matching two moments does not make a Gaussian target-shaped, and finite pilot batches add Monte Carlo noise to both center and covariance. In NLME calculations, the target integrand can change with the individual and with population-parameter values. An adapted proposal for one individual or parameter setting should not be assumed adequate for another without checking.

Common Mistakes

  • Using unweighted pilot moments. The proposal draws describe \(q_t\), not the normalized integrand. Use weights \(g(z_j)/q_t(z_j)\) when estimating target moments.
  • Forgetting the proposal density in the weight. A large value of \(g(z_j)\) alone is not the importance weight; the density \(q_t(z_j)\) must also be included.
  • Treating a covariance update as automatically valid. Check symmetry and positive definiteness. A singular estimate requires regularization before it defines a Gaussian proposal.
  • Interpreting high effective sample size as proof of correctness. A proposal can give balanced weights among sampled points while missing an entire important region.
  • Assuming adaptation changes the integral. It changes where draws are placed. The target \(g\) and its integral \(Z\) remain the same.
Key takeaway: Weighted moment updates offer a practical way to move a Gaussian proposal toward the normalized integrand. Use normalized importance weights, damp and regularize unstable updates, and validate the adapted proposal with fresh draws and checks for missed regions before relying on a likelihood estimate.

Check Your Understanding

Answer each question, then compare your reasoning with the solutions.

  1. Why are weighted pilot moments preferable to the ordinary mean and covariance of the proposal draws?
  2. For normalized weights \((0.2,0.3,0.5)\), calculate the weighted mean of points \((0,1,3)\).
  3. What does damping with \(\lambda=0.25\) do to an update, and why might it help?
  4. Why can a high effective sample size fail to reveal a missed mode?

Solutions. (1) Proposal draws reflect \(q\), whereas weighting by \(g/q\) makes their normalized averages estimate moments of the target integrand's normalized distribution. (2) The mean is \(0.2(0)+0.3(1)+0.5(3)=1.8\). (3) It moves one quarter of the way from the current proposal parameters toward the weighted estimates, retaining three quarters of the current values; this can limit the effect of noisy pilot moments. (4) Effective sample size is calculated from weights of sampled points. If the proposal never samples a separate target mode, that mode contributes no weights and the diagnostic cannot detect its absence.