3
$\begingroup$

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.

$\endgroup$
2
$\begingroup$

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.

Formulae

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.

Empirical Stuff

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.

enter image description here

| cite | improve this answer | |
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.