smc_samplers

Classical and waste-free SMC samplers.

Overview

This module implements SMC samplers, that is, SMC algorithms that may sample from a sequence of arbitrary probability distributions (and approximate their normalising constants). Applications include sequential and non-sequential Bayesian inference, rare-event simulation, etc. For more background on (standard) SMC samplers, see Chapter 17 (and references therein). For the waste-free variant, see Dau & Chopin (2020).

The following type of sequences of distributions are implemented:

  • SMC tempering: target distribution at time t as a density of the form mu(theta) L(theta)^{gamma_t}, with gamma_t increasing from 0 to 1.
  • IBIS: target distribution at time t is the posterior of parameter theta given data Y_{0:t}, for a given model.
  • SMC^2: same as IBIS, but a for state-space model. For each theta-particle, a local particle filter is run to approximate the likelihood up to time t; see Chapter 18 in the book.

SMC samplers for binary distributions (and variable selection) are implemented elsewhere, in module binary_smc.

Before reading the documentation below, you might want to have a look at the following notebook tutorial, which may be a more friendly introduction.

Target distribution(s)

If you want to use a SMC sampler to perform Bayesian inference, you may specify your model by sub-classing StaticModel, and defining method logpyt (the log density of data Y_t, given previous datapoints and parameter values) as follows:

class ToyModel(StaticModel):
    def logpyt(self, theta, t):  # density of Y_t given parameter theta
        return -0.5 * (theta['mu'] - self.data[t])**2 / theta['sigma2']

In this example, theta is a structured array, with fields named after the different parameters of the model. For the sake of consistency, the prior should be a distributions.StructDist object (see module distributions for more details), whose inputs and outputs are structured arrays with the same fields:

from particles import distributions as dists

prior = dists.StructDist(mu=dists.Normal(scale=10.),
                         sigma2=dists.Gamma())

Then you can instantiate the class as follows:

data = np.random.randn(20)  # simulated data
my_toy_model = ToyModel(prior=prior, data=data)

This object may be passed as an argument to the FeynmanKac classes that define SMC samplers, see below.

Under the hood, class StaticModel defines methods loglik and logpost which computes respectively the log-likelihood and the log posterior density of the model at a certain time.

What if I don’t want to do Bayesian inference

This is work in progress, but if you just want to sample from some target distribution, using SMC tempering, you may define your target as follows:

class ToyBridge(TemperingBridge):
    def logtarget(self, theta):
        return -0.5 * np.sum(theta**2, axis=1)

and then define:

base_dist = dists.MvNormal(scale=10., cov=np.eye(10))
toy_bridge = ToyBridge(base_dist=base_dist)

Note that, this time, we went for standard, bi-dimensional numpy arrays for argument theta. This is fine because we use a prior object that also uses standard numpy arrays.

FeynmanKac objects

SMC samplers are represented as FeynmanKac classes. For instance, to perform SMC tempering with respect to the bridge defined in the previous section, you may do:

fk_tpr = AdaptiveTempering(model=toy_bridge, len_chain=100)
alg = SMC(fk=fk_tpr, N=200)
alg.run()

This piece of code will run a tempering SMC algorithm such that:

  • the successive exponents are chosen adaptively, so that the ESS between two successive steps is cN, with c=1/2 (use parameter ESSrmin to change the value of c).
  • the waste-free version is implemented; that is, the actual number of particles is 100 * 200, but only 200 particles are resampled at each time, and then moved through 100 MCMC steps (parameter len_chain) (set parameter wastefree=False to run a standard SMC sampler).
  • the default MCMC strategy is random walk Metropolis, with a covariance proposal set to a fraction of the empirical covariance of the current particle sample. See next section for how to use a different MCMC kernel.

To run IBIS instead, you may do:

fk_ibis = IBIS(model=toy_model, len_chain=100)
alg = SMC(fk=fk_ibis, N=200)

Again see the notebook tutorials for more details and examples.

Under the hood

ThetaParticles

In a SMC sampler, a particle sample is represented as a ThetaParticles object X, which contains several attributes such as, e.g.:

  • X.theta: a structured array of length N, representing the N particles (or alternatively a numpy arrray of shape (N,d))
  • X.lpost: a numpy float array of length N, which stores the log-target density of the N particles.
  • X.shared: a dictionary that contains meta-information on the N particles; for instance it may be used to record the successive acceptance rates of the Metropolis steps.

Details may vary in a given algorithm; the common idea is that attribute shared is the only one which not behave like an array of length N. The main point of ThetaParticles is to implement fancy indexing, which is convenient for e.g. resampling:

from particles import resampling as rs

A = rs.resampling('multinomial', W)  # an array of N ints
Xp = X[A]  # fancy indexing

MCMC schemes

A MCMC scheme (e.g. random walk Metropolis) is represented as an ArrayMCMC object, which has two methods:

  • self.calibrate(W, x): calibrate (tune the hyper-parameters of) the MCMC kernel to the weighted sample (W, x).
  • self.step(x): apply a single step to the ThetaParticles object x, in place.

Furthermore, the different ways one may repeat a given MCMC kernel is represented by a MCMCSequence object, which you may pass as an argument when instantiating the FeynmanKac object that represents the algorithm you want to run:

move = MCMCSequenceWF(mcmc=ArrayRandomWalk(), len_chain=100)
fk_tpr = AdaptiveTempering(model=toy_bridge, len_chain=100, move=move)
# run a waste-free SMC sampler where the particles are moved through 99
# iterations of a random walk Metropolis kernel

Such objects may either keep all intermediate states (as in waste-free SMC, see sub-class MCMCSequenceWF) or only the states of the last iteration (as in standard SMC, see sub-class AdaptiveMCMCSequence).

The bottom line is: if you wish to implement a different MCMC scheme to move the particles, you should sub-class ArrayMCMC. If you wish to implement a new strategy to repeat several MCMC steps, you should sub-cass MCMCSequence (or one of its sub-classes).

References

Dau, H.D. and Chopin, N (2020). Waste-free Sequential Monte Carlo, arxiv:2011.02328

Module summary

StaticModel([data, prior]) Base class for static models.
ThetaParticles([shared]) Base class for particle systems for SMC samplers.
IBIS([model, wastefree, len_chain, move])
TemperingBridge([base_dist])
AdaptiveTempering([model, wastefree, …]) Feynman-Kac class for adaptive tempering SMC.
SMC2([ssm_cls, prior, data, smc_options, …]) Feynman-Kac subclass for the SMC^2 algorithm.
ArrayMCMC() Base class for a (single) MCMC step applied to an array.
MCMCSequence([mcmc, len_chain]) Base class for a (fixed length or adaptive) sequence of MCMC steps.
MCMCSequenceWF([mcmc, len_chain]) MCMC sequence for a waste-free SMC sampler (keep all intermediate states).