Retrieve particles at a particular time

[hollosigergely]

I am experimenting with pfilter(), however, I cannot find a decent method to get the particles for a given timestamp. The saved.states() method returns an array of weights, but I do not know the coordinates of the particles in state space, only the weight. Or I might be wrong? As far as I know, a particle is a state space point (with it's coordinates) and a weight. I see the weights, but the coordinates are missing. I would like to watch the - approximate - distribution of the latent variable with the help of the particles, however, the weights alone are not enough for me.

[hollosigergely]

Okay, I might have found the solution, I misunderstood the algorithm. The returned states are the latent space coordinates of the particles after resampling, so the weights can be easily calculated from $f_{Y_n|X_n}(.)$. I hope I got it well.

[kingaa]

The particles retrieved by saved.states applied to a pfilterd.pomp object are stored after resampling. Therefore, they can be treated as equally weighted.

[hollosigergely]

Thank You! I suggest this to be emphasized in the manual, as well, however, after some insight one can easily infer this information.

[kingaa]

Good idea. I have added the line:

These are the unweighted particles, saved after resampling.

to the saved.states manual page. If you can think of a way of making this clearer, please suggest it!

By the way, I wonder whether it would be useful to save the weighted particles and their weights.... What do you think?

[hollosigergely]

Well, the theory is about particles at sampling points with weights, while resampling is a method to avoid particles with small weights, so I think it is better to stick to the theory, and provide the particles with weights. Nevertheless, this returned states also allow to reconstruct the weights (and also the states are minimal representations), but we need to apply the same equations as the pfilter function. If someone starts to use the library, she will expect to find the particles accompanied with the weights, I think (as well as me :) ).

[kingaa]

Interesting perspective. Can you tell me what you intend to do with the particles with weights?

[hollosigergely]

While the main purpose of the library is to select the proper hypothesis, I used the library to access the distribution of the latent variable, i.e. to approximate as $p(x|y_{1:t})\sim \sum_i w^{(i)} \delta (x-x^{(i)})$. That's all, it is just for get some insight of the behaviour of the latent variable to get intution.

[kingaa]

Just pushed some changes that make it easier to retrieve the weighted particles. At your leisure, have a look. Comments and suggestions (and of course pull requests) are welcome!

[hollosigergely]

Very well, I've checked it, but I do not get it, please help me understand. Here is the algorithm as in the PDF (to speak about the very same SMC).

In the unweighted case, we get the resampled $X^F_{n,j}$ particles with equal weights, that is OK.

In the weighted case, I would expect to have $X^P_{n,j}$ at each $n$ steps. Now I have $J$ particles at sample points $X^P_{n,j}$ as a result of line 3 in the algorithm. I would also expect to have the weights $w$ (or the normalized weights $\tilde{w}$) for each sample point as in line 4 or 5. This is clearly before resampling.

Now, I have two issues. First, weights do not add up to 1. As I see, the returned weights are log weights, but that does not work either (e.g. sum(exp(saved.states(pomp)$weights$'20')) results in other than 1). So weights are not normalized, i.e. they are from line 4 in the algorithm?
The second issue bothers me also. The states shall be the sample points $X^P_{n,j}$. However, I see a lot of particles with the same value. After sampling at line 3 in the algorithm, it has very low probability to have two particles with the same value.

Can You confirm this? Am I wrong somewhere?

[kingaa]

Good catch! I had checked the results, but not thoroughly enough. This was a bug, fixed with a one-character edit. I've put in some more units tests (in tests/pfilter2.R.

The weights returned, as you surmised, are the log unnormalized weights (per line 4). Since normalization introduces noise, I think it better to return the unnormalized weights. At any rate, they are normalized easily enough.

[hollosigergely]

Sorry to disturb you again, but still I have the problem of duplicated particles. I have the new version:

> packageVersion('pomp')
[1] ‘4.4.1.1’

However, running e.g. this test:

ou2() -> ou2
pf <- pfilter(ou2,Np=1000,save.states='weighted')
s<-saved.states(pf,format='d')
p<-s[s$variable=='x1' & s$time==10,]$value
plot(p)

I still have a lot of particles at the same sampling point at step 10. It still seems to be the resampled particles, I think. I got a plot like this:

