Does this Monte Carlo Technique Have a Name?

Question

I sketched this algorithm out the other night. I am sure it has a name, I just do not know what it is yet. It would be helpful if someone could point me in the right direction for research. I provide motivation and then define the technique.

Motivation

I am in a situation where I have to recursively estimate Bayesian time series models. That is, I have to estimate a posterior for my data over a time periods $[1,2,...,t]$, then $[1,2,...,t,t+1]$, then $[1,2,...,t,t+1,t+2]$, going all the way toward the end of my sample.

I am using STAN to estimate posterior samples. Recursively estimating the posteriors through STAN is excessively burdensome. The motivation for the below MC technique is to speed up computation time by taking weighted re-samples from the posterior sample estimated over $[1,...,t]$ to approximate the posteriors for $[1,..,t+1]$,$[1,..,t+2]$, up to some time period $[1,..,t+h]$ where the weights are functions of the predictive likelihoods generated at time $t$.

The general logic is that the posterior for $[1,...,t]$ and $[1,...,t+h]$ should be "similar" if $t$ is large and $h$ is relativity small.

The MC Algorithm

Consider two conditional pdfs, $\pi(\theta|M_1)$ and $\pi(\theta|M_2)$. Assuming we have a sample $\theta_1, \theta_2,...,\theta_G$ where each $\theta_g$ was generated by one of the above distributions. The posterior odds ratio can be calculated as $$ \frac{\pi(\theta|M_1)}{\pi(\theta|M_2)} \times r $$ where $r=\frac{\pi(M_1)}{\pi(M_2)}$ is the prior odds ratio. Using this we know $$ Pr(\theta_g \sim \pi(\theta|M_1)) = w_g = 1-\frac{1}{\frac{\pi(\theta|M_1)}{\pi(\theta|M_2)}r +1} $$ We can then take a weighted sample of $\theta_1,...,\theta_G$, re-sampling each $\theta_g$ with probability $w_g$ to approximate $\pi(\theta|M_1)$.

Now consider the case where we new before hand that all the $\theta_g$ where sampled from $\pi(\theta|M_2)$. Then $w_g=r=0$. However, this also means that as $r \rightarrow 0$ the pdf of the weighted sample approaches $\pi(\theta|M_1)$. Thus, one could approximate $P(\theta|M_1)$ using a reasonably small $r$ (I have ran simulations for simple examples that seem to support this claim).

Application

The above example is trivial, the following application better exemplifies my motives for using this technique. Consider the time series data $y_1,...,y_T$. Let $Y_{n:m}=y_n,y_{n+1},...,y_{m}$, for $m>n$. Denote the non-normalized posterior distribution as $P(\theta|Y_{n:m})$, the likelihood as $f(Y_{n:m}|\theta)$, and the normalization constant (or marginal likelihood) as $m(Y_{n:m})$.

In this example I have a sample $\theta_1,...,\theta_G$ which is drawn exclusively from $P(\theta|Y_{1:t})$ and which I would like to re-sample to approximate $P(\theta|Y_{1:(t+h)})$. Following the procedure in the above section and using Corollary 1 and 2 given below I calculate the posterior odds ratio :

$$ \frac{P(\theta|Y_{1:(t+h)})/m(Y_{1:(t+h)})}{P(\theta|Y_{1:t})/m(Y_{1:t})} \times r=\frac{P(\theta|Y_{1:t})f(Y_{(t+1):(t+h)}|\theta)/m(Y_{1:(t+h)})}{P(\theta|Y_{1:t})/m(Y_{1:t})} \times r = $$ $$ \frac{f(Y_{(t+1):(t+h)}|\theta)}{E_{1:t}[ f(Y_{(t+1):(t+h)}|\theta)]} \times r $$

Which makes the re-sampling weights $$ w_g = 1-\frac{1}{\frac{f(Y_{(t+1):(t+h)}|\theta_g)}{E_{1:t}[ f(Y_{(t+1):(t+h)}|\theta)]}r + 1} $$

I am hoping that with a small enough choice of $r$, I can estimate the majority of posteriors in this fashion rather than re-estimating all the posteriors through STAN.

---

Corollary 1: $P(\theta|Y_{1:(t+h)}) = P(\theta|Y_{1:t})f(Y_{(t+1):(t+h)}|\theta)$.

Corollary 2: It can be shown that $\frac{m(Y_{(t+1):(t+h)})}{m(Y_{1:t})}=E_{1:t}[ f(Y_{(t+1):(t+h)}|\theta)]$, which is a joint predictive likelihood.

Answer 1

This is a sequential Monte Carlo problem where one considers a sequence of posteriors $\pi_1,\pi_2,\ldots,\pi_T$ that are "reasonably" close from one another for the simulations from $\pi_t$ to be interesting proposals for $\pi_{t+1}$.

The solution you propose is not valid as such because of the reliance on this proportion $r$ that is actually equal to zero when running the sequential Monte Carlo scheme. (Or, in other words, there is no prior ratio for the marginal distributions of $y_{1:t}$ and of $y_{1:t+1}$.) A very similar resolution is sequential importance sampling where the current sample from $\pi_t$ $\theta_1^t,...,\theta_G^t$ is reweighted by $$\{\pi_{t+1}/\pi_t\}(\theta_g^t)$$ to become a sample from $\pi_t$. In the case of a state-space model where $$\pi_{t+1}(\theta)\propto\pi_t(\theta) f(y_{t+1}|\theta,y_{1:t})$$ the importance weight reduces to $$\{\pi_{t+1}/\pi_t\}(\theta)\propto f(y_{t+1}|\theta,y_{1:t})$$

While valid per se (in the importance sampling sense that the expectation of the weighted average produces the correct expectation under $\pi_{t+1}$), this naïve type of sequential importance sampling is dominated by more advanced particle filters where the sample itself is modified before being reweighted. This avoids in particular the degeneracy phenomenon, namely that the weights$$\prod_{j=1}^h f(y_{t+j}|\theta,y_{1:(t+j-1)})$$degenerate to zero almost surely.

I suggest you read the Wikipedia entry on particle filters, as it is quite detailed and definitely rigorous, given that one of the leaders in the field, Pierre Del Moral, has contributed to its improvement.