# Sequential Monte Carlo (particle filter) with Metropolis-Hastings weighting

Let's say we are interested in approximating the following expectation:

$$\mathbb{E}[h(x)] = \int h(x)\pi(x) dx$$

Where $h(x)$ is an arbitrary function and $\pi(x)$ is a distribution known only up to a normalizing constant. We could approximate this expectation using a Metropolis-Hastings sampler to draw samples $\{x_t\}_{t=1}^N$ from a Markov chain with stationary distribution $\pi(x)$, we first define the MH weight function:

$$w(x,x') = \frac{\pi(x')T(x',x)}{\pi(x)T(x,x')}$$

• sample $x_0$ from arbitrary proposal
• for $t = 1...N$
• propose $x'$ by sampling from $T(x_t,\cdot)$
• if $w(x_t,x') > r$ (where $r \sim \mathcal{U}(0,1)$) then $x_{t+1} = x'$ else $x_{t+1} = x_t$
• $S \leftarrow S + h(x_{t+1})$

Giving the final approximation: $$\mathbb{E}[h(x)] \approx \frac{S}{N}$$

The question is, can we do a the same 'trick' (sampling from a Markov chain with stationary distribution $\pi$) using sequential Monte Carlo (particle filter) and use weights instead of a simple accept/reject rule. The comparison, in theory, is similar to that of ordinary rejection sampling vs. importance sampling. To do this, consider the following procedure: first sample a set of N particles: $\{x_0^{(i)}\}_{i=1}^N$ from arbitrary proposal, then execute the following loop L times:

• for each $x_t^{(i)}$ propose $x_{t+1}^{(i)}$ by sampling from $T(x_t^{(i)},\cdot)$
• weight each $x^{(i)}$ by $w(x_t^{(i)},x_{t+1}^{(i)})$
• $S \leftarrow S + \sum_{i=1}^N w(x_t^{(i)},x_{t+1}^{(i)})h(x_{t+1}^{(i)})$
• resample N new particles $x_{t+1}^{(i)}$ proportional to weights

Giving the final approximation: $$\mathbb{E}[h(x)] \approx \frac{S}{N L}$$

This procedure is different than typical sequential importance sampling, or adaptive importance sampling, techniques in a couple of ways. First of all, the proposal distribution ($T(\cdot,\cdot)$) does not change as new samples arrive, also, it is not a proposal that is attempting to be anywhere close to the optimal IS distribution $\frac{|h(x)|\pi(x)}{Z}$. Instead, the proposal is better seen as a Markov transition kernel and must be combined with the appropriate weighting function to ensure that the transitions are consistent with the desired Markov chain. In fact, this procedure is very similar (but not the same) to running $N$ independent Markov chains.

The question comes down to a couple issues. First, is this method even correct in it's stated goals, that is: if particles $\{x_t^{(i)}\}_{i=1}^N$ are distributed according to $\pi(\cdot)$, after one loop of the above procedure, will $\{x_{t+1}^{(i)}\}_{i=1}^N$ be distributed according to $\pi(\cdot)$? Also, if this is true, does this approach offer any advantages over the basic Metropolis-Hastings sampler presented first? Is anyone able to point toward any papers relating specifically to this idea: SMC to sample from stationary Markov chains? I suppose this is technically a MCMC method so, is there such thing as "weighted MCMC" or "weighted Metropolis"? Or is this similar to the idea of "waste recycling" for MCMC methods (I'm not very familiar with this)?

Later edits:

