Introduction
In the previous page, we developed Monte Carlo methods for approximating posterior
summaries such as credible intervals and HPD regions, and introduced MCMC as a way to sample from posteriors that cannot be evaluated
in closed form. These methods all work with samples drawn from the target distribution itself. But what if sampling from the target is
prohibitively expensive, or we want to estimate the probability of a rare event in the distribution's tail?
Importance sampling addresses precisely these challenges by sampling from a different, more convenient distribution
and correcting for the mismatch.
Suppose we want to compute the expectation of a target function
\(\varphi(\mathbf{x})\) under a target distribution \(\pi(\mathbf{x})\):
\[
I = \mathbb{E}_{\pi}[\varphi(\mathbf{x})] = \int \varphi(\mathbf{x}) \, \pi(\mathbf{x}) \, d\mathbf{x}.
\]
Convention. Throughout this page, \(\pi\) denotes the target distribution and
\(q\) the proposal distribution. The policy notation \(\pi_\theta(a \mid s)\) appearing in the
RLHF example below is a separate, conventional reuse of the symbol \(\pi\) in reinforcement learning.
Standard Monte Carlo estimation would draw independent samples
\(\mathbf{x}_1, \ldots, \mathbf{x}_{N_s} \sim \pi(\mathbf{x})\) and
approximate:
\[
I \approx \hat{I}_{\text{MC}} = \frac{1}{N_s} \sum_{n=1}^{N_s} \varphi(\mathbf{x}_n).
\]
However, this approach fails when sampling from \(\pi(\mathbf{x})\) is computationally expensive or infeasible (as is common in
Bayesian inference), or when \(\varphi(\mathbf{x})\) has large values only in regions
where \(\pi(\mathbf{x})\) assigns low probability (rare events).
Importance sampling addresses these problems by introducing a proposal distribution (or
importance distribution) \(q(\mathbf{x})\) from which we can easily sample, then reweighting samples to account for
the distribution mismatch. In its self-normalized form, it requires \(\pi(\mathbf{x})\) only up to its normalizing constant, which is
often unknown in Bayesian inference. The idea of reusing samples collected under one distribution to estimate expectations under
another is one of the most broadly applicable principles in computational statistics.
Importance sampling is used throughout machine learning, in particular in the following domains:
Large language models (LLMs) fine-tuned by PPO-based RLHF rely on importance sampling. The
Proximal Policy Optimization (PPO) algorithm uses importance sampling to enable multiple gradient updates from the
same batch of collected data:
\[
L(\theta) = \mathbb{E}_{(s,a) \sim \mu}\left[\min\left(\frac{\pi_\theta(a \mid s)}{\mu(a \mid s)} A(s,a), \operatorname{clip}\left(\frac{\pi_\theta(a \mid s)}{\mu(a \mid s)}, 1-\epsilon, 1+\epsilon\right) A(s,a)\right)\right]
\]
where \(\mu\) is the behavior policy (data collection policy), \(\pi_\theta\) is the current policy being optimized, and \(A(s,a)\) is
an estimate of the advantage of action \(a\) in state \(s\).
The importance ratio \(\pi_\theta/\mu\) allows us to reuse data collected under \(\mu\) for multiple updates to \(\pi_\theta\). PPO's
clipping mechanism (typically with \(\epsilon = 0.2\)) removes the incentive to move this ratio outside
\([1 - \epsilon, 1 + \epsilon]\) and thereby stabilizes training. Without importance sampling, every gradient step would require
fresh on-policy data.
2. Off-Policy Reinforcement Learning
In robotics, autonomous systems, and game AI, collecting on-policy data can be dangerous or expensive. Importance sampling enables
learning from historical data, human demonstrations, or safer exploration policies. It underlies off-policy policy evaluation and
off-policy policy-gradient methods.
3. Variational Inference and Deep Generative Models
Importance-weighted autoencoders (IWAE) use importance sampling to obtain tighter lower bounds on the log-evidence than the standard
evidence lower bound (ELBO). The tighter bound can
yield better generative models and more accurate uncertainty estimates in Bayesian deep learning.
4. Imbalanced Data and Rare Events
In applications such as fraud detection, medical diagnosis of rare diseases, and AI safety, the events we care most about are extremely
rare. Importance sampling allows us to oversample these critical cases and reweight appropriately, so that models do not ignore rare
but important scenarios.
5. Computing Normalizing Constants
Many probabilistic models in AI (Boltzmann machines, energy-based models, certain Bayesian posteriors) have intractable
normalizing constants. Importance sampling, particularly annealed importance sampling, provides practical methods to
estimate these constants, which is essential for likelihood-based model comparison.
Direct Importance Sampling
The key mathematical insight behind importance sampling is the change of measure. Suppose \(q(\mathbf{x})\) is a proposal distribution
satisfying the support condition: \(q(\mathbf{x}) \gt 0\) whenever \(\pi(\mathbf{x}) \gt 0\). Then we can rewrite the
expectation under \(\pi(\mathbf{x})\) as an expectation under \(q(\mathbf{x})\):
\[
\begin{align*}
I &= \int \varphi(\mathbf{x}) \, \pi(\mathbf{x}) \, d\mathbf{x} \\\\
&= \int_{\{\pi \gt 0\}} \varphi(\mathbf{x}) \, \frac{\pi(\mathbf{x})}{q(\mathbf{x})} \, q(\mathbf{x}) \, d\mathbf{x} \\\\
&= \mathbb{E}_{q}\!\left[\varphi(\mathbf{x}) \, \frac{\pi(\mathbf{x})}{q(\mathbf{x})}\right].
\end{align*}
\]
The second equality uses the support condition, which guarantees that the ratio \(\pi/q\) is well-defined on the integration domain.
The third equality extends the integration to all of \(\{q \gt 0\}\) by setting the integrand to zero where \(\pi = 0\). This
redefinition does not change the integral, since \(\varphi \cdot \pi = 0\) on that set.
Definition: Importance Weights
Given a target distribution \(\pi(\mathbf{x})\) and a proposal distribution \(q(\mathbf{x})\) with \(q(\mathbf{x}) \gt 0\) whenever
\(\pi(\mathbf{x}) \gt 0\), the importance weight function is defined as
\[
w(\mathbf{x}) = \frac{\pi(\mathbf{x})}{q(\mathbf{x})}.
\]
For a sample \(\mathbf{x}_n \sim q\), the weight
\[
\tilde{w}_n := w(\mathbf{x}_n) = \pi(\mathbf{x}_n)/q(\mathbf{x}_n)
\]
measures how much more (or less) likely \(\mathbf{x}_n\) is under the target compared to the proposal. Note that while \(\pi\) and
\(q\) are normalized distributions, the weights \(\tilde{w}_n\) do not sum to any particular value.
Assume that the normalized \(\pi(\mathbf{x})\) can be evaluated pointwise, but that we cannot sample from it.
If we draw \(N_s\) samples \(\mathbf{x}_1, \ldots, \mathbf{x}_{N_s} \sim q(\mathbf{x})\),
the direct importance sampling estimator is:
\[
\hat{I}_{\text{IS}} := \frac{1}{N_s} \sum_{n=1}^{N_s} \tilde{w}_n \, \varphi(\mathbf{x}_n).
\]
This estimator is unbiased whenever \(I\) exists. (Absolute integrability of \(\varphi\) under \(\pi\) is the same as
absolute integrability of \(w\varphi\) under \(q\), by the change of measure above.) Since \(\mathbf{x}_1, \ldots, \mathbf{x}_{N_s}\)
are i.i.d. samples from \(q\),
linearity of expectation (iterated over
the \(N_s\) terms) gives
\[
\begin{align*}
\mathbb{E}_q[\hat{I}_{\text{IS}}] &= \frac{1}{N_s}\sum_{n=1}^{N_s} \mathbb{E}_q\!\left[w(\mathbf{x}_n)\,\varphi(\mathbf{x}_n)\right] \\\\
&= \mathbb{E}_q\!\left[w(\mathbf{X})\,\varphi(\mathbf{X})\right] \\\\
&= I,
\end{align*}
\]
where the second equality uses that all terms have the same distribution (we write \(\mathbf{X} \sim q\) for a single representative),
and the third invokes the change-of-measure identity above.
Theorem: Variance of the Direct IS Estimator
Under the support condition, with \(\mathbf{x}_1, \ldots, \mathbf{x}_{N_s} \overset{\text{i.i.d.}}{\sim} q\),
let \(\mathbf{X} \sim q\) denote a single representative sample. If
\(\mathbb{E}_q[w(\mathbf{X})^2 \varphi(\mathbf{X})^2] \lt \infty\), then
\[
\operatorname{Var}_q[\hat{I}_{\text{IS}}] = \frac{1}{N_s}\!\left(\mathbb{E}_q[w(\mathbf{X})^2 \varphi(\mathbf{X})^2] - I^2\right).
\]
Proof.
Let \(Y_n := w(\mathbf{x}_n)\varphi(\mathbf{x}_n)\). Since the \(\mathbf{x}_n\) are i.i.d. under \(q\), so are the \(Y_n\) (we
take for granted that measurable functions of i.i.d. samples are again i.i.d.). Hence
\[
\begin{align*}
\operatorname{Var}_q[\hat{I}_{\text{IS}}] &= \operatorname{Var}_q\!\left[\frac{1}{N_s} \sum_{n=1}^{N_s} Y_n\right] \\\\
&= \frac{1}{N_s^2} \sum_{n=1}^{N_s} \operatorname{Var}_q[Y_n] \\\\
&= \frac{1}{N_s} \operatorname{Var}_q[Y_1],
\end{align*}
\]
where the third equality uses that the \(Y_n\) are identically distributed. For the second, each \(Y_n\) has mean \(I\) and
finite second moment \(\mathbb{E}_q[Y_n^2] = \mathbb{E}_q[w(\mathbf{X})^2 \varphi(\mathbf{X})^2]\) by hypothesis. Since
\(\sum_n Y_n - N_s I = \sum_n (Y_n - I)\), expanding the square and using linearity of expectation over the finitely many terms
(each finite, since \(|ab| \leq (a^2 + b^2)/2\)) gives
\[
\begin{align*}
\operatorname{Var}_q\!\left[\sum_{n=1}^{N_s} Y_n\right] &= \sum_{n=1}^{N_s} \mathbb{E}_q\bigl[(Y_n - I)^2\bigr] \\\\
&\quad + \sum_{n \neq m} \mathbb{E}_q\bigl[(Y_n - I)(Y_m - I)\bigr].
\end{align*}
\]
Each diagonal term is \(\operatorname{Var}_q[Y_n]\) by definition, and each cross term is \(\operatorname{Cov}_q[Y_n, Y_m]\),
which vanishes because distinct \(Y_n\) and \(Y_m\) are independent (part of the i.i.d. property just taken for granted) and
independence implies zero covariance.
That corollary is proved there only for pairs with a joint mass function or density, which \(Y_n\) and \(Y_m\) need not have, and
as with the variance identity below we take for granted that it holds for any independent pair with finite second moments. Only
this pairwise independence is used. The sum therefore has finite variance \(\sum_n \operatorname{Var}_q[Y_n]\), and the
variance of a linear transformation
by \(1/N_s\) supplies the factor \(1/N_s^2\). The single-sample variance follows from the
computational identity for variance
\(\operatorname{Var}_q[Y_1] = \mathbb{E}_q[Y_1^2] - (\mathbb{E}_q[Y_1])^2\) together with unbiasedness. The identity is stated
there for a random variable with a mass function or density, and we take for granted that it holds for any random variable with
finite second moment, such as \(Y_1\), which has the distribution of \(w(\mathbf{X})\varphi(\mathbf{X})\):
\[
\begin{align*}
\operatorname{Var}_q[Y_1] &= \mathbb{E}_q[w(\mathbf{X})^2 \varphi(\mathbf{X})^2] - \left(\mathbb{E}_q[w(\mathbf{X})\varphi(\mathbf{X})]\right)^2 \\\\
&= \mathbb{E}_q[w(\mathbf{X})^2 \varphi(\mathbf{X})^2] - I^2.
\end{align*}
\]
If the importance weights vary dramatically, a few samples often carry most of the total weight, and unless \(\varphi\) offsets the
weights, the variance is typically high. This degeneracy is quantified by the effective sample size (ESS), a
heuristic computed from the weights alone for the number of independent samples from \(\pi\) that would give comparable statistical
efficiency. The ESS can be computed in two equivalent ways:
Definition: Effective Sample Size
The effective sample size (ESS) of a weighted sample \(\{(\mathbf{x}_n, \tilde{w}_n)\}_{n=1}^{N_s}\) with
nonnegative weights, at least one of them positive, is defined as
\[
\begin{align*}
\operatorname{ESS} &:= \frac{\left(\sum_{n=1}^{N_s} \tilde{w}_n\right)^2}{\sum_{n=1}^{N_s} \tilde{w}_n^2} \\\\
&= \frac{1}{\sum_{n=1}^{N_s} W_n^2},
\end{align*}
\]
where \(W_n := \tilde{w}_n / \sum_{n'=1}^{N_s} \tilde{w}_{n'}\) are the normalized weights. The second equality follows by dividing
numerator and denominator of the first expression by \(\left(\sum_n \tilde{w}_n\right)^2\).
The ESS ranges from 1 (maximum degeneracy, where all weight concentrates on a single sample) to \(N_s\) (no degeneracy, where
all weights are equal). A good rule of thumb is that if \(\operatorname{ESS} \lt N_s/2\), the proposal distribution
\(q(\mathbf{x})\) should be reconsidered, as the normalized weights then satisfy \(\sum_n W_n^2 \gt 2/N_s\), twice the value
\(1/N_s\) attained by equal weights.
Requirements for Direct Importance Sampling
For direct importance sampling to be valid and practical:
- Support condition. \(q(\mathbf{x}) \gt 0\) whenever \(\pi(\mathbf{x}) \gt 0\)
- Normalization. Both \(\pi(\mathbf{x})\) and \(q(\mathbf{x})\) must be properly normalized probability distributions
- Evaluability. We must be able to evaluate \(\pi(\mathbf{x})\) and \(q(\mathbf{x})\) for any \(\mathbf{x}\)
- Sampling. We must be able to draw samples from \(q(\mathbf{x})\) efficiently
Self-Normalized Importance Sampling
Direct importance sampling, as developed above, requires evaluating the normalized target density \(\pi(\mathbf{x})\) pointwise.
In many real-world settings this is precisely what we cannot do. The next method removes this requirement by working with
unnormalized densities throughout.
Often we can evaluate only the unnormalized target distribution:
\[
\tilde{\gamma}(\mathbf{x}) = Z \pi(\mathbf{x})
\]
where the normalization constant
\[
Z = \int \tilde{\gamma}(\mathbf{x})d\mathbf{x}
\]
is unknown or intractable. This situation is ubiquitous in Bayesian inference (where \(Z\) is the marginal likelihood), energy-based
models (where \(Z\) is the partition function), and many other machine learning applications.
Self-normalized importance sampling (SNIS) addresses this by working with unnormalized weights and normalizing them
within the estimator itself. The key identity rewrites the expectation as a ratio whose numerator and denominator can both be
approximated by Monte Carlo:
\[
\begin{align*}
\mathbb{E}_\pi[\varphi(\mathbf{x})] &= \int \varphi(\mathbf{x}) \, \pi(\mathbf{x}) \, d\mathbf{x} \\\\
&= \frac{\int \varphi(\mathbf{x}) \, \tilde{\gamma}(\mathbf{x}) \, d\mathbf{x}}
{\int \tilde{\gamma}(\mathbf{x}) \, d\mathbf{x}} \\\\
&= \frac{\mathbb{E}_q\!\left[\varphi(\mathbf{x}) \, \tilde{\gamma}(\mathbf{x})/q(\mathbf{x})\right]}
{\mathbb{E}_q\!\left[\tilde{\gamma}(\mathbf{x})/q(\mathbf{x})\right]}.
\end{align*}
\]
The second equality uses \(\tilde{\gamma} = Z\pi\). The factor \(Z\) appears in both the numerator and the denominator (since
\(\int \tilde{\gamma} \, d\mathbf{x} = Z\)) and therefore cancels. The third equality applies the change-of-measure identity to both
integrals and therefore requires the support condition \(q(\mathbf{x}) \gt 0\) whenever \(\tilde{\gamma}(\mathbf{x}) \gt 0\).
Sampling \(\mathbf{x}_n \sim q(\mathbf{x})\) for \(n = 1, \ldots, N_s\), we redefine the importance weights using
the unnormalized target:
\[
\tilde{w}_n := \frac{\tilde{\gamma}(\mathbf{x}_n)}{q(\mathbf{x}_n)}.
\]
These weights differ from the importance weights defined earlier by a factor of \(Z\). The SNIS estimator below is invariant under this
rescaling because the factor cancels between numerator and denominator.
Definition: Self-Normalized Importance Sampling
The self-normalized importance sampling (SNIS) estimator is
\[
\begin{align*}
\hat{I}_{\text{SNIS}} &:= \frac{\sum_{n=1}^{N_s} \tilde{w}_n \, \varphi(\mathbf{x}_n)}{\sum_{n=1}^{N_s} \tilde{w}_n} \\\\
&= \sum_{n=1}^{N_s} W_n \, \varphi(\mathbf{x}_n),
\end{align*}
\]
where \(\tilde{w}_n = \tilde{\gamma}(\mathbf{x}_n)/q(\mathbf{x}_n)\) are the unnormalized importance weights and
\(W_n := \tilde{w}_n / \sum_{n'=1}^{N_s} \tilde{w}_{n'}\) are the normalized weights.
This estimator is in general biased for finite \(N_s\) (the expectation of a ratio is not generally the ratio of
expectations), but it is consistent:
Theorem: Consistency of the SNIS Estimator
Assume the support condition \(q(\mathbf{x}) \gt 0\) whenever \(\tilde{\gamma}(\mathbf{x}) \gt 0\),
\(Z := \int \tilde{\gamma}(\mathbf{x}) \, d\mathbf{x}\) satisfies \(0 \lt Z \lt \infty\), and
\(\mathbb{E}_q[|w(\mathbf{X}) \varphi(\mathbf{X})|] \lt \infty\), where
\(w(\mathbf{x}) := \tilde{\gamma}(\mathbf{x})/q(\mathbf{x})\). Then
\[
\hat{I}_{\text{SNIS}} \xrightarrow{a.s.} \mathbb{E}_\pi[\varphi(\mathbf{X})] = I \quad \text{as } N_s \to \infty.
\]
Proof.
Write \(\hat{I}_{\text{SNIS}} = A_{N_s} / B_{N_s}\), where
\[
\begin{align*}
A_{N_s} &:= \frac{1}{N_s}\sum_{n=1}^{N_s} \tilde{w}_n \, \varphi(\mathbf{x}_n), \\\\
B_{N_s} &:= \frac{1}{N_s}\sum_{n=1}^{N_s} \tilde{w}_n.
\end{align*}
\]
(Multiplying numerator and denominator of \(\hat{I}_{\text{SNIS}}\) by \(1/N_s\) preserves the ratio.)
Because the \(\mathbf{x}_n\) are i.i.d. under \(q\), the summands of \(A_{N_s}\) and of \(B_{N_s}\) are again i.i.d. (as in the
variance proof above), and their means are finite because \(\mathbb{E}_q[|w(\mathbf{X}) \varphi(\mathbf{X})|] \lt \infty\) by
assumption and \(\mathbb{E}_q[w(\mathbf{X})] = Z \lt \infty\). The
Strong Law of Large Numbers therefore gives
\[
\begin{align*}
A_{N_s} &\xrightarrow{a.s.} \mathbb{E}_q\!\left[w(\mathbf{X}) \, \varphi(\mathbf{X})\right], \\\\
B_{N_s} &\xrightarrow{a.s.} \mathbb{E}_q\!\left[w(\mathbf{X})\right].
\end{align*}
\]
Both expectations evaluate explicitly via change of measure:
\[
\begin{align*}
\mathbb{E}_q[w(\mathbf{X})] &= \int \tilde{\gamma}(\mathbf{x}) \, d\mathbf{x} = Z, \\\\
\mathbb{E}_q[w(\mathbf{X}) \, \varphi(\mathbf{X})] &= \int \varphi(\mathbf{x}) \, \tilde{\gamma}(\mathbf{x}) \, d\mathbf{x} = Z \cdot I,
\end{align*}
\]
where the second uses \(\tilde{\gamma} = Z\pi\) and \(I = \int \varphi \, \pi \, d\mathbf{x}\).
The two convergences hold on events of probability one, and the intersection of two such events again has probability one,
so \((A_{N_s}, B_{N_s}) \xrightarrow{a.s.} (Z I, Z)\). With \(Z \gt 0\), the function \((a, b) \mapsto a/b\) is continuous
at \((Z I, Z)\), so on that intersection \(A_{N_s}/B_{N_s} \to I\) pointwise. This is the almost-sure, vector-valued form of
the continuous mapping theorem (the
linked statement covers the scalar case in probability, and the pointwise argument just given is all that the almost-sure
case needs). Hence
\[
\hat{I}_{\text{SNIS}} = \frac{A_{N_s}}{B_{N_s}} \xrightarrow{a.s.} \frac{Z \cdot I}{Z} = I.
\]
The estimator converges almost surely to
the true expectation as \(N_s \to \infty\). The bias that SNIS can have at fixed \(N_s\) does not prevent this convergence.
Despite the bias, SNIS is widely used in practice. It requires evaluating only the unnormalized densities
\(\tilde{\gamma}(\mathbf{x})\) and \(q(\mathbf{x})\), so it applies whenever the normalization constant \(Z\) is intractable. SNIS
can reduce the dispersion (variance) of the estimator in certain cases, but it does not universally reduce variance compared to
direct importance sampling.
Annealed Importance Sampling (AIS)
Self-normalized importance sampling resolves the difficulty of unknown normalization, but it does not address a different failure mode.
When the target and proposal distributions have a substantial mismatch, importance sampling, direct or self-normalized, fails because
the importance weights have extremely high variance. Most samples from the proposal will have near-zero weight under the target, while
a few rare samples will have enormous weights. This phenomenon is known as weight degeneracy. The effective sample
size is then low. Despite the many samples drawn, few contribute meaningfully to the estimate, and the results are unreliable.
Annealed importance sampling (AIS) mitigates this limitation by introducing a sequence of intermediate distributions
that smoothly bridge the gap between the proposal and target distributions. The importance weights constructed along this sequence
are typically far more stable.
Suppose we wish to sample from a complicated target distribution:
\[
p_0(\mathbf{x}) \propto f_0(\mathbf{x}),
\]
where sampling is difficult because \(p_0\) may be high-dimensional or multimodal. Here \(f_0\) is the unnormalized target (the same
role played by \(\tilde{\gamma}\) in the previous section). We adopt the AIS-specific notation \(f_0, f_n\) to emphasize the indexed
family of intermediate distributions introduced below.
Suppose also that we have an easily sampled proposal distribution (for example, the prior in a Bayesian setting):
\[
p_n(\mathbf{x}) \propto f_n(\mathbf{x}).
\]
AIS constructs a sequence of \(n+1\) distributions \(p_j \propto f_j\), \(j = 0, 1, \ldots, n\), interpolating between
\(p_n\) and \(p_0\):
\[
f_j(\mathbf{x}) = f_0(\mathbf{x})^{\beta_j} f_n(\mathbf{x})^{1 - \beta_j},
\]
where \(1 = \beta_0 \gt \beta_1 \gt \cdots \gt \beta_n = 0\) defines an annealing schedule. At the
endpoints, \(\beta_0 = 1\) gives \(f_0\) (the target) and \(\beta_n = 0\) gives \(f_n\) (the proposal). Intermediate values
produce geometric mixtures.
Each \(\beta_j\) acts as an inverse temperature in the statistical-physics sense. Small \(\beta_j\) corresponds to
"high temperature" (closer to the easy proposal), and large \(\beta_j\) to "low temperature" (closer to the complicated target). The
AIS sampling procedure traverses these distributions in reverse, starting from \(p_n\) and ending at \(p_0\). Throughout this section
we assume that every \(f_j\) is strictly positive on a common domain in which all the kernels below move, so that every ratio of the
\(f_j\) that appears is well-defined. This holds, for instance, when \(f_0\) and \(f_n\) are both positive everywhere.
Typical schedules include:
- Linear. \(\beta_j = \frac{n - j}{n}\)
- Power. \(\beta_j = \left( \frac{n - j}{n} \right)^{\alpha}\), for some \(\alpha \gt 0\)
- Adaptive. Chosen to maintain approximately constant effective sample size (ESS) between steps
To sample from each intermediate distribution \(p_j\), we use a
Markov chain transition kernel \(T_j(\mathbf{x}, \mathbf{x}')\), the density of the next
state \(\mathbf{x}'\) given the current state \(\mathbf{x}\), which leaves \(p_j\) invariant, that is,
\[
\int p_j(\mathbf{x})T_j(\mathbf{x}, \mathbf{x}')d\mathbf{x} = p_j(\mathbf{x}').
\]
Common choices for \(T_j\) include Metropolis-Hastings or Hamiltonian Monte Carlo, applied for a few iterations to adapt
samples toward \(p_j\).
The AIS procedure proceeds as follows:
- Sample \(\mathbf{v}_n \sim p_n\) from the proposal distribution.
- Apply successive Markov transitions to obtain \(\mathbf{v}_0\):
\[
\begin{align*}
\mathbf{v}_{n-1} &\sim T_{n-1}(\mathbf{v}_n, \cdot), \\\\
\mathbf{v}_{n-2} &\sim T_{n-2}(\mathbf{v}_{n-1}, \cdot), \quad \ldots, \\\\
\mathbf{v}_0 &\sim T_0(\mathbf{v}_1, \cdot).
\end{align*}
\]
Each AIS run yields a trajectory \((\mathbf{v}_n, \mathbf{v}_{n-1}, \ldots, \mathbf{v}_0)\) with an associated importance weight:
\[
w =
\frac{f_{n-1}(\mathbf{v}_n)}{f_n(\mathbf{v}_n)}
\frac{f_{n-2}(\mathbf{v}_{n-1})}{f_{n-1}(\mathbf{v}_{n-1})}
\cdots
\frac{f_1(\mathbf{v}_2)}{f_2(\mathbf{v}_2)}
\frac{f_0(\mathbf{v}_1)}{f_1(\mathbf{v}_1)}. \tag{1}
\]
Each ratio is close to 1 when successive \(\beta_j\) are close, and with well-mixing transitions \(T_j\) this keeps the weights more
stable than in one-step importance sampling. The following theorem justifies calling \(w\) an importance weight.
Theorem: Validity of the AIS Weights
Assume, as above, that every \(f_j\) is strictly positive on the common domain, and in addition that every \(f_j\) has a finite
integral. Let each \(T_j\), \(j = 0, \ldots, n-1\), be a transition density \(T_j(\mathbf{x}, \mathbf{x}')\) as above that leaves
\(p_j\) invariant, with the invariance identity holding for every \(\mathbf{x}'\), and regard the weight \(w\) of Equation (1) as
a function of the trajectory \(\mathbf{v} = (\mathbf{v}_0, \ldots, \mathbf{v}_n)\). Then there is a nonnegative function
\(\tilde{f}\) on the extended space of trajectories with the following properties. (a)
\(\int \tilde{f}(\mathbf{v}) \, d\mathbf{v}_1 \cdots d\mathbf{v}_n = f_0(\mathbf{v}_0)\), so the \(\mathbf{v}_0\)-marginal of the
extended target \(p \propto \tilde{f}\) is \(p_0\). (b) For every \(\mathbf{v}\),
\[
\tilde{f}(\mathbf{v}) = w \, f_n(\mathbf{v}_n) \, T_{n-1}(\mathbf{v}_n, \mathbf{v}_{n-1}) \cdots T_0(\mathbf{v}_1, \mathbf{v}_0),
\]
where the product after \(w\) is proportional to the joint density of the trajectory generated by the AIS procedure, so \(w\) is
an importance weight of \(p\) against that trajectory distribution. (c) If
\(\mathbb{E}_{p_0}[|\varphi(\mathbf{X})|] \lt \infty\), the SNIS estimator built from \(S\) independent AIS runs, with values
\(\varphi(\mathbf{v}_0)\) and weights \(w\), converges almost surely to \(\mathbb{E}_{p_0}[\varphi(\mathbf{X})]\) as
\(S \to \infty\).
Proof.
The key technical device is to define the extended target distribution \(\tilde{f}\) on the product space of
trajectories as a chain that starts from \(f_0\) and runs backward through reversals of the kernels \(T_j\).
Consider the following unnormalized density on the extended space:
\[
\begin{align*}
p(\mathbf{v}) \propto \tilde{f}(\mathbf{v})
&= f_0(\mathbf{v}_0) \, \tilde{T}_0(\mathbf{v}_0, \mathbf{v}_1) \, \tilde{T}_1(\mathbf{v}_1, \mathbf{v}_2) \cdots \tilde{T}_{n-1}(\mathbf{v}_{n-1}, \mathbf{v}_n) \\\\
&\propto p(\mathbf{v}_0) \, p(\mathbf{v}_1 \mid \mathbf{v}_0) \cdots p(\mathbf{v}_n \mid \mathbf{v}_{n-1}),
\end{align*}
\]
where \(\tilde{T}_j\) is the reversal of \(T_j\) with respect to \(p_j\), defined by
\[
\begin{align*}
\tilde{T}_j(\mathbf{v}, \mathbf{v}') &= T_j(\mathbf{v}', \mathbf{v}) \, p_j(\mathbf{v}') \, / \, p_j(\mathbf{v}) \\\\
&= T_j(\mathbf{v}', \mathbf{v}) \, f_j(\mathbf{v}') \, / \, f_j(\mathbf{v}).
\end{align*}
\]
The second equality holds because the normalization constant \(Z_j\) of \(p_j = f_j/Z_j\) cancels in the ratio.
We write \(\tilde{f}(\mathbf{v})\) for the unnormalized extended density to distinguish it from the unnormalized densities
\(f_0, \ldots, f_n\) on the original space.
The invariance of \(p_j\) with respect to \(T_j\) ensures that each \(\tilde{T}_j\) is a valid Markov
transition kernel. Explicitly,
\[
\begin{align*}
\int \tilde{T}_j(\mathbf{v}, \mathbf{v}') \, d\mathbf{v}' &= \frac{1}{f_j(\mathbf{v})} \int T_j(\mathbf{v}', \mathbf{v}) \, f_j(\mathbf{v}') \, d\mathbf{v}' \\\\
&= \frac{f_j(\mathbf{v})}{f_j(\mathbf{v})} = 1,
\end{align*}
\]
where the second equality uses the invariance condition
\(\int T_j(\mathbf{v}', \mathbf{v}) \, p_j(\mathbf{v}') \, d\mathbf{v}' = p_j(\mathbf{v})\), rescaled by \(Z_j\)
to remove normalization.
For (a), that is, to see that the \(\mathbf{v}_0\)-marginal of \(p(\mathbf{v})\) is \(p_0\), we integrate out
\(\mathbf{v}_1, \ldots, \mathbf{v}_n\):
\[
\int \tilde{f}(\mathbf{v}) \, d\mathbf{v}_1 \cdots d\mathbf{v}_n = f_0(\mathbf{v}_0),
\]
where successive applications of \(\int \tilde{T}_j(\cdot, \mathbf{v}_{j+1}) \, d\mathbf{v}_{j+1} = 1\) collapse the chain from
the right. Hence sampling from \(p(\mathbf{v})\) yields, as a byproduct, samples from \(p_0(\mathbf{x})\) via
\(\mathbf{x} := \mathbf{v}_0\).
The AIS procedure samples on this extended state space using a forward chain starting from \(p_n\). This corresponds
to the proposal:
\[
\begin{align*}
q(\mathbf{v}) \propto g(\mathbf{v}) &= f_n(\mathbf{v}_n) \, T_{n-1}(\mathbf{v}_n, \mathbf{v}_{n-1}) \cdots T_1(\mathbf{v}_2, \mathbf{v}_1) \, T_0(\mathbf{v}_1, \mathbf{v}_0) \\\\
&\propto q(\mathbf{v}_n) \, q(\mathbf{v}_{n-1} \mid \mathbf{v}_n) \cdots q(\mathbf{v}_0 \mid \mathbf{v}_1).
\end{align*}
\]
The importance weight on the extended space is then
\[
w = \frac{\tilde{f}(\mathbf{v}_0, \ldots, \mathbf{v}_n)}{g(\mathbf{v}_0, \ldots, \mathbf{v}_n)}.
\]
Substituting the definition of \(\tilde{T}_j\) into \(\tilde{f}\) gives
\[
\tilde{f}(\mathbf{v}) = f_0(\mathbf{v}_0) \prod_{j=0}^{n-1} T_j(\mathbf{v}_{j+1}, \mathbf{v}_j) \, \frac{f_j(\mathbf{v}_{j+1})}{f_j(\mathbf{v}_j)}.
\]
Every forward kernel \(T_j(\mathbf{v}_{j+1}, \mathbf{v}_j)\) also appears in \(g\). Dividing by \(g\) therefore cancels the
kernels (first line below), and \(f_0(\mathbf{v}_0)\) cancels against the \(j = 0\) denominator (second line):
\[
\begin{align*}
w &= \frac{f_0(\mathbf{v}_0)}{f_n(\mathbf{v}_n)} \prod_{j=0}^{n-1} \frac{f_j(\mathbf{v}_{j+1})}{f_j(\mathbf{v}_j)} \\\\
&= \prod_{j=1}^{n} \frac{f_{j-1}(\mathbf{v}_j)}{f_j(\mathbf{v}_j)},
\end{align*}
\]
which is exactly Equation (1). Multiplying back, \(\tilde{f}(\mathbf{v}) = w \, g(\mathbf{v})\) holds for every \(\mathbf{v}\),
including where \(g(\mathbf{v}) = 0\), since the kernels enter \(\tilde{f}\) and \(g\) as the same factors; this is (b). The
weight does not involve \(\mathbf{v}_0\), so the final transition \(T_0\) leaves it unchanged.
For (c), integrating \(g\) successively over \(\mathbf{v}_0, \ldots, \mathbf{v}_{n-1}\) gives \(f_n(\mathbf{v}_n)\), hence
\(\int g(\mathbf{v}) \, d\mathbf{v} = Z_n\), so we may apply the SNIS consistency theorem on the extended space with
\(\tilde{\gamma} = \tilde{f}\) and \(q = g / Z_n\). Its support condition holds by (b), (a) gives
\(\int \tilde{f}(\mathbf{v}) \, d\mathbf{v} = Z_0 \in (0, \infty)\), and (a) also gives
\(\mathbb{E}_q[(\tilde{f}/q) \, |\varphi(\mathbf{v}_0)|] = \int |\varphi(\mathbf{x})| \, f_0(\mathbf{x}) \, d\mathbf{x}\), which
equals \(Z_0 \, \mathbb{E}_{p_0}[|\varphi(\mathbf{X})|] \lt \infty\). The SNIS weights \(\tilde{f}/q = Z_n w\) differ from \(w\)
by a constant factor, which cancels in the estimator.
A key application of AIS is to estimate the ratio of partition functions. Since part (a) of the theorem gives
\(\int \tilde{f}(\mathbf{v}) \, d\mathbf{v}_1 \cdots d\mathbf{v}_n = f_0(\mathbf{v}_0)\), integrating once more yields
\[
Z_0 = \int f_0(\mathbf{x}) \, d\mathbf{x} = \int \tilde{f}(\mathbf{v}) \, d\mathbf{v},
\]
and an analogous identity for the proposal side gives
\[
Z_n = \int f_n(\mathbf{x}) \, d\mathbf{x} = \int g(\mathbf{v}) \, d\mathbf{v}.
\]
Therefore:
\[
\begin{align*}
\frac{Z_0}{Z_n} &= \frac{\int \tilde{f}(\mathbf{v}) \, d\mathbf{v}}{\int g(\mathbf{v}) \, d\mathbf{v}} \\\\
&= \frac{\int \frac{\tilde{f}(\mathbf{v})}{g(\mathbf{v})} g(\mathbf{v}) \, d\mathbf{v}}{\int g(\mathbf{v}) \, d\mathbf{v}} \\\\
&= \mathbb{E}_g\!\left[\frac{\tilde{f}(\mathbf{v})}{g(\mathbf{v})}\right] \\\\
&\approx \frac{1}{S}\sum_{s=1}^S w_s,
\end{align*}
\]
where \(w_s = \tilde{f}(\mathbf{v}_s) / g(\mathbf{v}_s)\) is the AIS weight from Equation (1) for the \(s\)-th run, and \(S\) is the
number of independent AIS runs.
By abuse of notation, \(\mathbb{E}_g\) here denotes the expectation under the normalized distribution \(q \propto g\). The
normalization constant \(Z_n\) cancels between numerator and denominator in the ratio above, so the formula is unaffected by which
version of \(g\) is used.
A common application is Bayesian model evidence estimation. Here the proposal is the prior,
\(f_n(\boldsymbol{\theta}) = p(\boldsymbol{\theta})\) with known normalization \(Z_n = 1\), and the target is the unnormalized
posterior \(f_0(\boldsymbol{\theta}) = p(\boldsymbol{\theta}) \, p(\mathcal{D} \mid \boldsymbol{\theta})\), whose normalization
constant is the marginal likelihood \(Z_0 = p(\mathcal{D})\). The ratio formula above then yields
\[
p(\mathcal{D}) = Z_0 = \frac{Z_0}{Z_n} \approx \frac{1}{S}\sum_{s=1}^S w_s,
\]
so the same AIS runs whose samples, together with their weights \(w_s\), approximate the posterior also provide a direct estimator of
the marginal likelihood.
The variance of AIS can be reduced by increasing the number of intermediate distributions \(n\). In practice, \(n = 100\) to
\(n = 10000\) steps is common. The optimal choice balances computational cost (proportional to \(n\)) against variance reduction.
Several practical considerations affect the performance of AIS:
1. Choosing the Annealing Schedule
The choice of \(\{\beta_j\}\) is crucial. Several designs have been proposed:
- Moment matching. Choose \(\beta_j\) such that the second moments of successive distributions are approximately equal
- Constant ESS. Adapt \(\beta_j\) to maintain ESS \(\approx 0.8 \times S\) between steps
- Sigmoid schedules. \(\beta_j = \frac{\sigma(a(1 - 2j/n)) - \sigma(-a)}{\sigma(a) - \sigma(-a)}\) with \(\sigma\) the sigmoid and \(a \gt 0\) its steepness, rescaled to \(\beta_0 = 1\), \(\beta_n = 0\)
2. Transition Kernels
The choice and tuning of \(T_j\) significantly affects efficiency:
- Hamiltonian Monte Carlo (HMC). Often most effective for continuous distributions, especially in high dimensions
- Number of MCMC steps. Usually 1-10 steps per temperature is sufficient. More steps increase cost with diminishing returns
- Step size adaptation. The HMC step size is critical and can be adapted based on the acceptance rate at each temperature
3. Applications in Deep Learning
AIS is used in deep learning for:
- Evaluating generative models. Estimating log-likelihoods of VAEs and energy-based models
- Bayesian deep learning. Computing model evidence for neural network architectures
- Differentiable AIS. Gradient-based optimization through AIS estimates, using techniques such as the reparameterization trick
Diagnostics and Debugging
Importance sampling can fail silently, and the estimates it produces may have high bias or variance. Here are essential
diagnostics to monitor:
1. Weight Diagnostics
The following metrics monitor the health of the importance weights at a glance. The effective sample size quantifies overall weight
degeneracy, while the maximum normalized weight, coefficient of variation, and entropy capture complementary aspects of how
concentrated or dispersed the weight distribution is.
Key Diagnostic Metrics
- Effective Sample Size (ESS). Values below \(N_s/2\) call for a better proposal (rule of thumb). A large ESS does not by itself certify the estimate
- Maximum normalized weight. \(\max_n W_n\) should be small (typically \(\lt 0.1\) is recommended for \(N_s \gt 10\), since \(\max_n W_n \geq 1/N_s\)) to avoid single-sample domination
- Coefficient of variation. \(\operatorname{CV} = \frac{\operatorname{std}(\tilde{w})}{\operatorname{mean}(\tilde{w})}\) should be \(\lt 1\) (rule of thumb)
- Weight entropy. \(H(W) = -\sum_n W_n \log W_n\) (\(0 \log 0 := 0\)) measures weight uniformity (higher is better). Here \(\log\) is the natural logarithm, so \(H\) is measured in nats
2. Common Failure Modes and Solutions
| Symptom |
Likely Cause |
Solution |
| \(\operatorname{ESS} \ll N_s\) |
Poor proposal choice |
Use heavier-tailed proposal or AIS |
| High variance across runs |
Weight degeneracy |
Increase samples or improve proposal |
| Biased estimates |
Violated support condition |
Ensure \(q(\mathbf{x}) \gt 0\) wherever \(\pi(\mathbf{x}) \gt 0\) |
| Numerical instability |
Extreme weight ratios |
Work in log-space and use the log-sum-exp trick |
3. Validation Techniques
To verify an importance sampling implementation:
- Known ground truth. Test on problems with analytical solutions
- Convergence check. Estimate should stabilize as \(N_s\) increases
- Multiple proposals. Different valid proposals should give consistent estimates
- Bootstrap confidence intervals. Resample the pairs \((\mathbf{x}_n, \tilde{w}_n)\) to assess uncertainty
Importance sampling, SNIS, and AIS complete our treatment of Monte Carlo methods for Bayesian computation. Throughout the Monte Carlo
and importance sampling material, we have worked with parametric models where uncertainty is over a finite-dimensional parameter
\(\boldsymbol{\theta}\). In Gaussian Processes, we take a fundamentally
different perspective. Instead of placing priors over parameters, we place priors directly over functions, and obtain
nonparametric Bayesian models with built-in uncertainty quantification.