Course:CPSC522/Sequential Monte Carlo samplers

Sequential Monte Carlo samplers

This page first discusses the Sequential Monte Carlo technique for performing "online" Bayesian inference - inference performed continuously as observations are made. We then describe Sequential Monte Carlo samplers, which build on this technique to develop an algorithm for efficient "offline" inference.

Principal Author: William Harvey
Collaborators:

Abstract

Sequential Monte Carlo allows for fast online inference by representing the posterior distribution at each step with a set of samples. Since each distribution in the sequence is typically closely related to the previous distribution, it is possible to use the set of samples from the posterior at time $t$ as a starting point for estimating the posterior at time $t+1$. This is the key to the Sequential Monte Carlo approach. A more recent development concerns the use of this algorithm for offline inference: when we wish to approximately sample from a single distribution. By defining this distribution as simply the last distribution in a sequence, it is possible to leverage Sequential Monte Carlo as a powerful and general-purpose inference algorithm.

Builds on

Probability is fundamental as background for Sequential Monte Carlo (SMC). Building on probability, Markov chain Monte Carlo methods are used as an important building block of the SMC algorithm. Also, note that SMC is a generalisation of particle filtering (which considers only state-space models) to any dynamic system.

Related Pages

As a powerful general method for Bayesian inference, SMC can be used in models such as HMMs.

Sequential Monte Carlo for Online Inference

Sequential Monte Carlo (SMC) has been developed as a general method for inference in models where information is obtained sequentially. It allows a posterior distribution over unobserved variables to be estimated as the information is obtained, and efficiently updated for each new observation. As with a particle filter, the core idea is to represent the posterior with a set of approximate samples from it. This idea allows the representation of arbitrarily complex and multi-modal posteriors over the states of systems with arbitrarily complex dynamics. It is defined for systems which can be specified as follows:

Model Definition

The model can be described as a sequence of probability distributions $\pi _{t}(\mathbf {x} _{t})$ for $t=0,1,2,\ldots$ Each of these refers to the posterior distribution over the latent variables, $\mathbf {x} _{t}$, given the observations made up until time $t$. Note that the dimension of $\mathbf {x} _{t}$ can change with $t$ (as in the example below, where it increases with $t$ to allow it to store all states $1,\ldots ,t$).

State-space example

To give a concrete example of such a model fitting the above definition, we now describe a state-space model. This has a $d$-dimensional latent state, $x_{t}$, which evolves over time with known stochastic dynamics, $p(x_{t+1}|x_{t})$. Observations of the system, $y_{t}$, are made according to $p(y_{t}|x_{t})$, and used to condition the estimate of $x$. $x$ is initially distributed according to $p(x_{0})$. We define $\mathbf {x} _{t}$ to be the set of all $x_{i}$ for $i=0,1,\ldots ,t$, so that $\mathbf {x_{t}}$ essentially contains the "history" of the particle. Note that this is a special case and, in SMC in general, $\mathbf {x_{t}}$does not "remember" all previous states. Similarly, we denote $\mathbf {y} _{t}$ to be the set of all $y_{i}$ for $i=0,1,\ldots ,t$. Note that this means that the dimension of $\mathbf {x} _{t}$ increases by $d$ with each increase in $t$. This gives the following expression for the joint distribution of $\mathbf {x} _{t}$ and $y_{1:t}$:

$p(\mathbf {x} _{t},\mathbf {y} _{t})={\frac {1}{Z}}p(x_{0})\prod _{i=1}^{t}p(y_{i}|x_{i})p(x_{i}|x_{i-1})$

Since the posterior $\pi _{t}(\mathbf {x} _{t})$ is proportional to the joint probability (down to a normalising constant independent of $\mathbf {x} _{t}$), we can write it as:

$\pi _{t}(\mathbf {x} _{t})={\frac {1}{Z}}p(x_{0})\prod _{i=1}^{t}p(y_{i}|x_{i})p(x_{i}|x_{i-1})$

In the case that each distribution is Gaussian, with variance independent of $x_{t}$, and the means of $p(y_{t}|x_{t})$ and $p(x_{t+1}|x_{t})$ are affine transformations of $x_{t}$, it is possible to calculate this posterior analytically with a Kalman filter. When these conditions are not met, SMC can be used to instead compute set of samples, with which expectations over the posterior can be approximated. Note that, in this example, SMC is equivalent to using a particle filter.

