Importance Sampling

Importance sampling estimates an expectation under one distribution from samples of another. Each sample is weighted by the ratio of the two densities. Resampling in proportion to the weights gives an unweighted sample set that follows the first distribution. The example in this section is one-dimensional.

%config InlineBackend.figure_formats = ["svg"]

import matplotlib.pyplot as plt
import numpy as np
from cycler import cycler

plt.rcParams.update(
    {
        "figure.figsize": (6.5, 3.5),  # matches the 6.5 in text column
        "font.size": 9,
        "axes.titlesize": 10,
        "axes.labelsize": 9,
        "xtick.labelsize": 8,
        "ytick.labelsize": 8,
        "legend.fontsize": 8,
        "axes.prop_cycle": cycler(
            color=[
                "#0072B2",
                "#D55E00",
                "#009E73",
                "#cb77a6",
                "#c28600",
                "#1d99e0",
                "#7f7f7f",
            ]
        ),
        "axes.grid": True,
        "grid.alpha": 0.3,
        "svg.fonttype": "path",
    }
)

# named colors from the palette above
BLUE, VERMILLION, GREEN, PINK, GOLD, SKY, GRAY = (
    "#0072B2",
    "#D55E00",
    "#009E73",
    "#cb77a6",
    "#c28600",
    "#1d99e0",
    "#7f7f7f",
)

np.random.seed(13)

Target and proposal

Suppose you want the expectation of a function

h

under a density

f

\mathbb{E}_f[h(x)] = \int h(x) \, f(x) \, dx.

With samples

x^{[1]}, \ldots, x^{[M]}

drawn from

f

, the average of

h(x^{[m]})

estimates this integral. Often there is no direct way to sample from

f

, for example when

f

is the posterior belief of a Bayes filter. Importance sampling draws the samples from another density

g

, the proposal distribution.

f

is the target distribution.

g

must be positive wherever

f

is positive, so that every value

f

can produce can also be drawn from

g

. Multiplying and dividing by

g

gives

\mathbb{E}_f[h(x)] = \int h(x) \, \frac{f(x)}{g(x)} \, g(x) \, dx = \mathbb{E}_g\left[h(x) \, w(x)\right], \qquad w(x) = \frac{f(x)}{g(x)}.

The ratio

w(x)

is the importance weight. With samples

x^{[m]} \sim g

, the estimate is

\mathbb{E}_f[h(x)] \approx \frac{1}{M} \sum_{m=1}^{M} w(x^{[m]}) \, h(x^{[m]}).

Where

f > g

, the proposal undersamples and the weight is above 1. Where

f < g

, the weight is below 1. The example target has two modes, like a localization belief when two places match the robot’s observations:

f(x) = 0.3 \, \mathcal{N}(x; 2, 0.5^2) + 0.7 \, \mathcal{N}(x; 7, 1^2).

The proposal is a single wide Gaussian,

g(x) = \mathcal{N}(x; 5, 3^2)

. The mean of the target is

0.3 \cdot 2 + 0.7 \cdot 7 = 5.5

def normal_pdf(x, mean, std):
    return np.exp(-0.5 * ((x - mean) / std) ** 2) / (std * np.sqrt(2 * np.pi))


def target_pdf(x):
    """Two-mode target density f."""
    return 0.3 * normal_pdf(x, 2.0, 0.5) + 0.7 * normal_pdf(x, 7.0, 1.0)


target_mean = 0.3 * 2.0 + 0.7 * 7.0

g_mean, g_std = 5.0, 3.0
M = 1000
x = np.random.normal(g_mean, g_std, M)  # samples from the proposal g
w = target_pdf(x) / normal_pdf(x, g_mean, g_std)  # importance weights f/g

estimate = np.mean(w * x)
print(f"importance sampling estimate of the mean: {estimate:.3f}")
print(f"true mean of the target:                  {target_mean:.3f}")
print(f"unweighted mean of the samples:           {x.mean():.3f}")

importance sampling estimate of the mean: 5.630
true mean of the target:                  5.500
unweighted mean of the samples:           4.934

The samples in the top plot cover the range of the proposal, and their unweighted average estimates the mean of

g

. The weights in the bottom plot are large under the two modes of

f

and near zero elsewhere. The weighted average is 5.63, against the true mean of 5.5.

Self-normalized weights

In filtering, the target is known only up to a constant. Bayes’ rule gives the posterior as

\eta \, p(z \mid x) \, p(x)

, and the normalizer

\eta

is an integral over the state space that the filter does not compute. Write the target as

f(x) = \eta \tilde{f}(x)

, where

\tilde{f}

can be evaluated and

\eta

is unknown. Compute the weights from

\tilde{f}

and divide by their sum:

\tilde{w}^{[m]} = \frac{\tilde{f}(x^{[m]})}{g(x^{[m]})}, \qquad \mathbb{E}_f[h(x)] \approx \frac{\sum_{m=1}^{M} \tilde{w}^{[m]} \, h(x^{[m]})}{\sum_{m=1}^{M} \tilde{w}^{[m]}}.

\eta

cancels between numerator and denominator. This self-normalized estimate is slightly biased for finite

M

, and the bias vanishes as

M

grows. Its error typically shrinks as

1/\sqrt{M}

, with a constant that depends on how similar

g

is to

f

(Thrun, Burgard and Fox, 2005, Section 4.2.2). The code below scales the target by an arbitrary factor, as if

\eta

were unknown, and repeats the estimate 500 times for each sample size.

