Jason Blevins’s Notebook
Doucet de Freitas and Gordon (2001)

Citation

These notes are based on the following article:

Doucet, Arnaud, Nando de Frietas, and Neil Gordon (2001): “An Introduction to Sequential Monte Carlo Methods,” in Sequential Monte Carlo Methods in Practice, ed. by A. Doucet, N. de Freitas, and N. Gordon. New York: Springer-Verlag.

Introduction

Real-world data analysis often requires estimating unknown quantities given only a sequence of observations on some related observable quantity. In a Bayesian framework, one typically has a priori knowledge of the model: a prior distribution of the unobservable quantity of interest and likelihood functions which relate the observables to the unobservables. The resulting posterior distribution of the unobservables can be calculated using Bayes’ theorem. This allows one to conduct inference about the unobserved quantities.

In some cases, it is natural to process observations sequentially. These cases are the focus of this introductory article and the rest of the book. For example, on-line applications such as radar tracking or estimating financial instruments where new data become available in real time, it is easier to update the previously formed posterior distribution than to recalculate it from scratch.

Sequential Monte Carlo methods are simulation-based methods for calculating approximations to posterior distributions. They avoid making linearity or normality assumptions required by related methods such as the Kalman filter.

Model

Let $\{x_t{\}}_{t=0}^{\mathrm{\infty }}$ denote the sequence of unobserved states, with $x_t\in 𝒳$. The authors consider only the case when $\{x_t\}$ is a Markov process with initial distribution $p(x_0)$. The Markov assumption is not needed and is only for tractability. Let $p(x_t\mid x_{t-1})$ denote the transition equation. Let $\{y_t{\}}_{t=1}^{\mathrm{\infty }}$ denote a sequence of observations with $y_t\in 𝒴$. Suppose that, given $\{x_t\}$, the observations are conditionally independent. Let $p(y_t\mid x_t)$ denote the marginal distribution, or the observation equation.

Thus, the model is defined by the following distributions:

Define $x_{0:t}\equiv \{x_0,\dots ,x_t\}$ and $y_{0:t}\equiv \{y_1,\dots ,y_t\}$. Our goal is to estimate the posterior distribution $p(x_{0:t}\mid y_{1:t})$ recursively. We also may care about the marginal distribution $p(x_t\mid y_{1:t})$ or an expectation such as

for some $h_t:{𝒳}^t\to {ℝ}^n$ that is integrable with respect to $p(x_{0:t}\mid y_{1:t})$.

Bayes’ theorem gives us the posterior distribution at any time $t$:

There is an inherent recursive relationship for updating the posterior distribution. We can express $p(x_{0:t+1}\mid y_{0:t+1})$ in terms of $p(x_{0:t}\mid y_{1:t})$:

where the last equality uses the definition of the conditional density, the conditional independence assumption, and the Markov assumption. The expression in the denominator is constant with respect to $x_{0:t+1}$.

Perfect Monte Carlo Sampling

In the ideal situation, we could draw a sample of size $N$ from $p(x_{0:t}\mid y_{1:t})$, the posterior distribution of interest. Call this sample $\{x_{0:t}^{(i)}{\}}_{i=1}^N$. We can form an estimate of $p(x_{0:t}\mid y_{1:t})$ using the empirical distribution of the sample:

From here, we can estimate the integral $I(h_t)$ by

Typically it is impossible to get such a sample since $p(x_{0:t}\mid y_{1:t})$ is multivariate, known only up to a constant of proportionality, and non-standard. Markov chain Monte Carlo methods may be used in similar situations, but they are not well-suited to recursive problems such as the ones considered here.

Importance Sampling

Since we can’t draw from $p(x_{0:t}\mid y_{1:t})$ directly, we can use importance sampling to draw from $p(x_{0:t}\mid y_{1:t})$ indirectly using an importance sampling density $\pi (x_{0:t}\mid y_{1:t})$ and a corresponding importance weight for each draw. The importance sampling density must be chosen so that its support includes the support of $p(x_{0:t}\mid y_{1:t})$. We have

where

is the importance weight.

Thus, to approximate the expectation $I(h_t)$, we simply need to draw a sample $\{x_{0:t}^{(i)}{\}}_{i=1}^N$ from $\pi (x_{0:t}\mid y_{1:t})$ and calculate

where

are the normalized importance weights.

This procedure can be interpreted as a sampling method with an approximate posterior given by

We can then approximate $I(h_t)$ by

Although this solves the problem of being able to obtain draws, this approach is still not well-suited for our recursive problem.

Sequential Importance Sampling (SIS)

The goal of Sequential Importance Sampling (SIS) is to compute an estimate of $p(x_{0:t}\mid y_{1:t})$ using the past simulated values $\{x_{0:t-1}^{(i)}{\}}_{i=1}^N$. We have the following recursive representation of the importance sampling distribution:

This recursive relationship allows us to calculate the importance weights recursively:

Hence,

Notice that the $p(y_t\mid y_{1:t-1})$ term is constant for all $x_{0:t}$ and cancels out.

A special case occurs when

so that

and

SIS is usually inefficient for high-dimensional integrals because as $t\to \mathrm{\infty }$, the importance weights for some particles quickly approaches zero.

The Bootstrap Filter

To prevent particle degeneracy, the bootstrap filter introduces a resampling step which eliminates particles with low importance weights. Instead of using the importance-weighted empirical distribution ${\hat{P}}_N$ we use a uniformly-weighted distribution

where $N_t^{(i)}$ is the number of offspring of the particle $x_{0:t}^{(i)}$ which is to be determined by some branching mechanism. The most common such mechanism involves resampling $N$ times from ${\hat{P}}_N$ (Gordon et al., 1993). We require ${\sum }_i{N}_t^{(i)}=N$ for all $t$. If $N_t^{(j)}=0$, the particle $x_{0:t}^{(j)}$ dies. We want to choose $N_t^{(i)}$ such that

The surviving particles, those with $N_t^{(i)}>0$, are approximately distributed according to $p(x_{0:t}\mid y_{1:t})$.

Algorithm

1. Initialization ($t=0$).

• For $i=1,\dots ,N$, draw $x_0^{(i)}\sim p(x_0)$.
• Set $t\leftarrow 1$.

1. Importance Sampling Step

• For $i=1,\dots ,N$, draw ${\tilde{x}}_t^{(i)}\sim p(x_t\mid x_{t-1}^{(i)})$ and set ${\tilde{x}}_{0:t}^{(i)}\leftarrow (x_{0:t-1}^{(i)},{\tilde{x}}_t^{(i)})$.

• For $i=1,\dots ,N$, calculate the importance weights ${\tilde{w}}_t^{(i)}=p(y_t\mid {\tilde{x}}_t^{(i)})$ and normalize them.

1. Resampling Step

• Take $N$ draws $\{x_{0:t}^{(i)}{\}}_{i=1}^N$ with replacement from the set $\{{\tilde{x}}_{0:t}^{(i)}{\}}_{i=1}^N$ with weights $\{{\tilde{w}}_t^{(i)}{\}}_{i=1}^N$.

• Set $t\leftarrow t+1$ and go to step 2.