Inference

SMC involves performing inference by maintaining a set of $m$ samples, $\mathbf {x} _{t}^{(j)}$ for $j=0,\ldots ,m$, which approximate expectations over $\pi _{t}(\mathbf {x} _{t})$. At the beginning of inference, these are initialised, e.g. by sampling from a prior before making any observations. Thereafter, they are updated at each time-step to approximate $\pi _{t+1}(\mathbf {x} _{t+1})$ rather than $\pi _{t}(\mathbf {x} _{t})$. In this section, we discuss two different procedures for updating the samples, although many more exist. A rough outline of the procedure is as follows:

1. Intialise samples (e.g. by sampling from a prior over the initial state)
2. for each timestep:
   1. transition approximate samples of $\pi _{t-1}(\mathbf {x} _{t-1})$ to samples of $\pi _{t}(\mathbf {x} _{t})$
   2. optional: resample $\pi _{t}(\mathbf {x} _{t})$ to increase effective sample size

Sequential Importance Sampling

For the discussion of sequential importance sampling (SIS), we will consider the case where $\mathbf {x} _{t}$ grows in dimension with $t$, and so $\mathbf {x} _{t+1}$ can be written as $\{x_{t+1},\mathbf {x_{t}} \}$. Liu et al. advocate sampling each $x_{t+1}^{(j)}$ from $p(x_{t+1}|\mathbf {y} _{t+1}\mathbf {x} _{t}^{(j)})$ for each $j$, where $\mathbf {x} _{t}^{(j)}$ is the $j$th sample from $\pi _{t}(\mathbf {x} _{t})$. This yields samples $\mathbf {x} _{t+1}^{(j)}=\{x_{t+1}^{(j)},\mathbf {x_{t}} ^{(j)}\}$ from the distribution:

$\pi _{t}(\mathbf {x} _{t})p(x_{t+1}|\mathbf {y} _{t+1},\mathbf {x} _{t})$

However, we want samples which approximate the distribution:

$\pi _{t+1}(\mathbf {x} _{t+1})=\pi _{t+1}(\mathbf {x} _{t})p(x_{t+1}|\mathbf {y} _{t+1},\mathbf {x} _{t})$

Intuitively, our samples do not condition any of $\mathbf {x} _{t}^{(j)}$ on the new observation, $y_{t+1}$. To correct for this, the samples are weighted. By weighting samples, we can approximate expectations over a distribution, $p(x)$, with weighted samples from $q(x)$. The approximation is accurate in the sense that, for properly chosen weights:

$\mathbb {E} _{p(x)}\left[f(x)\right]=\mathbb {E} _{q(x)}\left[w(x)f(x)\right]$

It can be shown that the weights which give this property are given by $w(x)={\frac {p(x)}{q(x)}}$, so can be used whenever the probability density functions of the true and approximating function are known. If expectations are normalised by the expected value of the weights, the weights can be calculated using unnormalised densities:

$\mathbb {E} _{p(x)}\left[f(x)\right]={\frac {\mathbb {E} _{q(x)}\left[w(x)f(x)\right]}{\mathbb {E} \left[w(x)\right]}}$

Since we want to approximate $\pi _{t+1}(\mathbf {x} _{t+1})$ with samples from $\pi _{t}(\mathbf {x} _{t})p(x_{t+1}|\mathbf {y} _{t+1},\mathbf {x} _{t})$, we must weight them by:

$w(\mathbf {x} _{t+1})={\frac {\pi _{t+1}(\mathbf {x} _{t+1})}{\pi _{t}(\mathbf {x} _{t})p(x_{t+1}|\mathbf {y} _{t+1},\mathbf {x} _{t})}}={\frac {\pi _{t+1}(\mathbf {x} _{t})p(x_{t+1}|\mathbf {y} _{t+1},\mathbf {x} _{t})}{\pi _{t}(\mathbf {x} _{t})p(x_{t+1}|\mathbf {y} _{t+1},\mathbf {x} _{t})}}={\frac {\pi _{t+1}(\mathbf {x} _{t})}{\pi _{t}(\mathbf {x} _{t})}}$

