# Particle filter (sequential Monte Carlo) for a non-Gaussian hierarchical model

I have the following, which I am attempting to model with a particle filter.

\begin{align*} y_{i,t}&\sim\mathrm{Poisson}\left(\lambda_{i,j,t}\right)\\ y_{j,t}&\sim\mathrm{Poisson}\left(\mu_{i,j,t}\right) \end{align*}

\begin{align*} \lambda_{i,j,t}&=\delta_1\exp\left(\alpha_{i,t}-\alpha_{j,t}\right)\\ \mu_{i,j,t}&=\delta_2\exp\left(\alpha_{j,t}-\alpha_{i,t}\right) \end{align*}

\begin{align*} \alpha_{i,t}&\sim\mathrm{N}\left(\alpha_{i,t-1},\sigma_{\alpha}^2\right)\\ \alpha_{j,t}&\sim\mathrm{N}\left(\alpha_{j,t-1},\sigma_{\alpha}^2\right)\\ \end{align*}

Where:

t is a time index, $t=1,\ldots,5320$;

$i,j\in\left\lbrace 1,\ldots,40\right\rbrace$, $i\ne j$ (so 40 states altogether but only 2 update at each time $t$);

the $\alpha$ states are unobserved;

and the $y$'s are observed for every $t$.

For simplicity, assume the values of $\delta_1$, $\delta_2$, and $\sigma_{\alpha}^2$ are both known and constant.

My idea was to have a set of particles for each state. Then update them (the 2 sets involved at time $t$) with Gaussian noise, using the equations for the states. Given these values, the lambdas can be calculated and then the conditional distributions for $y_{i,t}$ and $y_{j,t}$, given the states, are known. However, on calculating the importance weights from the conditional distributions $y_{i,t}|\alpha_{i,t},\alpha_{j,t}$ and $y_{j,t}|\alpha_{j,t},\alpha_{i,t}$, you get two sets of importance weights from \begin{align*} y_{i,t}&|(\alpha_{i,t}-\alpha_{j,t})\\ y_{j,t}&|(\alpha_{j,t}-\alpha_{i,t}). \end{align*}

My question is how then to resample particles for $\alpha_{i,t}$ and $\alpha_{j,t}$, or possibly how can I do this differently to overcome the problem. I am fairly new to sequential Monte Carlo methods, so any other advice or suggestions on tackling the problem would be appreciated.

I will give you two references of possible interest and an alternative method to your (non trivial) problem.

References:

An alternative for sampling from your model consists of using a Metropolis within Gibbs sampler. There is a very nice R package which implements an adaptive version of this method which doesn't require tunning parameters: spBayes. The only thing you have to implement is the joint posterior (parameters + random effects), which should be relatively easy. I have tried it, and it works very well. It might be slow in high dimensions, though.

• The references were helpful, thanks. The first reference is good for the idea of multiple states but not so much for dealing with the specific problem I have of only certain states updating. For the second reference and your alternative suggestion, I am hoping for a filter that I can code myself and I need to be able to explain the steps involved as well as get the results, but they are interesting alternative strategies. Also, note that I have modified my original problem with a slight but acceptable simplification after much thought, it may make the problem easier to deal with.
– MLaz
Jul 25, 2014 at 16:53

Your model does not make sense. A model consists of a FULL state transition distribution, and a FULL observation distribution (you can't ignore some of the observation vector at any time). Also, you are not saying how $i$ and $j$ are chosen at each time; I assume it's non-randomly. Also you might have superfluous subscripts for $\mu$ and $\lambda$.

Your full transition distribution for $(\alpha_{1,t},\ldots, \alpha_{40,t})'$ given $(\alpha_{1,t-1},\ldots, \alpha_{40,t-1})'$ could be

$$p(\mathbf{\alpha_{t}}|\mathbf{\alpha_{t-1}}) = \mathcal{N}(\alpha_{i,t};\alpha_{i,t-1},\sigma^2_a)\mathcal{N}(\alpha_{j,t};\alpha_{j,t-1},\sigma^2_a)\prod_{k\neq i,j}1(\alpha_{k,t} = \alpha_{k,t-1}),$$ assuming WLOG that $i < j$. Not updating the parts that aren't changing is the same as updating with respect to that distribution.

Then you need a full observation density. It does not make sense to ignore all but two elements of the obserations. It could be something like this: $$p(\mathbf{y_{t}}|\mathbf{\alpha_{t}}) = \text{Poisson}(y_{i,t};\lambda_{i,t})\text{Poisson}(y_{j,t};\mu_{j,t}) \prod_{k\neq i,j}\text{(missing stuff here)}.$$

Once you have these two things fully defined, every particle has a weight, and resampling is then done the same way it always is. I cannot outline a full algorithm because you have told us certain need-to-know details about how you are carrying all of this out. But I suspect that you're running into trouble about how to implement this because you do not fully understand what a state space model is.

Edit: didn't realize this was asked so long ago...