def unnormalized_target(x):
    """The target scaled by a constant the estimator does not know."""
    return 17.0 * target_pdf(x)


def self_normalized_mean(M, g_mean, g_std):
    x = np.random.normal(g_mean, g_std, M)
    w = unnormalized_target(x) / normal_pdf(x, g_mean, g_std)
    w = w / w.sum()
    return np.sum(w * x)


for M in [10, 100, 1000, 10000]:
    estimates = np.array([self_normalized_mean(M, g_mean, g_std) for _ in range(500)])
    rmse = np.sqrt(np.mean((estimates - target_mean) ** 2))
    print(f"M = {M:>5d}   mean of estimates = {estimates.mean():.3f}   RMSE = {rmse:.3f}")

M =    10   mean of estimates = 5.510   RMSE = 1.085
M =   100   mean of estimates = 5.512   RMSE = 0.336
M =  1000   mean of estimates = 5.507   RMSE = 0.111

M = 10000   mean of estimates = 5.500   RMSE = 0.035

Each tenfold increase in

M

reduces the error by about

\sqrt{10} \approx 3.2

. The factor 17 does not change the estimate.

Choosing the proposal

When

g

is close to

f

, the weights are nearly equal and every sample contributes. When

g

puts little mass where

f

is large, the few samples that land there get very large weights, and the estimate depends mostly on them. The effective sample size, computed with normalized weights,

\hat{N}_{\mathrm{eff}} = \frac{1}{\sum_{m=1}^{M} \left(w^{[m]}\right)^2},

measures this. It equals

M

when all weights are equal and 1 when a single sample holds all the weight. The code compares three proposals with

M = 1000

samples each: the wide Gaussian used so far, a narrower Gaussian centered between the modes, and a Gaussian that covers only the right mode.

proposals = {
    "wide, N(5, 3^2)": (5.0, 3.0),
    "narrow, N(5, 1.5^2)": (5.0, 1.5),
    "right mode only, N(7, 1^2)": (7.0, 1.0),
}

M = 1000
for name, (mean, std) in proposals.items():
    estimates, ess = [], []
    for _ in range(500):
        x = np.random.normal(mean, std, M)
        w = unnormalized_target(x) / normal_pdf(x, mean, std)
        w = w / w.sum()
        estimates.append(np.sum(w * x))
        ess.append(1.0 / np.sum(w**2))
    rmse = np.sqrt(np.mean((np.array(estimates) - target_mean) ** 2))
    print(f"{name:<28s} mean ESS = {np.mean(ess):6.1f}   mean estimate = {np.mean(estimates):.3f}   RMSE = {rmse:.3f}")

wide, N(5, 3^2)              mean ESS =  501.2   mean estimate = 5.502   RMSE = 0.109
narrow, N(5, 1.5^2)          mean ESS =  288.5   mean estimate = 5.505   RMSE = 0.176
right mode only, N(7, 1^2)   mean ESS =  906.1   mean estimate = 6.867   RMSE = 1.496

The wide proposal has an effective sample size of about 500 of 1000, and its estimate is close to 5.5. The narrow proposal has little mass near

x = 2

, where the left mode is, so the few samples there get large weights. Its effective sample size is about 290 and its error is larger. The proposal centered on the right mode almost never produces a sample near

x = 2

, so the left mode is missing from the weighted set. The weights of its samples are nearly equal, the effective sample size is about 900, and the estimate is 6.9 instead of 5.5. The effective sample size measures how even the weights of the drawn samples are. It cannot reveal a region of the target that the proposal never reaches.

Resampling

Some algorithms need an unweighted sample set, for example to pass the samples through another sampling step. Resampling produces one from the weighted set. Draw

M

samples with replacement from the weighted set, each with probability equal to its normalized weight. Samples with large weights are copied several times and samples with small weights tend to disappear. The new set is approximately distributed according to

f

, and every sample in it has the same weight.

M = 5000
x = np.random.normal(g_mean, g_std, M)
w = unnormalized_target(x) / normal_pdf(x, g_mean, g_std)
w = w / w.sum()

idx = np.random.choice(M, size=M, replace=True, p=w)
x_resampled = x[idx]

print(f"distinct samples kept after resampling: {len(np.unique(idx))} of {M}")
print(f"mean of the resampled set: {x_resampled.mean():.3f}")

distinct samples kept after resampling: 2296 of 5000
mean of the resampled set: 5.491

The histogram of the resampled set follows the target. Fewer than half of the original samples survive, many of them as several copies. The mean of the resampled set is no more accurate than the weighted mean before it, and the random draw adds variance. Resampling moves the samples into the regions where the target has mass, where a later step can use them.

Importance ratios in policy optimization

PPO collects actions with an old policy

\pi_{\text{old}}

and updates a new policy

\pi_\theta

. The ratio

\pi_\theta(a_t \mid s_t) / \pi_{\text{old}}(a_t \mid s_t)

in its objective is the importance weight that corrects for the samples coming from the old policy, and the clipping limits how large that weight can get.

References

Thrun, S., Burgard, W., and Fox, D. (2005). Probabilistic Robotics. MIT Press.

Found an error on this page? Report an issue.

State Estimation & Mapping

State Estimation

SLAM

Importance Sampling

Target and proposal

Self-normalized weights

Choosing the proposal

Resampling

Importance ratios in policy optimization

References

​Target and proposal

​Self-normalized weights

​Choosing the proposal

​Resampling

​Importance ratios in policy optimization

​References

Target and proposal

Self-normalized weights

Choosing the proposal

Resampling

Importance ratios in policy optimization

References