We then have a set of weighted samples, $\{\mathbf {x} _{t+1}^{(j)},w_{t+1}^{(j)}\}$ for $j=0,\ldots ,m$, which approximate $\pi _{t+1}(\mathbf {x} _{t+1})$. These weights can be used when it is possible to sample exactly from $p(x_{t+1}|\mathbf {y} _{t+1}\mathbf {x} _{t}^{(j)})$. If we are not able to, we can instead sample from an approximation, $q(x_{t+1}|\mathbf {y} _{t+1}\mathbf {x} _{t}^{(j)})$, and calculate weights so that the approximation of $\pi _{t+1}(\mathbf {x} _{t+1})$ is still exact (in the limit as $m\rightarrow \infty$). This leads to weights of the form:

$w(\mathbf {x} _{t+1})={\frac {\pi _{t+1}(\mathbf {x} _{t})p(x_{t+1}|\mathbf {y} _{t+1},\mathbf {x} _{t})}{\pi _{t}(\mathbf {x} _{t})q(x_{t+1}|\mathbf {y} _{t+1},\mathbf {x} _{t})}}$

These weights may be kept and multiplied by later weights to approximate the posterior at later timesteps. Alternatively, the posterior can be resampled (as discussed later on), allowing the weights to be discarded.

Rejection Sampling

Rejection sampling is based on a similar idea to SIS but, instead of storing weights, samples are accepted with a probability proportional to the weight. Samples are continually drawn until $m$ have been accepted. The procedure is as follows:

1. Sample an index, $j$, with probability proportional to $w_{t}^{(j)}$ (or from a uniform distribution if the samples are not weighted)
2. Draw $x_{t+1}^{(j)}$ from $g(x_{t+1}|\mathbf {x} _{t}^{(j)},\mathbf {y} _{t+1})$, which should approximate $\pi _{t+1}(x_{t+1}|\mathbf {x} _{t}^{(j)})$ as well as possible for efficiency
3. Accept ($j$, $x_{t+1}^{(j)}$) with probability $p={\frac {\pi _{t+1}(\mathbf {x} _{t}^{(j)},x_{t+1})}{c_{t+1}\pi _{t}(\mathbf {x} _{t}^{(j)})g(x_{t+1}|\mathbf {x} _{t},\mathbf {y} _{t+1})}}$ where $c_{t+1}$ is a constant chosen to ensure that $p\leq 1$.

This procedure is performed repeatedly until $m$ samples have been accepted. Variations of this method exist which can improve efficiency (e.g. the proportion of samples of accepted) with techniques such as Rao-Blackwellisation.

Resampling

One issue with the sequential importance sampling approach (and other, related approaches), is that the variance of the weights can grow large when they are continually changed over many time steps. Therefore, although we may have $m$ samples of $\pi _{t}(\mathbf {x} _{t})$, many of them may have a very small weight, and so our approximation could be dominated by very few samples (or even just one). One method to avoid this problem is to resample them from an empirical distribution given by:

${\hat {\pi }}_{t}(\mathbf {x} _{t})=\sum _{j=1}^{m}{\frac {w_{t}^{(j)}}{\sum _{j'=1}^{m}w_{t}^{(j')}}}\delta _{\mathbf {x} _{t}^{(j)}}(\mathbf {x} _{t})$

where $\delta$ is the Kronecker delta. This is simple to sample from by sampling from a categorical distribution over $j$ with probabilities proportional to the weights, and maintains the property that $\mathbf {x} _{t}^{(j)}$ are exact samples from $\pi _{t}(\mathbf {x} _{t})$in the limit as $m\rightarrow \infty$. Since this distribution, can be sampled from exactly, the resulting samples do not have weights, and so the resampling step gets rid of the weights, increasing the effective sample size. This procedure effectively gets rid of samples with small weights, and duplicates samples with high weights. The duplicated samples will diverge at the next timestep, so this procedure allows the most "likely" regions of sample space to be explored more. However, the resampling procedure also introduces noise to the procedure, and so it is typically not recommended to perform it at every step. Liu et al. suggest either resampling at fixed intervals (with a size specific to the problem), or checking the effective sample size (a measure of the degeneracy of the weights) after each transition, and resampling when this decreases below a threshold. Also, note that the resampling step is the simplest technique, but other techniques exist which can add less noise.

