smoothing

This module implements:

  1. particle history classes, which store the full or partial history of a SMC algorithm.
  2. off-line smoothing algorithms as methods of these classes.

For on-line smoothing, see instead the collectors module.

History classes

A SMC object has a hist attribute, which is used to record at certain times t:

  • the N current particles \(X_t^n\);
  • their weights;
  • (optionally, see below), the ancestor variables \(A_t^n\).

The frequency at which history is recorded depends on option store_history of class SMC. Possible options are:

  • True: records full history (at every time t);
  • False: no history (attribute hist set to None);
  • callable f: history is recorded at time t if f(t) returns True
  • int k: records a rolling window history of length k (may be used to perform fixed-lag smoothing)

This module implements different classes that correspond to the different cases:

All these classes provide a similar interface. If smc is a SMC object, then:

  • smc.hist.X[t] returns the N particles at time t
  • smc.hist.wgts[t] returns the N weights at time t (see resampling.weights)
  • smc.hist.A[t] returns the N ancestor variables at time t

Partial History

Here are some examples on one may record history only at certain times:

# store every other 10 iterations
smc = SMC(fk=fk, N=100, store_history=lambda t: (t % 10) == 0)

# store at certain times given by a list
times = [10, 30, 84]
smc = SMC(fk=fk, N=100, store_history=lambda t: t in times)

Once the algorithm is run, smc.hist.X and smc.hist.wgts are dictionaries, the keys of which are the times where history was recorded. The ancestor variables are not recorded in that case:

smc.run()
smc.hist.X[10]  # the N particles at time 10
smc.hist.A[10]  # raises an error

Full history, off-line smoothing algorithms

For a given state-space model, off-line smoothing amounts to approximate the distribution of the complete trajectory \(X_{0:T}\), given data \(y_{0:T}\), at some fixed time horizon T. The corresponding algorithms take as an input the complete history of a particle filter, run until time T (forward pass). Say:

# forward pass
fk = ssm.Bootstrap(ssm=my_ssm, data=y)
pf = particles.SMC(fk=fk, N=100, store_history=True)
pf.run()

Then, pf.hist is an instance of class particles.smoothing.ParticleHistory, which has the following methods:

  • backward_sampling: implements O(N) and O(N^2) FFBS algorithms, which generates smoothing trajectories from the history of the forward pass;
  • backward_sampling_qmc: same as above, but for when the forward pass was based on QMC (quasi-Monte Carlo).
  • two_filter_smoothing: to estimate expectations of marginal smoothing distributions (using the two-filter smoothing approach).

For more details, see the documentation of particles.smoothing.ParticleHistory, the ipython notebook on smoothing, and Chapter 12 of the book.

Warning

the complete history of a particle filter may take a lot of memory.

Rolling history, Fixed-lag smoothing

To obtain a rolling window (fixed-length) history:

smc = SMC(fk=fk, N=100, store_history=10)
smc.run()

In that case, fields smc.hist.X, smc.hist.wgts and smc.hist.A are deques of max length 10. Using negative indices:

smc.hist.X[-1]  # the particles at final time T
smc.hist.X[-2]  # the particles at time T - 1
# ...
smc.hist.X[-10] # the N particles at time T - 9
smc.hist.X[-11] # raises an error

Note that this type of history makes it possible to perform fixed-lag smoothing as follows:

B = smc.hist.compute_trajectories()
# B[t, n] is index of ancestor of X_T^n at time t
phi = lambda x: x  # any test function
est = np.average(phi(smc.hist.X[-10][B[-10, :]]), weights=smc.W)
# est is an estimate of E[ phi(X_{T-9}) | Y_{0:T}]

Note

recall that it is possible to run SMC algorithms step by step, since they are iterators. Hence it is possible to do fixed-lag smoothing step-by-step as well.

Module summary

ParticleHistory(fk, qmc) Particle history.
PartialParticleHistory(func) Partial history.
RollingParticleHistory(length) Rolling window history.
smoothing_worker([method, N, fk, fk_info, …]) Generic worker for off-line smoothing algorithms.