Sequential Monte Carlo

A hidden Markov model

A 1-D random walk seen through noise: a latent state \(X_t\) takes a Gaussian step each tick, and each step emits one noisy reading \(y_t\).

model = """
(let [x1 (sample (normal 0 1))   o1 (observe (normal x1 0.5) -0.7638)
      x2 (sample (normal x1 1))  o2 (observe (normal x2 0.5)  0.0611)
      x3 (sample (normal x2 1))  o3 (observe (normal x3 0.5) -0.4634)
      ...                        ...] x30)              # T = 30 readings, one per step
"""

This is a hidden Markov model: the latent states form a Markov chain and each reading depends only on the current state, so the joint factorizes as

\[p(X_{1:T}, y_{1:T}) = p(X_1)\,\prod_{t=2}^{T} p(X_t \given X_{t-1})\;\prod_{t=1}^{T} p(y_t \given X_t).\]

Single-site MH samples it by editing one \(X_t\) of a finished trace and re-running:

chain, acc = single_site_mh(model, rng, S=20000)   # run-to-completion evaluator (companion notebook)
# E[x_30] = -12.34   (exact: -12.60)   effective sample size = 33 of 20000

\(20{,}000\) edits \(\approx 33\) independent draws. The factors arrive in sequence: what would it take to use that?

One hard problem, a sequence of easy ones

Truncate that factorization at step \(t\): the first \(t\) states and readings define an intermediate target, the joint of everything seen so far.

\[\gamma_1\subseteq\gamma_2\subseteq\dots\subseteq\gamma_T=p(Y,X),\qquad \gamma_t(X_{1:t})=p(y_{1:t},\,X_{1:t}).\]

• \(\gamma_1\) is low-dimensional and easy to hit; \(\gamma_T\) is the full posterior we want. The cuts interpolate between them.
• So one hard problem becomes a sequence of easy ones: advance many runs together and correct each after every reading, while there is still time to act, instead of nudging one finished trace forever.

The bootstrap particle filter

Carry \(L\) partial execution particles. Resume them together until the program ends; at each observe:

1. Advance: each draws its new state \(X_t\) from the transition prior \(p(X_t \given X_{t-1})\) and stops at the reading.
2. Weight by the one new reading, \(W_t^{(\ell)} = p(y_t \given X_t^{(\ell)})\) (the prior draw of \(X_t\) cancels).
3. Resample draws \(L\) ancestors with probability \(\propto W_t\). Heavy particles duplicate, near-zero ones drop.
4. Resume, and stop once every particle reaches the end.

def smc(program, L, seed=0):
    g = np.random.default_rng(seed)
    ms = [M(program, np.random.default_rng(int(s))) for s in g.integers(1, 2**62, L)]
    while True:
        outs = [resume(m) for m in ms]                    # advance each to its next observe
        if outs[0][0] == "done":                          # programs finished: return final states
            return np.array([o[1] for o in outs])
        inc = np.array([o[1] for o in outs])              # incremental weights: log p(y_t | X_t)
        W = np.exp(inc - inc.max()); W /= W.sum()         # normalize
        idx = g.choice(L, size=L, p=W)                    # resample: ancestors proportional to weight
        ms = [ms[a].fork(np.random.default_rng(int(s)))   # duplicate survivors, fresh rng each
              for a, s in zip(idx, g.integers(1, 2**62, L))]

The driver uses two operations the language must provide: resume (advance a particle to its next observe, then pause) and m.fork (duplicate a paused particle).

Question 1

In step 2 a particle has just drawn its new state \(X_t\) from the transition prior and reached the \(t\)-th observe.

What is particle \(\ell\)’s incremental weight \(W_t^{(\ell)}\) there?

Vote: https://pe.app/ppls

Question 2

After a few steps one particle holds almost all the weight and the rest are near zero. A classmate says: “just duplicate the heavy particle \(L\) times to refill the set, then keep going.”

Why does duplication alone not fix the degeneracy?

Vote: https://pe.app/ppls

Question 3

Make execution resumable as an explicit stack machine: a control stack C (operations left to do) and a value stack V (results). resume pops C until an observe should pause it; by then the distribution and the observed value are on V.

def resume(m):                                 # advance until the next observe pauses the machine
    C, V = m.C, m.V
    while C:
        t, *rest = C.pop()
        ...
        elif t == "observek":                  # V = [..., dist, value]
            ...                                # <- your code: weight by this reading, then pause
    return ("done", V[-1])                     # control stack empty: the program finished

In your notebook, write the observek case so the machine pauses at the reading and can be resumed later. Hint: pop value and dist, add dist.log_prob(value) to m.log_w, push the value back, then return so the leftover C and V are the rest of the computation.

Next class

Answer key