Importance Sampling and SMC

This notebook demonstrates the importance sampling (IS) and Sequential Monte Carlo (SMC) tools in psiphy.mcmc.importance_sampling, applied to the built-in toy models.

What is importance sampling?

We want to draw samples from a target \(p(\theta | x)\) but can only evaluate it up to a constant. IS draws samples from a tractable proposal \(q(\theta)\) and corrects with weights:

\[w_i = \frac{p(\theta_i | x)}{q(\theta_i)}, \qquad \theta_i \sim q\]

The effective sample size (ESS) measures how much information the weighted set carries:

\[\text{ESS} = \frac{1}{\sum_i w_i^2} \in [1, N]\]

Low ESS means the proposal is a poor match for the posterior — most weight falls on very few particles.

SMC fixes this by annealing from the prior to the posterior through a sequence of tempered targets, resampling and applying MCMC moves when ESS drops too low.

Output format

importance_sampling, sequential_importance_sampling and sequential_monte_carlo (SMC) all return a dict:

Key	Meaning
samples, weights	weighted particles and their normalised weights (aligned)
equal_weight_samples	particles resampled by weight — use for plain means, stds, predictive draws
ess	effective sample size
log_evidence	log marginal likelihood (SMC; IS when log_proposal_fn is given), else None

Always pass samples together with weights to weighted statistics and plots.

import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import norm

from psiphy.mcmc import (
    UniformPrior,
    eff_sample_size,
    importance_sampling,
    sequential_importance_sampling,
    sequential_monte_carlo,
)
from psiphy.toy_models import GaussianSignal, NoisyLine

np.random.seed(42)

---

1 UniformPrior helper

UniformPrior wraps a box prior and provides sample(n) and log_prob(theta) — the two methods expected by importance_sampling and sequential_monte_carlo.

prior = UniformPrior({'mean': (-5, 5), 'sigma': (0.1, 5)})
print('ndim  :', prior.ndim)
print('sample:', prior.sample(3))
print('log_prob in support  :', prior.log_prob([0.0, 1.0]))
print('log_prob out of support:', prior.log_prob([10.0, 1.0]))

ndim  : 2
sample: [[-1.25459881  4.7585001 ]
 [ 2.31993942  3.03342657]
 [-3.4398136   0.86437315]]
log_prob in support  : -3.891820298110627
log_prob out of support: -inf

---

2 Importance Sampling — GaussianSignal

Target: posterior \(p(\mu, \sigma \mid x)\) for a Gaussian signal with \(N=40\) observations.

gs = GaussianSignal(N=40)
theta_true = np.array([1.5, 0.8])
x_obs = gs.simulate(theta_true)

print(f'x_obs  mean={x_obs.mean():.3f}  std={x_obs.std():.3f}  (true: {theta_true})')

x_obs  mean=1.250  std=0.734  (true: [1.5 0.8])

# Broad uniform prior — deliberate mismatch to show ESS degradation
prior_broad = UniformPrior([(-5, 5), (0.1, 5)])
log_posterior = lambda t: gs.log_prob(t, x_obs) + prior_broad.log_prob(t)

proposal_broad = prior_broad.sample(5000)
res = importance_sampling(log_posterior, proposal_broad)
s_broad, w_broad = res['samples'], res['weights']
ess_broad = eff_sample_size(w_broad)
print(f'Broad proposal  — ESS = {ess_broad:.1f} / 5000 ({100*ess_broad/5000:.1f}%)')

Broad proposal  — ESS = 11.4 / 5000 (0.2%)

# Tighter prior centred near the posterior — higher ESS
prior_tight = UniformPrior([(0.0, 3.0), (0.3, 2.0)])
log_posterior_tight = lambda t: gs.log_prob(t, x_obs) + prior_tight.log_prob(t)

proposal_tight = prior_tight.sample(5000)
res = importance_sampling(log_posterior_tight, proposal_tight)
s_tight, w_tight = res['samples'], res['weights']
ess_tight = eff_sample_size(w_tight)
print(f'Tight proposal  — ESS = {ess_tight:.1f} / 5000 ({100*ess_tight/5000:.1f}%)')

Tight proposal  — ESS = 114.1 / 5000 (2.3%)

# Analytic posterior on mean for comparison (sigma plugged in as sample std)
post_mean, post_sigma = gs.analytic_posterior(x_obs)

fig, axes = plt.subplots(1, 2, figsize=(11, 4))

