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?
