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*}


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.


2 Answers 2


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


  1. Generalised linear mixed model analysis via sequential Monte Carlo sampling

  2. SMCTC : sequential Monte Carlo in C++

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.

  • $\begingroup$ 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. $\endgroup$
    – 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...


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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