# Does the Auxiliary Particle Filter not ordinarily yield estimates for $p(y_t|y_{1:t-1})$

## Background

Say you're implementing basic importance sampling. To approximate a distribution $$p(x)$$, you utilize an importance density $$q(x)$$ that you can sample from. Also, assume that you can only evaluate $$p_u(x) = C_p p(x)$$, the target density up to a constant of proportionality.

So you sample $$X^i \sim q(x)$$ for $$i=1,\ldots,n$$ and you weight each of these samples with $$W^i = p_u(X^i)/q(X^i)$$, where $$p_u(x) = C_p p(x)$$. By the law of large numbers, a consistent estimator of $$C_p$$ is $$\hat{C}_p = \frac{1}{N}\sum_i W^i$$.

However when you can't evaluate $$q(x)$$ exactly either, but instead can only evaluate $$q_u(x) = C_q q(x)$$, your weights become $$W^i_a = p_u(X^i)/q_u(X^i)$$, and your average weight is no longer a consistent estimate of $$C_p$$

## Auxiliary Particle Filter decomposition

It seems to me that with auxiliary particle filters, the weight updates are sort of like this second situation. Here's a hopefully useful decomposition of the time $$t$$ smoothed distribution.

\begin{align*} p(x_{1:t},k|y_{1:t}) &= C^{-1} \frac{g(y_t|x_t)f(x_t|x_{t-1},k)p(k|x_{1:t-1},y_{1:t-1}) }{q_t(x_t,k|y_{t},x_{t-1})} \\ & \hspace{20mm} \times \frac{p(x_{r+1:t-1},y_{r+:t-1}|x_r)}{q(x_{r+1:t-1}|y_{r+1:t-1})} \\ & \hspace{20mm} \times q_t(x_t,k|y_{1:t},x_{1:t-1}) q(x_{r+1:t-1}|y_{r+1:t-1})p(x_{1:r}|y_{1:r}) \end{align*}

Notation: this is a state space model with observations $$y_t$$, and hidden states $$x_t$$. $$g(y_t|x_t)$$ is the observation density, $$f(x_t|x_{t-1})$$ is the state transition density, $$p(k|x_{1:t-1},y_{1:t-1})$$ is the mass function of the previous time's particle index. The time $$r$$ is the time of the last resampling step. And $$q_t(x_t,k|y_{1:t},x_{t-1})$$ is your proposal distribution. Sorry that this notation isn't standard.

## My question

This $$q_t(x_t,k|y_{t},x_{t-1})$$, it seems like everything I read describes an implementation where the user is not not evaluating this exactly, just up to proportionality. Right now I have

$$q^{\text{apf}}_t(x_t,k|y_{1:t},x_{1:t-1}) = \frac{f(x_t|x_{t-1},k)g(y_t|\mu_t[x_{t-1},k])p(k|x_{1:t-1}, y_{1:t-1})}{\sum_{k'} g(y_t|\mu_t[x_{t-1},k'])p(k'|x_{1:t-1}, y_{1:t-1})},$$ exactly. So the weight update at time $$t$$ for particle $$i$$ $$\frac{g(y_t|X^i_t)f(X_t^i|x_{t-1}^i,k^i)p(k^i|x_{1:t-1},y_{1:t-1}) }{q_t(x_t^i,k^i|y_{t},x^i_{t-1})}$$ simplifies to $$g(y_t|X_t^i)/g(y_t|\mu_t[x_{t-1}^i, k^i])$$ if you use the non-standardized evaluation of $$q_t^{\text{apf}}$$. For example wiki seems to have it this way. . These are the so-called second stage weights from Pitt and Shephard's original paper. But it simplifies to

$$\frac{ g(y_t|X^i_t)}{g(y_t|\mu_t[x^i_{t-1},k^i])} \sum_{k'} g(y_t|\mu_t[x_{t-1},k'])p(k'|x_{1:t-1}, y_{1:t-1}) ,$$ if I use the exact $$q_t$$ and spend that extra time normalizing it.

If I want estimates of $$p(y_t|y_{1:t-1})$$ at every time I should use the latter, correct? Since there are only two possible answers, I could code this up and find out that way, but if anyone wants to weigh in on this in the meantime, I'm all ears.

I was right about having to normalize the function that evaluates the importance/proposal density. At least if you want to use the usual formula to estimate at each time point $p(y_t|y_{1:t-1})$. For what it's worth, here's an answer for other people to reference.
Evaluate the normalized importance density so you can get these weight updates: $$\alpha_t^i \overset{\text{def}}{=} \frac{ g(y_t|X^i_t)}{g(y_t|\mu_t[x_{t-1},k^i])}\sum_{k} g(y_t|\mu_t[x_{t-1},k])p(k|x_{1:t-1}, y_{1:t-1}) ,$$ then use that to get estimates of the conditional likelihoods with the formula $$\hat{p}(y_t|y_{1:t-1}) = \sum_i \widetilde{W}^i_{r+1:t-1}\alpha_t^i.$$ where $\widetilde{W}^i_{r+1:t-1}$ are the normalized particle weights you had sitting in memory before you updated.
The proof pretty much follows the first section of the question. Also I tried coding both on simulated 200 observations from the following model: $$X_t = .91\cdot X_{t-1} + W_t$$ and $$Y_t = X_t + V_t$$ with $\{W_t\} \overset{iid}{\sim} \text{Normal}(0, 1)$ and $\{V_t\} \overset{iid}{\sim} \text{Normal}(0, 1.5^2)$ and $X_1 \sim \text{Normal}(0, \frac{1}{1-.91^2})$. I ran an APF, a regular bootstrap version of SISR, and a Kalman filter, the simulation-based filters each with 1500 particles and resampling at every time point. Here's a plot showing that the log of the conditional likelihoods all line up.