particles.collectors¶
Objects that collect summaries at each iteration of a SMC algorithm.
Overview¶
This module implements “summary collectors”, that is, objects that collect at every time t certain summaries of the particle system. Important applications are fixed-lag smoothing and on-line smoothing. However, the idea is a bit more general that that. Here is a simple example:
import particles
from particles import collectors as col
# ...
# define some_fk_model
# ...
alg = particles.SMC(fk=some_fk_model, N=100,
collect=[col.Moments(), col.Online_smooth_naive()])
alg.run()
print(alg.summaries.moments) # list of moments
print(alg.summaries.naive_online_smooth) # list of smoothing estimates
Once the algorithm is run, the object alg.summaries
contains the computed
summaries, stored in lists of length T (one component for each iteration t).
Note that:
argument
collect
expects a list of Collector objects;the name of the collector classes are capitalised, e.g.
Moments
;by default, the name of the corresponding summaries are not, e.g.
pf.summaries.moments
.
Default summaries¶
By default, the following summaries are collected (even if argument collect
is not used):
ESSs
: ESS (effective sample size) at each iteration;
rs_flags
: whether resampling was triggered or not at each time t;
logLts
: log-likelihood estimates.
For instance:
print(alg.summaries.ESSs) # sequence of ESSs
You may turn off summary collection entirely:
alg = particles.SMC(fk=some_fk_model, N=100, collect='off')
This might be useful in very specific cases when you need to keep a large number of SMC objects in memory (as in SMC^2). In that case, even the default summaries might take too much space.
Computing moments¶
To compute moments (functions of the current particle sample):
def f(W, X): # expected signature for the moment function
return np.average(X, weights=W) # for instance
alg = particles.SMC(fk=some_fk_model, N=100,
collect=[Moments(mom_func=f)])
Without an argument, i.e. Moments()
, the collector computes the default
moments defined by the FeynmanKac
object; for instance, for a FeynmanKac
object derived from a state-space model, the default moments at time t consist
of a dictionary, with keys 'mean'
and 'var'
, containing the particle
estimates (at time t) of the filtering mean and variance.
It is possible to define different defaults for the moments. To do so,
override method default_moments
of the considered FeynmanKac class:
from particles import state_space_models as ssms
class Bootstrap_with_better_moments(ssms.Bootstrap):
def default_moments(W, X):
return np.average(X**2, weights=W)
# ...
# define state-space model my_ssm
# ...
my_fk_model = Bootstrap_with_better_moments(ssm=my_ssm, data=data)
alg = particles.SMC(fk=my_fk_model, N=100, moments=True)
In that case, my_fk_model.summaries.moments
is a list of weighed averages
of the squares of the components of the particles.
Fixed-lag smoothing¶
Fixed-lag smoothing means smoothing of the previous h states; that is, computing (at every time t) expectations of
for a fixed integer h (at times t <= h; if t<h, replace h by t).
This requires keeping track of the h previous states for each particle;
this is achieved by using a rolling window history, by setting option
store_history
to an int equals to h+1 (the length of the trajectories):
alg = particles.SMC(fk=some_fk_model, N=100,
collect=[col.Fixed_lag_smooth(phi=phi)],
store_history=3) # h = 2
See module smoothing
for more details on rolling window and other types of
particle history. Function phi takes as an input the N particles, and returns
a numpy.array:
def phi(X):
return np.exp(X - 2.)
If no argument is provided, test function \(\varphi(x)=x\) is used.
Note however that X is a deque of length at most h; it behaves like a list,
except that its length is always at most h + 1. Of course this function
could simply return its arguments W
and X
; in that case you simply
record the fixed-lag trajectories (and their weights) at every time t.
On-line smoothing¶
On-line smoothing is the task of approximating, at every time t, expectations of the form:
On-line smoothing is covered in Sections 11.1 and 11.3 in the book. Note that on-line smoothing is typically restricted to additive functions \(\phi\), see below.
The following collectors implement online-smoothing algorithms:
Online_smooth_naive
: basic forward smoothing (carry forward full trajectories); cost is O(N) but performance may be poor for large t.Online_smooth_ON2
: O(N^2) on-line smoothing. Expensive (cost is O(N^2), so big increase of CPU time), but better performance.Paris
: on-line smoothing using Paris algorithm. (Warning: current implementation is very slow, work in progress).
These algorithms compute the smoothing expectation of a certain additive function, that is a function of the form:
The elementary function \(\psi_t\) is specified by defining method
add_func
in considered state-space model. Here is an example:
class BootstrapWithAddFunc(ssms.Bootstrap):
def add_func(self, t, xp, x): # xp means x_{t-1} (p=past)
if t == 0:
return x**2
else:
return (xp - x)**2
The reason why additive functions are specified in this way is that
additive functions often depend on fixed parameters of the state-space model
(which are available in the closure of the StateSpaceModel
object, but
not outside).
The two first algorithms do not have any parameter, the third one (Paris) have one (default: 2). To use them simultaneously:
alg = particles.SMC(fk=some_fk_model, N=100,
collect=[col.Online_smooth_naive(),
col.Online_smooth_ON2(),
col.Paris(Nparis=5)])
Variance estimators¶
The variance estimators of Chan & Lai (2013), Lee & Whiteley (2018), etc., are
implemented as collectors in module variance_estimators
; see the
documentation of that module for more details.
User-defined collectors¶
You may implement your own collectors as follows:
import collectors
class Toy(collectors.Collector):
# optional, default: toy (same name without capital)
summary_name = 'toy'
# signature of the __init__ function (optional, default: {})
signature = {phi=None}
# fetch the quantity to collect at time t
def fetch(self, smc): # smc is the particles.SMC instance
return np.mean(self.phi(smc.X))
Once this is done, you may use this new collector exactly as the other ones:
alg = particles.SMC(N=30, fk=some_fk_model, collect=[col.Moments(), Toy(phi=f)])
Then pf.summaries.toy
will be a list of the summaries collected at each
time by the fetch
method.
Classes
|
Base class for collectors. |
|
|
|
Compute some function of fixed-lag trajectories. |
|
|
|
Collects empirical moments (e.g. |
|
Mix-in for on-line smoothing algorithms. |
|
|
|
|
|
Hybrid version of the Paris algorithm. |
|
|
|
Class to store and update summaries. |