> ## Documentation Index
> Fetch the complete documentation index at: https://aegean.ai/llms.txt
> Use this file to discover all available pages before exploring further.

# Importance Sampling

> Expectations from weighted samples of a proposal distribution, and resampling.

[![Open In Colab](https://colab.research.google.com/assets/colab-badge.svg)](https://colab.research.google.com/github/pantelis/eng-ai-agents/blob/main/notebooks/rse/importance-sampling/index.ipynb)

*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.

```python theme={null}
%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$.

```python theme={null}
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}")
```

```output theme={null}
importance sampling estimate of the mean: 5.630
true mean of the target:                  5.500
unweighted mean of the samples:           4.934
```

![Output from cell 3](https://mintcdn.com/aegeanaiinc/izJXkwFvXGX6Sn6q/aiml-common/lectures/rse/importance-sampling/images/cell_3_output_1.svg?fit=max\&auto=format\&n=izJXkwFvXGX6Sn6q\&q=85\&s=86974815ed5b97b1962d3a1b50d403a5)

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.

```python theme={null}
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}")
```

```output theme={null}
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
```

```output theme={null}
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.

```python theme={null}
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}")
```

```output theme={null}
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.

```python theme={null}
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}")
```

```output theme={null}
distinct samples kept after resampling: 2296 of 5000
mean of the resampled set: 5.491
```

![Output from cell 7](https://mintcdn.com/aegeanaiinc/izJXkwFvXGX6Sn6q/aiml-common/lectures/rse/importance-sampling/images/cell_7_output_1.svg?fit=max\&auto=format\&n=izJXkwFvXGX6Sn6q\&q=85\&s=8ac11e8961b8094b5b0a27863acad519)

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](https://aegean.ai/aiml-common/lectures/reinforcement-learning/policy-based-algorithms/PPO/ppo-explained-minimal/ppo_explained_minimal.md) 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](https://github.com/aegean-ai/eaia-feedback/issues/new?template=page-issue.yml\&page=/aiml-common/lectures/rse/importance-sampling).


This documentation is built and hosted on [Mintlify](https://mintlify.com), a developer documentation platform.