On further thought it looks like the sequential samples are not distributed according to $\pi(\cdot)$ because the weighting function does not meet the detailed balance condition unless $\pi(x)$ happens to be uniform. This can be seen by considering symmetric transition function $T(x,x')$

$$\begin{eqnarray*} \pi(x)w(x,x') = \pi(x')w(x'x) \\ \pi(x)\frac{\pi(x')}{\pi(x)} = \pi(x')\frac{\pi(x)}{\pi(x')} \\ \pi(x') = \pi(x) \end{eqnarray*}$$

This detail could be fixed by using the exact weighting function $\min(1,w(x,x'))$ and ensuring that each of the previous states is included in the weighted set with weight equal to $1-w(x,x')$. If this is done the method is very similar to just running $N$ Markov chains in parallel. These chains would not be independent per say, as there would be some complex interactions between chains (due to the combined resampling stage); I'm not able to see what the effects of this would be, if any. However, there still might be some merit to doing this because the weighted samples could still be used in the computation of the final expectation. That is, it seems this method still might achieve the goal of a rejection-free MCMC sampler.

More edits:

I should make another distinction clear. The generalized importance sampling method covered in 14.2 of Monte Carlo Statistical Methods appears to be very close to the SMC procedure I wrote out. However, they are actually quite different. First of all in this method, as with MH, the support for the proposal distribution $T(x,x')$ need not contain the support of the target proposal $\pi(x)$ (as required by Lemma 14.1), rather the support requirements fall out of the ergodicity of the simulated Markov chain. This makes $T(x,x')$ significantly different than the combination of Kernel function $K(x,x')$ and proposal $g(x)$ found in GIS.

The population Monte Carlo method, however, seems to subsume the procedure I've mentioned here. The difference being I have not adapted my proposal over time where a PMC method allows a framework for doing so (among other potential advantages). So this looks to be the answer to my question, Yes you can do this and in fact, you can be even more general about your choice of transition kernel (proposal function) than you can with MH.

If this is correct, I'm wondering why Metropolis-Hastings and related approaches remain so popular if they can simply be subsumed by the more general PMC method?

• This is not a Metropolis-Hastings scheme given that you resample and weight at each step $t$. It is thus validated by a sequential importance sampling argument.. See, e.g., Del Moral, Doucet and Jasra, 2006. – Xi'an Dec 9 '11 at 13:12
• For Metropolis Hastings you are missing the acception/rejection step, which is necessary to ensure that the samples are distributed according to $\pi$ – Thies Heidecke Dec 10 '11 at 10:47
• @Theis: Yes, but my question is whether you can do SIS to also get points distributed according to $\pi$. In my above procedure I would expect that after some 'burn in' time the particles would be drawn from the stationary distribution equal to $\pi$. The difference then, between Metropolis, is that the expectation is given by a weighted distribution rather than an accept/reject distribution (likely with many duplicates). This difference would be in some ways similar to vanilla importance sampling vs. vanilla rejection sampling. – fairidox Dec 10 '11 at 21:08

The answer is in the title: this is called sequential Monte Carlo or particle filtering or population Monte Carlo. It is validated in wider generality as an iterated importance sampling scheme where each importance sample is used to generate the following sample. This is for instance covered in Chapter 14 of our book Monte Carlo Statistical Methods.

The specific issue of using the whole sequence of simulation is found in the literature, but not in the direct way you propose: using all samples at once with the same weights does not behave nicely when some of the weights are huge (as when one starts with a poor guess). This is covered in the fantastic multiple mixture paper by Owen and Zhou (2000, JASA) and in our more recent adaptive version (when $T$ depends on the iteration $t$ and on the past simulations) called AMIS.

• Thanks for the answer, I know this is just SIS, I'm really just asking if SIS can be used to create a no-rejection MH sampler. Essentially, by using weights instead of a simple accept/reject. The difference between regular SIS is that it doesn't use the whole sequence of points but the set of weighted points at each step (since they are presumably from the stationary distribution). Either way, I implemented this just to see, it appears that the expectation is correct. – fairidox Dec 10 '11 at 20:37
• However, as you point out, in my implementation some of the weights ended up being fantastically large even for relatively straightforward distributions. Consequently, I could not get the algorithm to provide a better estimate than vanilla Metropolis-Hastings. – fairidox Dec 10 '11 at 20:42
• I think you've answered my question in Chapter 14 of your book. So thank you again. I've edited my response (a couple of times) as it looks like what I've described is a less general version of population Monte Carlo. – fairidox Dec 11 '11 at 23:56
• @anonymous_4322: thank you for your remarks and further editing. To answer your final and most recent question as to why still MH if PMC is available, I would think PMC deteriorates as the dimension or the number of homogeneous groups of parameters (as in hierarchical models) increases. – Xi'an Dec 12 '11 at 9:57