particles.smc_samplers¶
Classical and wastefree 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 nonsequential Bayesian inference, rareevent simulation, etc. For more background on (standard) SMC samplers, see Chapter 17 (and references therein). For the wastefree variant, see Dau & Chopin (2022).
The following type of sequences of distributions are implemented:
SMC tempering: target distribution at time t has 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 statespace model. For each thetaparticle, a local particle filter is run to approximate the likelihood up to time t; see Chapter 18 in the book.
Nested sampling: target at time t is prior constrained to likelihood being above threshold l_t, with an increasing sequence of l_t’s (inspired by Salomone et al, 2018).
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 subclassing 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 loglikelihood 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, bidimensional 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 wastefree 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). See Dau & Chopin (2022) for more details (or the notebook on SMC samplers).
if you use the wastefree version, you may also compute a singlerun estimate of the variance of a particular estimate; see
var_wf
,Var_logLt
andVar_phi
for more details.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 array of shape (N,d))X.lpost
: a numpy float array of length N, which stores the logtarget density of the N particles.X.shared
: a dictionary that contains metainformation 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 hyperparameters of) the MCMC kernel to the weighted sample (W, x).self.step(x)
: apply a single step to theThetaParticles
objectx
, 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 wastefree 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 wastefree SMC, see
subclass MCMCSequenceWF
) or only the states of the last iteration (as in
standard SMC, see subclass AdaptiveMCMCSequence
).
The bottom line is: if you wish to implement a different MCMC scheme to move
the particles, you should subclass ArrayMCMC
. If you wish to implement a new
strategy to repeat several MCMC steps, you should subcass MCMCSequence (or one
of its subclasses).
References
Dau, H.D. and Chopin, N. Wastefree Sequential Monte Carlo, J. R. Stat. Soc. Ser. B. Stat. Methodol. 84, 1 (2022), 114–148, arxiv:2011.02328, doi:10.1111/rssb.12475
Functions

Returns the list [l[i] for i in idx] When needed, objects l[i] are replaced by a copy, to make sure that the elements of the list are all distinct. 



Find next annealing exponent by solving ESS(exp(lw)) = alpha * N. 

Turns record array arr into a dict 

Singlerun variance estimate for wastefree SMC. 

Returns a view to record array theta which behaves like a (N,d) float array. 
Classes

MCMC sequence for a standard SMC sampler (keep only final states). 

FeynmanKac class for adaptive tempering SMC. 

Independent Metropolis (Gaussian proposal). 

Base class for a (single) MCMC step applied to an array. 

Base class for Metropolis steps (whatever the proposal). 

Gaussian random walk Metropolis. 

Base FeynmanKac class for SMC samplers. 

A list that implements fancy indexing, and forces elements to be distinct. 



Importance sampler. 

Base class for a (fixed length or adaptive) sequence of MCMC steps. 

MCMC sequence for a wastefree SMC sampler (keep all intermediate states). 

FeynmanKac subclass for the SMC^2 algorithm. 

Base class for static models. 

FeynmanKac class for tempering SMC (fixed exponents). 



Base class for particle systems for SMC samplers. 

Collects singlerun estimates of the variance of log L_t (normalising cst) from a wastefree sampler (Dau & Chopin, 2022). 

Collects singlerun estimates of the variance of a given estimate (for a certain function phi) from a wastefree sampler (Dau & Chopin, 2022). 