Sequential Monte Carlo samplers

As described in the previous section, SMC provides an efficient method for sampling from a sequence of related distributions, making it suitable for situations such as performing inference as observations are made. Del Moral et al. suggest using SMC even when only a single distribution is to be inferred. Whereas Liu et al. considered a sequence of $\pi _{t}(\mathbf {x} _{t})$ (and focused on the case where $\mathbf {x} _{t}$increases in dimension with $t$), Del Moral et al. focus on the case where we only wish to infer $\pi (\mathbf {x} )$ (which we can evaluate down to a normalising constant). To allow SMC to be used, an artificial sequence of $p$ distributions is constructed from $\pi (\mathbf {x} )$. We shall denote this sequence as $\pi _{1}(\mathbf {x} ),\ldots ,\pi _{p}(\mathbf {x} )$. Note that $\mathbf {x}$ is now constant in dimension across all distributions. Several methods are discussed for creating this sequence:

1. If $\pi (\mathbf {x} )$ is a posterior distribution conditioned on $p$ variables, $y_{1,\ldots ,p}$, we can construct a series of distributions by setting $\pi _{1}(\mathbf {x} )$ to the prior and then conditioning each distribution on one more observation than the previous: $\pi _{n+1}(\mathbf {x} )=p(\mathbf {x} |y_{1,\ldots ,n})$ for $n=1,\ldots ,p$
2. In order to sample from an intractable distribution, $\pi (\mathbf {x} )$, we could construct a "smooth" path to it from $\mu (\mathbf {x} )$, some tractable initial distribution, as follows: $\pi _{n}(\mathbf {x} )=\mu (x)^{1-\phi _{n}(x)}\pi (x)^{\phi _{n}(x)}$where $0\leq \phi _{1}\leq \ldots \leq \phi _{p}=1$

Other sequences are suggested if we wish to find the mode of $\pi (\mathbf {x} )$ or estimate the probability of a rare event, rather than sampling from it:

1. To converge to the mode of $\pi (\mathbf {x} )$, we can use the sequence $\pi _{n}=\pi (x)^{\phi _{n}}$ where $\phi _{n}\rightarrow \infty$ as $n$ increases
2. If $\pi (\mathbf {x} )$ is simple to sample from and we wish to estimate $\pi (A)=\mathbb {P} (\mathbf {x} \in A)$ for a rare event, $A$, such that $\pi (A)\approx 0$, we can use the sequence $\pi _{n}(\mathbf {x} )=\pi (\mathbf {x} )\mathbf {1} _{A_{n}}(\mathbf {x} )$ where $A_{1}\supset \ldots \supset A_{n}=A$, and the estimate of $\pi (A)$ is obtained form the normalising constant during inference.

Transition Kernels

As in Liu et al., SMC samplers requires some kernel to transition from a sample from $\pi _{n}(\mathbf {x} )$ to a sample from $\pi _{n+1}(\mathbf {x} )$. In Liu et al., we discussed kernels consisting of SIS and rejection sampling. However, these were only suitable for the case when the dimension of $\mathbf {x}$ was increasing. Since the dimension of $\mathbf {x}$ is constant in the case of SMC samplers, these are no longer suitable. For the discussion of this paper, we denote the transition kernel by $K_{n}(\mathbf {x} ,\mathbf {x} ')$, which transitions from $\mathbf {x}$, a sample from $\pi _{n-1}$, to $\mathbf {x'}$, a sample from $\pi _{n}$. In order to minimise the variance of the weights, the optimal transition would be $K_{n}(\mathbf {x} ,\mathbf {x} ')=\pi _{n}(\mathbf {x} ')$. However, this is typically not possible and so it is common to use a local move from $\mathbf {x}$: Del Moral et al. suggest using either a Markov chain Monte Carlo (MCMC) transition or an approximate Gibbs step. We will now describe these options in more detail.

Markov chain Monte Carlo:

$K_{n}(\mathbf {x} ,\mathbf {x} ')$ may be set to be an MCMC kernel with $\pi _{n}(\mathbf {x} ')$ as a stationary distribution. This means that, if $\mathbf {x}$ is a sample from $\pi _{n}(\mathbf {x} )$, $\mathbf {x} '$ will also be a sample from $\pi _{n}(\mathbf {x} )$ (although it will not be independent of $\mathbf {x}$). Since $\mathbf {x}$ is a sample from $\pi _{n-1}(\mathbf {x} )$ instead, $\mathbf {x} '$ will not be a perfect sample from $\pi _{n}(\mathbf {x} )$, but will be "close" if $\pi _{n}$ is close to $\pi _{n-1}$ (and $\mathbf {x} '$ will be "closer" to $\pi _{n}$ than $\mathbf {x}$ was). The discrepancy is then accounted for in the weights (which are discussed later). A vast literature exists for constructing good MCMC kernels for arbitrary stationary distributions , and so this choice can take advantage of this literature to design an efficient kernel, and so approximate $\pi _{n}(\mathbf {x} )$ well.

Approximate Gibbs Step

This corresponds to the SIS and rejection sampling methods discussed above, where one variable (e.g. $x_{n}$) is sampled conditioned on all the other variables (e.g. $\mathbf {x} _{n-1}$). This is an exact Gibbs step if $x_{n}$ is sampled from $p(x_{t+1}|\mathbf {y} _{t+1}\mathbf {x} _{t}^{(j)})$ exactly, or an approximate step is SIS or an MCMC approach is used. Note that, since the dimension of $\mathbf {x}$ is constant in an SMC sampler, the sampled value will overwrite a previous value, rather than being appended to $\mathbf {x} _{n-1}$ as when this was discussed above.

Weights

Ideally, it would be possible to calculate weights for each $\mathbf {x} _{n}^{(j)}$ as ${\frac {\pi _{n}(\mathbf {x} _{n}^{(j)})}{\gamma _{n}(\mathbf {x} _{n}^{(j)})}}$, where $\gamma _{n}(\mathbf {x} _{n}^{(j)})$ is the density function of the "proposal" distribution for $\mathbf {x} _{n}^{(j)}$ given by:

$\gamma _{n}(\mathbf {x} _{n}^{(j)})=\int \gamma _{1}(\mathbf {x} _{1}^{(j)})\prod _{n'=2}^{n}K_{n'}(\mathbf {x} _{n'-1}^{(j)},\mathbf {x} _{n'}^{(j)})d\mathbf {x} _{1},\ldots ,d\mathbf {x} _{n-1}$.

However, the marginalisation over $\mathbf {x} _{1},\ldots ,\mathbf {x} _{n-1}$ is typically intractable and so it is impossible to calculate the weights in this form. To avoid this marginalisation, we consider the samples as samples of $\mathbf {x} _{1,\ldots ,n}$. from $\pi _{1,\ldots ,n}$, rather than just $\mathbf {x} _{n}$ from $\pi _{n}$. This means that the weights are given by:

$w_{1:n}^{(j)}={\frac {\pi _{1,\ldots ,n}(\mathbf {x} _{1,\ldots ,n}^{(j)})}{\eta _{1,\ldots ,n}(\mathbf {x} _{1,\ldots ,n}^{(j)})}}$,

where $\eta _{1,\ldots ,n}(\mathbf {x} _{1,\ldots ,n}^{(j)})=\eta _{1}(\mathbf {x} _{1}^{(j)})\prod _{n'=2}^{n}K_{n'}(\mathbf {x} _{n'-1}^{(j)},\mathbf {x} _{n'}^{(j)})$. Since $\pi _{n}$ is a marginal of $\pi _{1,\ldots ,n}$, sampling from the joint distribution will yield samples from $\pi _{n}$. Since it is not always possible to evaulate $K_{n}(\mathbf {x} _{n-1},\mathbf {x} _{n})$ (for example, the rejection probability of an MCMC kernel is usually intractable), Del Moral et al. propose to express $\eta _{1,\ldots ,n}(\mathbf {x} _{1,\ldots ,n}^{(j)})$ using artificial "backward" kernels $L_{n}$, such that:

$\eta _{1,\ldots ,n}(\mathbf {x} _{1,\ldots ,n}^{(j)})=\eta _{n}(\mathbf {x} _{n}^{(j)})\prod _{n'=2}^{n}L_{n'}(\mathbf {x} _{n'-1}^{(j)},\mathbf {x} _{n'}^{(j)})$

Computing $L_{n}$ such that the above equation gives the "true" proposal distribution is still not necessarily tractable. However, Del Moral et al. use the insight that the weights are only required to define an asymptotically correct marginal distribution over $\mathbf {x} _{n}$, and the distribution over $\mathbf {x} _{1:n-1}$ is unimportant. This means that $L_{n}$ need not be an accurate inverse of the transitions for the sampler to yield asymptotically consistent estimates (although the choice of $L_{n}$ will affect the variance of the weights). For a given $L_{n}$, the weights of particle $j$ is updated after step $n$ by multiplying by:

${\tilde {w}}(\mathbf {x} _{n-1}^{(j)},\mathbf {x} _{n}^{(j)})={\frac {\gamma _{n}(\mathbf {x} _{n}^{(j)})L_{n-1}(\mathbf {x} _{n-1}^{(j)},\mathbf {x} _{n}^{(j)})}{\gamma _{n-1}(\mathbf {x} _{n-1}^{(j)})K_{n}(\mathbf {x} _{n-1}^{(j)},\mathbf {x} _{n}^{(j)})}}$

where $\gamma _{i}(\mathbf {x} )$ is a (potentially) unnormalised density function of the $i$th distribution, $\pi _{i}(\mathbf {x} )$. Depending on the choice of $L_{n-1}(\mathbf {x} _{n-1},\mathbf {x} _{n})$, $K_{n}(\mathbf {x} _{n-1},\mathbf {x} _{n})$ may cancel out and so this technique can be used even when $K_{n}(\mathbf {x} _{n-1},\mathbf {x} _{n})$ is intractable.

Algorithm

The algorithm for the SMC sampler is defined as follows:

1. Define $L_{n-1}(x_{n-1},x_{n})$, $K_{n}(x_{n-1},x_{n})$ sequence of $\pi _{1}(\mathbf {x} ),\ldots ,\pi _{p}(\mathbf {x} )$
2. Initialization:
   • set $n=1$
   • for $j=1,\ldots ,m$ draw $\mathbf {x} _{1}^{(j)}\sim \eta _{1}$
   • set $\{w_{1}^{(j)}\}=\{{\frac {\pi _{1}(\mathbf {x} _{1}^{(j)})}{\eta _{1}(\mathbf {x} _{1}^{(j)})}}\}$ and normalise to obtain $\{W_{1}^{(j)}\}$
   • iterate through steps 2 and 4
3. Resampling:
   • If effective sample size less than some threshold, resample and set $W_{n}^{(j)}={\frac {1}{m}}$
4. Sampling:
   • set $n=n+1$
   • if $n>p$, stop
   • for $j=1,\ldots ,m$ draw $\mathbf {x} _{n}^{(j)}\sim K_{n}(\mathbf {x} _{n-1}^{(j)},.)$
   • set $\{w_{n}^{(j)}\}=\{W_{n-1}^{(j)}{\frac {\gamma _{n}(\mathbf {x} _{n}^{(j)})L_{n-1}(\mathbf {x} _{n-1}^{(j)},\mathbf {x} _{n}^{(j)})}{\gamma _{n-1}(\mathbf {x} _{n-1}^{(j)})K_{n}(\mathbf {x} _{n-1}^{(j)},\mathbf {x} _{n}^{(j)})}}\}$ and normalise to obtain $\{W_{n}^{(j)}\}$

My Thoughts on the Contributions

Del Moral et al. suggest that SMC can be used even when inference is only required over a single distribution, rather than a series of distributions. A very powerful application of this is for inference in well-defined models, in addition to the likes of Metropolis-Hastings and Gibbs sampling, and SMC can sometimes provide significantly improved inference compared to these methods. One potential issue is that it is not necessarily clear how best to select parameters such as the transition kernels, and how to define the sequence of distributions (which the quality of inference is presumably heavily dependent on). Despite this, I believe that the introduction of the framework is a valuable contribution.