mu_grid = np.linspace(0.5, 2.5, 300)

for ax, samples, w, label, color in [
    (axes[0], proposal_broad, w_broad, 'broad prior', 'C0'),
    (axes[1], proposal_tight, w_tight, 'tight prior', 'C1'),
]:
    ess = eff_sample_size(w)
    ax.hist(samples[:, 0], weights=w, bins=40, density=True,
            alpha=0.6, color=color, label=f'IS ({label}, ESS={ess:.0f})')
    ax.plot(mu_grid, norm.pdf(mu_grid, post_mean, post_sigma),
            'k--', lw=1.5, label='analytic posterior')
    ax.axvline(theta_true[0], c='red', lw=1.5, ls=':', label='true mean')
    ax.set(xlabel=r'$\mu$', title=f'IS — {label}')
    ax.legend(fontsize=8)

plt.suptitle('Importance Sampling: posterior on mean', y=1.01)
plt.tight_layout()
plt.show()

Key takeaway: the broad uniform proposal samples mostly from regions of low likelihood, yielding very few effective samples. A tighter proposal that overlaps the posterior better gives a far higher ESS and more accurate estimate.

Why is the posterior not centred on the true :math:`mu`? The posterior describes this data set. The 40 observed values have a sample mean that, by chance, lies about two standard errors below the true mean, and every correct method — including the analytic posterior — follows the data there. A single data set cannot reveal the truth more precisely than its own noise allows; whether a method is biased is a question about many data sets (Section 7).

se = x_obs.std(ddof=1) / np.sqrt(len(x_obs))
print(f'sample mean = {x_obs.mean():.3f}, true mean = {theta_true[0]}, '
      f'standard error = {se:.3f}  ->  data offset = {(x_obs.mean() - theta_true[0]) / se:+.1f} sigma')

sample mean = 1.250, true mean = 1.5, standard error = 0.118  ->  data offset = -2.1 sigma

---

3 Effective Sample Size as a diagnostic

# Show how ESS changes with proposal width
half_widths = np.linspace(0.3, 5.0, 20)
ess_list = []

for hw in half_widths:
    p = UniformPrior([
        (theta_true[0] - hw, theta_true[0] + hw),
        (max(0.1, theta_true[1] - hw), theta_true[1] + hw),
    ])
    prop = p.sample(2000)
    res = importance_sampling(lambda t: gs.log_prob(t, x_obs) + p.log_prob(t), prop)
    _, w = res['samples'], res['weights']
    ess_list.append(eff_sample_size(w) / 2000 * 100)

fig, ax = plt.subplots(figsize=(7, 3.5))
ax.plot(half_widths, ess_list, 'o-', c='steelblue')
ax.set(xlabel='Proposal half-width', ylabel='ESS / N  (%)',
       title='ESS vs proposal width (N=2000 particles)')
ax.axhline(100, ls='--', c='k', lw=0.8, label='perfect (ESS=N)')
ax.legend()
plt.tight_layout()
plt.show()

---

4 Importance Sampling — NoisyLine

Two parameters: slope and intercept.

nl = NoisyLine(Nx=40, noise_sigma=0.5)
theta_nl = np.array([-0.8, 5.0])
y_obs = nl.simulate(theta_nl)

prior_nl = UniformPrior([(-3, 3), (0, 10)])
log_post_nl = lambda t: nl.log_prob(t, y_obs) + prior_nl.log_prob(t)

prop_nl = prior_nl.sample(8000)
res = importance_sampling(log_post_nl, prop_nl)
s_nl, w_nl = res['samples'], res['weights']
ess_nl = eff_sample_size(w_nl)
print(f'NoisyLine IS  ESS={ess_nl:.1f}/8000  '
      f'slope={np.average(s_nl[:,0], weights=w_nl):.3f} (true {theta_nl[0]})  '
      f'intercept={np.average(s_nl[:,1], weights=w_nl):.3f} (true {theta_nl[1]})')

NoisyLine IS  ESS=4.2/8000  slope=-0.832 (true -0.8)  intercept=5.280 (true 5.0)

fig, axes = plt.subplots(1, 2, figsize=(10, 4))

for ax, idx, name, true_val in [
    (axes[0], 0, 'slope', theta_nl[0]),
    (axes[1], 1, 'intercept', theta_nl[1]),
]:
    ax.hist(s_nl[:, idx], weights=w_nl, bins=40, density=True, alpha=0.6, label='IS posterior')
    ax.axvline(true_val, c='red', lw=1.5, ls='--', label='true value')
    ax.set(xlabel=name, title=f'NoisyLine IS — {name}')
    ax.legend()

