The Particle Filter

This page is for the reader who wants the machinery. The implementation is clocks.inference.ParticleFilter: adaptive tempered resample-move sequential Monte Carlo (SMC) for a static parameter vector.

State and prior support

\(N\) particles represent complete hypotheses \(\theta_i\) with normalized weights \(w_i\). The public point-mass API uses a uniform rectangular prior conditioned on every actual support rule:

positions and masses lie inside the configured PriorConfig ranges;
multiple masses obey the canonical strict ordering by their first spatial coordinate, \(x_{1,1} < \cdots < x_{K,1}\); and
every candidate satisfies the project’s \(|2\Phi| \le 0.1\) weak-field policy at every clock.

Initialization rejection-samples from that normalized conditional prior. Proposals outside it have log density \(-\infty\); they are rejected rather than clipped, reflected, sorted, or repaired. A directly constructed ParticleFilter likewise requires a prior_sampler and log_prior_density that describe the same distribution. If they do not, posterior and evidence results are undefined.

One observation, several tempered stages

Suppose observations \(y_1,\ldots,y_{t-1}\) have been completed and \(y_t\) has just arrived. One public update() advances through the intermediate targets

\[ \pi_{t,\beta}(\theta) \propto p(\theta)\prod_{s<t}p(y_s\mid\theta) p(y_t\mid\theta)^\beta, \qquad 0 \le \beta \le 1. \]

At a stage with exponent \(\beta\), the filter chooses the largest next exponent \(\beta' \le 1\) whose reweighted effective sample size remains at the configured target ess_target * N; bisection finds \(\beta'\) when the full remaining likelihood would be too sharp. With \(\ell_i = \log p(y_t\mid\theta_i)\),

\[ \widetilde w_i = w_i\exp[(\beta'-\beta)\ell_i], \qquad \mathrm{ESS} = \frac{1}{\sum_i (\widetilde w_i/\sum_j\widetilde w_j)^2}. \]

All normalization is done in log space. Completed independent Gaussian observations are retained as sufficient statistics, so an MH target evaluation does not replay the entire observation history.

Resampling and invariant moves

Before any resampling, the filter records the stage’s normalizing-constant increment

\[ \log \widehat{Z}_{\beta:\beta'} = \log\sum_i w_i\exp[(\beta'-\beta)\ell_i]. \]

An intermediate stage is then resampled to equal weights. A final stage is also resampled if its ESS reaches the configured threshold. Systematic (the default), stratified, and residual resampling are available.

Every resampling is followed by rejuvenation_steps random-walk Metropolis-Hastings steps. For a stage, the code freezes one regularized empirical covariance, shared by all particles, and scales it by proposal_scale**2 / D. Equivalently, its Cholesky factor (the proposal’s standard-deviation scale) is multiplied by proposal_scale / sqrt(D). The resulting Gaussian proposal is symmetric, so for a proposal \(\theta'\) the acceptance probability is exactly

\[ \alpha(\theta,\theta') = \min\left\{1, \frac{\pi_{t,\beta'}(\theta')}{\pi_{t,\beta'}(\theta)} \right\}. \]

This accept/reject move leaves the current intermediate target invariant. Merely perturbing resampled copies without the MH correction would not.

Evidence

Summing the pre-resampling log increments across stages and observations gives the SMC estimate of the log normalizing constant. This is a consistent Monte Carlo evidence estimator as particle count grows; at any finite \(N\) it has sampling error. It is neither deterministic nor an exact finite-particle marginal likelihood. Model comparison combines each model’s estimate with its normalized model prior.

The filter exposes log_evidence, log_evidence_history, and the most recent last_log_evidence_increments. Per-update diagnostics report the number of tempering stages plus MH proposal and acceptance counts; public inference results retain those diagnostics in each history entry.

Watch tempering and ESS

Code

import matplotlib.pyplot as plt
import numpy as np
from clocks import (
    ClockArray, InferenceConfig, MassConfig, NoiseConfig,
    PriorConfig, SimulationConfig, build_particle_filter, simulate,
)

clock_array = ClockArray(positions=np.array([[-5.0], [0.0], [5.0]]), track_offset=1.0)
truth = MassConfig(positions=np.array([[2.5]]), masses=np.array([0.10]))
sim = simulate(SimulationConfig(
    clock_array=clock_array, ground_truth=truth,
    noise=NoiseConfig(observation_std=0.005), n_observations=25, seed=42,
))
pf = build_particle_filter(InferenceConfig(
    clock_array=clock_array, noise=NoiseConfig(observation_std=0.005),
    prior=PriorConfig(position_range=(-8.0, 8.0), mass_range=(0.005, 0.15)),
    n_particles=400, n_masses=1,
    ess_target=0.8, rejuvenation_steps=2, proposal_scale=2.38, seed=42,
))

ess_trace = []
stage_trace = []
for obs in sim.observations:
    state = pf.update(obs)
    ess_trace.append(1.0 / np.sum(state.weights**2))
    stage_trace.append(pf.last_diagnostics.tempering_stages)

steps = np.arange(1, len(ess_trace) + 1)
fig, (ax_ess, ax_stages) = plt.subplots(2, 1, sharex=True)
ax_ess.plot(steps, ess_trace, marker="o", ms=3)
ax_ess.axhline(
    pf.ess_target * pf.n_particles,
    ls="--", color="gray", label="ESS target",
)
ax_ess.set_ylabel("final ESS")
ax_ess.legend()
ax_stages.step(steps, stage_trace, where="mid")
ax_stages.set_xlabel("observation #")
ax_stages.set_ylabel("tempering stages")
plt.close(fig)
fig

Final ESS and internal tempering-stage count for each public observation update.

What the diagnostics can and cannot show

Weight collapse raises RuntimeError instead of propagating a NaN. Low MH acceptance or many tempering stages can diagnose poor proposal geometry or a sharp update, but no single diagnostic proves convergence. The corrected SMC controls and acceptance tolerances are frozen from the development seeds, and the reserved 400–411 block has now been run exactly once. That fixed-seed protocol is regression and calibration evidence for those particular cases, not a population reliability estimate.

With multiple masses, ordering removes label switching from the represented parameterization. Proposals crossing the strict ordering boundary are rejected. Near that boundary, sampling can still be difficult.