using particle filters

[Brook Milligan]

I’m trying to learn nimble particle filters for a time series and running into a problem. In this little example, the top level parameters are the initial condition, n0, the growth rate, lambda, and the standard deviation, sigma; latent variables are x[i], and data are y[i]. Based upon the nimbleSMC paper, my model is

MCMC_code <- nimbleCode({
n0 ~ dnorm(10,50)
lambda ~ dunif(0.9,1.1)
sigma ~ dunif(0,10)
x[1] ~ dnorm(mean=n0, sd=sigma)
for (i in 2:n) {
x[i] ~ dnorm(mean=lambda * x[i-1], sd=sigma)
}
for (i in 1:n) {
y[i] ~ dpois(lambda=x[i])
}
})

The samplers are defined with

auxpf <- buildAuxiliaryFilter(nimble_model, “x",
control = list(saveAll = FALSE, smoothing = FALSE, initModel = FALSE))
propSD <- c(0.089, 0.089, 0.039)
nimble_configured$addSampler(target = c("n0", "lambda", "sigma"),
type = "RW_PF_block", control = list(propCov = diag(propSD^2),
pf = auxpf, adaptive = TRUE, pfNparticles = MCMC_particles, latents = “x"))

The results are nonsensical (to me). In short, n0 goes straight to the minimum value and the trace plots flatline.

If I remove n0 from the likelihood altogether, and give x[1] a prior like dnorm(30,50), the model behaves much better, but still does not fit the data well.

I expect this is all me not understanding how to use the particle filters, especially to get the initial condition for a hidden Markov chain.

Any help is greatly appreciated.

Thanks a lot.

Cheers,
Brook

[Perry de Valpine]

Hi Brook,

Can you say more about the data (y) and whether you are experimenting with simulated data to get started?

I notice that your setup is somewhat different from some common time-series models. Typically one sees the auto-correlation parameter (lambda here, often called rho) constrained between -1 and 1, sometimes with a separate time trend if that is relevant. And then it is common to include an intercept (to lambda * x[i-1]).

Often with Poisson data one sees a log link (e.g. dpois(lambda=exp(x[i])), which avoids negative values and often works out reasonably. Maybe since you are starting at a positive n0 and limiting lambda between 0.9 and 1.1 you are not worried about that.

It could be useful for exploratory purposes to run some PF likelihoods for specific parameters, to understand if it really is the case that the minimum n0 gives higher likelihoods.

HTH

Perry

--
You received this message because you are subscribed to the Google Groups "nimble-users" group.
To unsubscribe from this group and stop receiving emails from it, send an email to nimble-users...@googlegroups.com.
To view this discussion visit https://groups.google.com/d/msgid/nimble-users/6C998BDB-717B-4183-9E97-7717A61F9018%40biology.nmsu.edu.

[Brook Milligan]

Hi Perry,

Thanks for responding.

I’m just trying to explore this with simple data. I’ve attached the full example file so you can see the details. The data are just some numbers that are declining pretty evenly; I’m mainly trying to determine the rate of decline and the initial value. I’m probably doing something wrong, but the results do not look plausible. Perhaps you will notice something that will help me better understand this.

As for your specific questions: in this case, the rate of decline is pretty close to one and I’m comfortable (for now) not worrying about anything that makes this less simple. The deviation around x[i] is also small enough that negative values should not occur; again, just trying to keep things simple. Given that I’m just starting with this, I don't understand how to get PF likelihoods.

Thanks a lot for your help.

Cheers,
Brook

> To view this discussion visit https://groups.google.com/d/msgid/nimble-users/CABeTmqhXmpO1UqtWoyaG5RC6K8GRONfvUVKPUirAbRig6wz0Wg%40mail.gmail.com.

[Perry de Valpine]

Hi Brook,

I explored your example for a while. I ran regular MCMC on the whole thing and also saw x0 centered around 10.

The problem seems to be that you might have meant "sd = 50" for the prior on x0, but you just had "50". The default second parameter of dnorm is precision, so you have a prior with mean 10 and sd = 1/sqrt(50) = 0.14. So the prior was strongly holding x0 near 10. When I change that to "sd = 50", results are more reasonable, although you'll want to experiment with the proposal std. dev. values.

HTH

Perry