plt.suptitle(f'NoisyLine: IS posterior marginals  (ESS={ess_nl:.0f}/8000)', y=1.01)
plt.tight_layout()
plt.show()

---

5 SMC — annealing from prior to posterior

SMC overcomes the low-ESS problem by bridging from the prior to the posterior through tempered targets \(p_t \propto \text{prior} \times \mathcal{L}^{\beta_t}\). When ESS drops below a threshold, particles are resampled and a Gaussian random-walk MCMC move is applied to maintain diversity.

# GaussianSignal with broad prior — compare IS vs SMC
prior_broad2 = UniformPrior([(-5, 5), (0.1, 5)])

# IS baseline
prop = prior_broad2.sample(1000)
res = importance_sampling(
    lambda t: gs.log_prob(t, x_obs) + prior_broad2.log_prob(t), prop
)
s_is, w_is = res['samples'], res['weights']
print(f'IS  ESS={eff_sample_size(w_is):.1f}/1000')

# SMC
res = sequential_monte_carlo(
    lambda t: gs.log_prob(t, x_obs),
    prior_broad2,
    n_particles=1000,
    n_steps=10,
    verbose=False,
)
s_smc, w_smc = res['samples'], res['weights']
print(f'SMC ESS={eff_sample_size(w_smc):.1f}/1000')

IS  ESS=1.5/1000
SMC ESS=762.0/1000

fig, axes = plt.subplots(1, 2, figsize=(11, 4))

for ax, name, idx, true_val in [
    (axes[0], r'$\mu$', 0, theta_true[0]),
    (axes[1], r'$\sigma$', 1, theta_true[1]),
]:
    ax.hist(s_is[:, idx], weights=w_is, bins=30, density=True,
            alpha=0.5, label=f'IS (ESS={eff_sample_size(w_is):.0f})', color='C0')
    ax.hist(s_smc[:, idx], weights=w_smc, bins=30, density=True,
            alpha=0.5, label=f'SMC (ESS={eff_sample_size(w_smc):.0f})', color='C2')
    ax.axvline(true_val, c='red', lw=1.5, ls='--', label='true value')
    ax.set(xlabel=name, title=f'IS vs SMC — {name}')
    ax.legend(fontsize=8)

plt.suptitle('Broad prior: IS (poor ESS) vs SMC (annealed, high ESS)', y=1.01)
plt.tight_layout()
plt.show()

---

6 2-D posterior scatter: SMC on NoisyLine

prior_nl2 = UniformPrior([(-3, 3), (0, 10)])

res = sequential_monte_carlo(
    lambda t: nl.log_prob(t, y_obs),
    prior_nl2,
    n_particles=1000,
    n_steps=10,
    verbose=False,
)

smc_nl, w_nl2 = res['samples'], res['weights']   # weighted particles
smc_eq = res['equal_weight_samples']              # resampled: use for plain mean/std
print(f'NoisyLine SMC  ESS={res["ess"]:.1f}/1000')
print(f'slope     = {smc_eq[:,0].mean():.3f} ± {smc_eq[:,0].std():.3f}  (true {theta_nl[0]})')
print(f'intercept = {smc_eq[:,1].mean():.3f} ± {smc_eq[:,1].std():.3f}  (true {theta_nl[1]})')

NoisyLine SMC  ESS=839.4/1000
slope     = -0.825 ± 0.024  (true -0.8)
intercept = 5.225 ± 0.147  (true 5.0)

fig, axes = plt.subplots(1, 2, figsize=(11, 4))

sc = axes[0].scatter(smc_nl[:, 0], smc_nl[:, 1],
                     c=w_nl2, cmap='viridis', s=8, alpha=0.7)
axes[0].scatter(*theta_nl, c='red', s=100, zorder=5, label='true theta')
plt.colorbar(sc, ax=axes[0], label='weight')
axes[0].set(xlabel='slope', ylabel='intercept',
            title='SMC posterior — NoisyLine')
axes[0].legend()

axes[1].plot(nl.x, y_obs, '.', alpha=0.6, label='data')
for s, i in smc_eq[np.random.choice(len(smc_eq), 80)]:
    axes[1].plot(nl.x, s * nl.x + i, alpha=0.05, c='C1', lw=0.8)
