Three algorithms, three evaluators

Likelihood weighting, Metropolis-Hastings, and sequential Monte Carlo were each written as their own evaluator. Compare what they do at each expression form:

Only sample and observe differ. Everything else is evaluated the same way, every time.

The message interface

resume(m) runs the machine until the next probabilistic effect, then returns a message:

return ('sample',  address, distribution, m)         # reached (sample d)
return ('observe', address, distribution, value, m)  # reached (observe d v)
return ('done',    value, m)                          # program finished

The controller replies by pushing a value, and resumes:

def send(m, value):
    m.V.append(value)
resume(m)

The control stack C and value stack V inside m are everything left to do: the machine is the continuation. Pausing is returning from resume; resuming is calling it again; forking is copying the stacks.

Same control state, two spellings

The book writes the continuation in continuation-passing style, as a function k. Our evaluator holds the same control state as explicit stacks.

The new idea is not how to evaluate let or if; it is where inference lives: outside the evaluator, as a controller over messages.

Question 1

LW, MH, and SMC were three separate evaluators; now one runtime serves all three. What makes that possible?

Question 2

Single-site MH and SMC can both be implemented by re-running the program from the start at each step, conditioning on stored values. What does pausing and forking a live execution mainly buy over that?

Why pause and fork, not re-run

A controller could instead re-run the program from the start at each step, conditioning on stored values. That is what made earlier implementations expensive:

single-site MH : re-run the whole program to change one site
SMC            : re-run from the top at every observe  ->  quadratic in observations

Pausing keeps the partial execution; forking copies it and continues both halves. Re-execution of the deterministic prefix disappears, and SMC becomes linear.

Activity 1: implement the SSMH controller

Single-site MH is a third controller over the same runtime. It keeps an address-keyed trace, resamples one site per step, and reuses the rest by address. The dimension-corrected ratio is given as mh_log_alpha(...): the point is not to rederive it, but to see that the controller has enough information to compute it from the two address-keyed traces.

def run(program, rng, x0, cache):
    m = initial_machine(program, rng); X, S, O = {}, {}, {}
    while True:
        msg = resume(m); tag = msg[0]
        if tag == 'sample':
            _, a, d, m = msg
            x = ...                          # resample the selected/new site, else reuse cache[a]
        elif tag == 'observe':
            _, a, d, y, m = msg
            O[a] = log_prob(d, y); m.V.append(y)
        elif tag == 'done':
            _, value, m = msg; return value, X, S, O

In your notebook, write the sample reply policy and the accept/reject loop.

x = d.sample(rng) if (a == x0 or a not in cache) else cache[a]         # reuse by address

value, X, S, O = run(program, rng, None, {})                           # initial trace
for _ in range(steps):
    a0 = list(X)[rng.integers(len(X))]                                 # pick one site to change
    v2, X2, S2, O2 = run(program, rng, a0, X)                          # propose; reuse the rest
    if np.log(rng.random()) < mh_log_alpha(X, X2, S, S2, O, O2, a0):
        value, X, S, O = v2, X2, S2, O2                                # accept

A model the first-order language cannot express

(defn geom []
  (if (sample (bernoulli 0.3))
      0
      (+ 1 (geom))))
(geom)

Each call draws one coin; a failure recurses and draws again. The number of sample calls is decided by the run itself, and it differs from run to run.

So the model has no fixed list of random variables: it denotes an unbounded number of them. The program can no longer be compiled to a graph, because the graph would need infinitely many nodes, and eager evaluation of the if would never terminate.

Why the message runtime still works

Graph compilation fails, but evaluation does not. The reason is simple:

Likelihood weighting needs only complete runs, so it transfers unchanged: draw at each sample reached, score at each observe, read the return value. Each terminating run generates its own finite trace. The SMC controller uses the same interface, but it assumes particles reach a compatible sequence of observe breakpoints.

A function value is a closure

Recursion needs defn to refer to itself; first-class functions go further, letting a function be returned and used elsewhere:

(let [make-shift (fn [mu] (fn [x] (+ x mu)))
      f          (make-shift 10)]
  (f 3))            ; => 13

The returned f still needs mu = 10 after make-shift has returned.

Addressing a moving target

A controller matches choices across runs by address. With dynamic control flow, a flat counter breaks: a different branch renumbers every site after it.

The runtime instead names each effect by the path the machine took to reach it: which binding, which branch, which argument, which call. That structural address keeps a site’s name stable across runs even when the trace changes shape, which is exactly what single-site MH needs to reuse and propose.

It also exposes a hazard: if a proposal makes the run take a different branch, a site that the current trace addressed may not be reached at all in the proposed run.

Question 4

The geom program can call sample an unbounded number of times, so the model denotes infinitely many possible random variables and cannot be compiled to a finite graph. Yet likelihood weighting runs on it unchanged. Why?

Question 5

(let [make-shift (fn [mu] (fn [x] (+ x mu)))
      f          (make-shift 10)]
  (f 3))          ; => 13, though mu is not in scope at the call site

(f 3) returns 13 even though mu is not bound where f is called. This works because:

Question 6

When applying a closure, which environment should be extended with the argument bindings?

Activity 2: implement closure application

The closure-application case of the evaluator is incomplete. new_env starts from the closure’s captured environment; the step that binds the formal parameters to the actual arguments is missing.

elif t == 'callk':
    _, n, addr = instr
    args = [V.pop() for _ in range(n)][::-1]; f = V.pop()
    if isinstance(f, Closure):
        ...                                 # <- your code
    else:
        V.append(f(*args))

In your notebook, write the binding step.

for p, arg in zip(f.params, args): new_env[p] = arg

Next class

Summary

Answer key