Moreover, the sampled distribution of $x_1$ variable plots like this, which is wrong (of course, without normalization):

ou2() -> ou2
pf <- pfilter(ou2,Np=1000,save.states='weighted')
s<-saved.states(pf,format='d')
p<-s[s$variable=='x1' & s$time==10,]$value
w<-s[s$variable=='weight' & s$time==10,]$value
plot(p,exp(w))

Can You confirm this?

[kingaa]

I don't get the same thing. I see

library(pomp)
library(tidyverse)
theme_set(theme_bw())
stopifnot(packageVersion("pomp")=="4.4.1.1")


ou2() |>
  pfilter(Np=1000,save.states='weighted') -> pf
pf |>
  saved.states(pf,format='d') |>
  pivot_wider(names_from=variable) -> s

s |>
  filter(time==10) |>
  pull(x1) |>
  plot()

Some overlap is evident, but the particles are actually all distinct:

s |>
  select(-.id) |>
  group_by(time) |>
  distinct() |>
  summarize(n=n()) |>
  ungroup() |>
  filter(n != 1000) |>
  nrow()==0

# [1] TRUE

Plotting unnormalized weights against $x_1$, I see

s |>
  filter(time==10) |>
  ggplot(aes(x=x1,y=exp(.log.weight)))+
  geom_point()

Plotting the particles in the $x_1$- $x_2$ plane as dots with size proportional to their weights, I see the following. The contours reflect the density of the particles, without consideration of their weights. The red crosshairs show the location of the data.

s |>
  filter(time==10) |>
  ggplot(aes(x=x1,y=x2))+
  geom_density2d()+
  geom_point(aes(size=exp(.log.weight)),alpha=0.4)+
  geom_vline(
    xintercept=pf |>
      as_tibble() |>
      filter(time==10) |>
      pull(y1),
    color='red'
  )+
  geom_hline(
    yintercept=pf |>
      as_tibble() |>
      filter(time==10) |>
      pull(y2),
    color='red'
  )

I don't understand the discrepancies between what you and I see. Can you reproduce what I see?

[hollosigergely]

Well, I forgot to clean up, so that was the issue I assume (you can see, I used the 'weight' attribute, not the '.log.weight', so sad :) ). Now, I can reproduce your plots, and yes, it seems to work well. Thank You!

[kingaa]

As of release 4.4.3.0, this functionality is now part of pomp.

[jeswheel]

Hi @kingaa, lately I have been playing around with the saved latent states obtained from the pfilter function. After reading the documentation for the save.states argument of pfilter and then comparing it to the actual code, I too was a bit confused with what was going on, and found this nice discussion on the topic. While the documentation of save_states does clear some of this up, I'm wondering if it worth revisiting the save.states argument and documentation. I'm happy to make changes and submit a pull-request, but wanted your thoughts first.

• save.states="unweighted": When I first read this, I thought unweighted particles referred to the prediction particles $X_{n, j}^P$ (line 3 of Algorithm 1), as these are effectively unweighted simulations from the process model, and $X_{n, j}^F$ are "re-weighted" simulation values. This is cleared up in the save_states documentation, which states: "These are the unweighted particles, saved after resampling.", but perhaps a better argument name is "resampled" or "filtered" rather than "unweighted"? A similar line of thinking might prefer something like "raw" or "prediction" to get the particles $X_{n, j}^P$ and their weights.
• There's a long vector of possibilities: c("no", "weighted", "unweighted", "FALSE", "TRUE"). Some I think are around just for historical purposes, but since several do the same thing, maybe reducing this to c("no", "weighted", "unweighted") makes sense, or c("FALSE", "weighted", "unweighted") since this may have more backwards compatibility.

[kingaa]

@jeswheel: thanks for weighing in on this. Yes, it should be possible to make this quite clear. Please do take a stab at it and submit a pull request!

You are correct that the FALSE, TRUE options are there for backwards compatibility. In the next six months or so, a breaking version will be released and it would be a good time to remove those.

[jeswheel]

In that case, I'll start working something up that will include a breaking change, and won't expect it to be merged until you're ready to do so.

[kingaa]

If it's not too much trouble, consider introducing a change that would not break existing code but would issue a warning message when the old options are used. This helps users prepare for the eventual breaking change.