axes[1].plot(nl.x, theta_nl[0] * nl.x + theta_nl[1], 'r-', lw=2, label='true line')
axes[1].set(xlabel='x', ylabel='y', title='Posterior predictive draws')
axes[1].legend()

plt.tight_layout()
plt.show()

---

7 Is the posterior biased? Coverage over many data sets

A single posterior can sit 1–2σ from the truth just because of noise in that data set. A method is unbiased and calibrated if, repeated over many independent data sets,

• the offset of the posterior mean from the truth, in units of posterior std (\(z = (\bar\theta_{\rm post} - \theta_{\rm true}) / \sigma_{\rm post}\)), averages to zero with a spread of about 1, and

• the 68% (95%) credible interval contains the truth in about 68% (95%) of the data sets.

sequential_monte_carlo with a vectorised log-likelihood (one call evaluates all particles) makes this cheap: 300 independent data sets take about a second.

n_datasets = 300
prior_cov = UniformPrior([(0.0, 3.0), (0.1, 2.5)])
rng = np.random.default_rng(1)

z, in68, in95 = [], [], []
for _ in range(n_datasets):
    x = rng.normal(theta_true[0], theta_true[1], size=gs.N)          # a fresh data set

    def loglike(th, x=x):                                             # th: (n_particles, 2)
        mu, sig = th[:, :1], th[:, 1]
        return -len(x) * np.log(sig) - 0.5 * np.sum((x[None, :] - mu) ** 2, axis=1) / sig**2

    res = sequential_monte_carlo(loglike, prior_cov, n_particles=1000, n_steps=10,
                                 verbose=False, vectorized=True)
    s = res['equal_weight_samples']
    z.append((s.mean(0) - theta_true) / s.std(0))
    lo68, hi68, lo95, hi95 = np.percentile(s, [16, 84, 2.5, 97.5], axis=0)
    in68.append((lo68 <= theta_true) & (theta_true <= hi68))
    in95.append((lo95 <= theta_true) & (theta_true <= hi95))
z, in68, in95 = map(np.array, (z, in68, in95))

print(f"{'':10s} {'mean z':>8s} {'std z':>7s} {'68% cover':>10s} {'95% cover':>10s}")
for i, name in enumerate(['mu', 'sigma']):
    print(f"{name:10s} {z[:, i].mean():+8.2f} {z[:, i].std():7.2f} {in68[:, i].mean():10.0%} {in95[:, i].mean():10.0%}")
print(f"(mean z uncertainty: ±{1 / np.sqrt(n_datasets):.2f})")

             mean z   std z  68% cover  95% cover
mu            -0.08    0.95        70%        96%
sigma         +0.10    0.97        68%        95%
(mean z uncertainty: ±0.06)

from scipy.stats import norm as _norm

fig, axes = plt.subplots(1, 2, figsize=(11, 3.5))
grid = np.linspace(-4, 4, 200)
for i, (ax, name) in enumerate(zip(axes, [r'$\mu$', r'$\sigma$'])):
    ax.hist(z[:, i], bins=30, range=(-4, 4), density=True, alpha=0.6, label=f'{n_datasets} data sets')
    ax.plot(grid, _norm.pdf(grid), 'k--', label='N(0, 1)')
    ax.set(xlabel=f'z for {name}  (posterior mean − truth) / posterior std')
axes[0].axvline((x_obs.mean() - theta_true[0]) / se, c='r', label='this tutorial\'s data set')
axes[0].legend(fontsize=8)
plt.tight_layout()
plt.show()

The offsets scatter around zero with unit spread and the credible intervals have the stated coverage: the method is calibrated. The data set used in the sections above (red line) is simply in the tail — about 1 in 20 data sets is at least that far off. The small positive mean \(z\) for \(\sigma\) is expected: with a flat prior the posterior of a scale parameter is right-skewed, so its mean sits slightly above the peak.

This repeated-data test is the idea behind simulation-based calibration (SBC), which is the standard check for likelihood-free and neural inference methods.

---

Summary

Method	When to use	ESS behaviour
importance_sampling	Proposal already close to posterior	Degrades quickly if proposal is too broad
sequential_importance_sampling	Moderate mismatch, iterative refinement	Better than plain IS; no annealing schedule
sequential_monte_carlo	Broad prior, multi-modal, or expensive likelihood	Maintains high ESS via annealing + MCMC move

All three accept any log_target_fn / log_likelihood_fn and a UniformPrior (or any object with .sample(n) and .log_prob(theta)).