PMCMC samples of latent states

167 views
Skip to first unread message

Tom Hull

unread,
Nov 17, 2019, 12:56:43 PM11/17/19
to nimble-users
Hi

Pretty sure I'm missing something obvious, but how do I extract the samples for the latent states when using the `RW_PF_block` MCMC sampler?

e.g.
timeModelCode <- nimbleCode({
  x[1] ~ dnorm(mu_0, 1)
  y[1] ~ dnorm(x[1], 1)
  for(i in 2:t){
    x[i] ~ dnorm(x[i-1] + c, pro_error)
    y[i] ~ dnorm(x[i], obs_error)
  }

  pro_error ~ dnorm(0, 1)
  obs_error ~ dnorm(0, 1)
  c ~ dnorm(1,1)
  mu_0 ~ dnorm(0, 1)
})

# simulate some data
seed = 159
t <- 25
mu_0 <- -10
pro_error <- 0.01
obs_error <- 0.5
x <- rnorm(1 ,mu_0, pro_error)
y <- rnorm(1, x, obs_error)
c <- 1;
for(i in 2:t){
  x[i] <- rnorm(1, x[i-1] + c, pro_error)
  y[i] <- rnorm(1, x[i], obs_error)
}

# build the model
rTimeModel <- nimbleModel(timeModelCode, constants = list(t = t),
                          data <- list(y = y), check = FALSE )

cTimeModel <- compileNimble(rTimeModel)

timeConf <- configureMCMC(rTimeModel)

# Add random walk PMCMC sampler with particle number optimization.
timeConf$addSampler(target = c("obs_error", "pro_error", "c", "mu_0"), type = "RW_PF_block",
                    control <- list(propCov= diag(4), adaptScaleOnly = FALSE,
                                    pfType = "auxiliary", pfControl = list(saveAll = T, smoothing = T),
                                    latents = "x", pfOptimizeNparticles = TRUE))
RtimeConf <- buildMCMC(timeConf)
CtimeConf <- compileNimble(RtimeConf, project=rTimeModel, resetFunctions=T)
CtimeConf$run(500)

CtimeConf$mvSamples # This doesn't contain "x"!

Many Thanks

Daniel Turek

unread,
Nov 17, 2019, 1:19:29 PM11/17/19
to Tom Hull, nimble-users
For recording MCMC samples, you'll need to add a monitor for "x" using the MCMC configuration object:

timeConf$addMonitors("x")

The options you tried, saveAll = TRUE, will control how the nested particle filter algorithm (inside the RW_PF_block MCMC sampler) will operate, but that's located much deeper within the MCMC that you're using.  It sounds like you just want the MCMC algorithm to monitor "x".



--
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 on the web visit https://groups.google.com/d/msgid/nimble-users/7a90d35e-75bf-478d-8f9c-af4fa4132055%40googlegroups.com.

Chris Paciorek

unread,
Nov 21, 2019, 11:21:22 AM11/21/19
to Tom Hull, nimble-users
Hi Tom,

One other note is that in the upcoming release (which should be done
in a couple weeks), we're fixing a bug affecting use of the auxiliary
particle filter with saveAll=FALSE.
If you're sampling the latent nodes in PMCMC using the auxiliary
filter, you should be fine, because in that case we set saveAll=TRUE.

But on the off-chance you didn't want to monitor the latent nodes (in
which case we set saveAll=FALSE), please update
to version 0.9.0 when it comes out (or in the interim, install nimble
off our 'devel' branch on github).

-chris
> To view this discussion on the web visit https://groups.google.com/d/msgid/nimble-users/CAKbe0hr9embgwz9jj1FhB%2BLBhYOGk3o3hGrWTtSvGfyr_E%2B98Q%40mail.gmail.com.

Tom Hull

unread,
Nov 21, 2019, 12:20:54 PM11/21/19
to nimble-users
Hi Chris

Thanks, how does the bug manifest?

I'm having some possibly related problems where what I thought was a subtle change to the model now means the model takes a very very long time to sample, whereas it took a few seconds before.
>> To unsubscribe from this group and stop receiving emails from it, send an email to nimble...@googlegroups.com.
>> To view this discussion on the web visit https://groups.google.com/d/msgid/nimble-users/7a90d35e-75bf-478d-8f9c-af4fa4132055%40googlegroups.com.
>
> --
> 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...@googlegroups.com.

Chris Paciorek

unread,
Nov 23, 2019, 3:32:19 PM11/23/19
to Tom Hull, nimble-users
Hi Tom,

I think what was happening was some incorrect copying of values during the algorithm. 
So numerical results would be incorrect. I haven't dug into the changes in the code, but 
I doubt it would affect sampling time.

Chris

To unsubscribe from this group and stop receiving emails from it, send an email to nimble-users...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/nimble-users/05ef22de-307b-49a6-b8d8-d67ab0ec0213%40googlegroups.com.

Tom Hull

unread,
Nov 25, 2019, 5:30:21 AM11/25/19
to nimble-users
Thanks Chris,

Yes I think it's to do with changing the priors for the initial state.

semi-related question, can the inbuilt auxiliary PF cope with missing data? if so, how would I code that up here?

Thanks

Chris Paciorek

unread,
Nov 26, 2019, 11:49:54 AM11/26/19
to Tom Hull, nimble-users
If you're using the PF as part of particle MCMC then NIMBLE will just sample the missing observations as part of the MCMC. So I think that should work fine.

If you're just applying a PF directly, not as part of an MCMC, that won't work as some likelihood terms would be NA. You might instead be able to set up the model
so there are no nodes for the missing observations, and I suspect that would work, but haven't thought carefully about it. Happy to discuss further if this is of interest.

-chris

--
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 on the web visit https://groups.google.com/d/msgid/nimble-users/213312d8-521f-4a1c-8962-12a2c1373e38%40googlegroups.com.

Tom Hull

unread,
Dec 1, 2019, 5:45:28 PM12/1/19
to nimble-users
Great, thanks Chris.

I've given it a go, but I think I'm missing something again. I can't find any similar examples to check against, but with the below model the MCMC is not doing a good job of estimating the process and observation errors.
The other two parameters are estimated reasonably well.
I'm not sure what I've done wrong that means I can't recover these top level parameters.

```
require(nimble)

# Building a simple linear state-space model.
# x is latent space, y is observed data
timeModelCode <- nimbleCode({
  x[1] <- x_0
  y[1] ~ dnorm(x[1], 0.3)

  for(i in 2:t){
    x[i] ~ dnorm(x[i-1] + c, pro_error)
    y[i] ~ dnorm(x[i], obs_error)
  }

  x_0 ~ dnorm(0, 1)
  obs_error ~ dunif(0, 10)
  pro_error ~ dunif(0, 10)
  c ~ dnorm(1,1)

})

# simulate some data
set.seed(125)
t <- 50
x_0 <- 0.5
pro_error <- 0.1
obs_error <- 0.3
x <- x_0

y <- rnorm(1, x, obs_error)
c <- 0.7;
for(i in 2:t){
  x[i] <- x[i-1] + c

  y[i] <- rnorm(1, x[i], obs_error)
}


# y[10:30] = NA # missing values testing, not used

plot(x, type="l")
points(y, col="red")


# build the model
rTimeModel <- nimbleModel(timeModelCode, constants = list(t = t),
                          data <- list(y = y), check = FALSE )
rTimeModel$setInits(list(c = 1, x_0 = 1))

cTimeModel <- compileNimble(rTimeModel)

timeConf <- configureMCMC(rTimeModel)

# Add random walk PMCMC sampler with particle number optimization.
timeConf$addSampler(target = c("c", "x_0", "obs_error", "pro_error"), type = "RW_PF_block",
                    control <- list(propCov= diag(4), adaptScaleOnly = T,
                                    pfType = "bootstrap", pfControl = list(saveAll = T, smoothing = T),

                                    latents = "x", pfOptimizeNparticles = TRUE))
timeConf$addMonitors("x")

RtimeConf <- buildMCMC(timeConf)
CtimeConf <- compileNimble(RtimeConf, project=rTimeModel, resetFunctions=T)
CtimeConf$run(1000)

hist(as.matrix(CtimeConf$mvSamples, "obs_error"), breaks=100)
abline(v=obs_error, col="red"
hist(as.matrix(CtimeConf$mvSamples, "pro_error"), breaks=100)
abline(v=pro_error, col="red")
hist(as.matrix(CtimeConf$mvSamples, "c"), breaks=100)
abline(v=c, col="red")
hist(as.matrix(CtimeConf$mvSamples, "x_0"), breaks=100)
abline(v=x_0, col="red")
```
To unsubscribe from this group and stop receiving emails from it, send an email to nimble...@googlegroups.com.

Chris Paciorek

unread,
Dec 2, 2019, 10:26:24 AM12/2/19
to Tom Hull, nimble-users
hi Tom, do you mean to specify the dnorm using the precision rather than standard deviation or variance?

-Chris

To unsubscribe from this group and stop receiving emails from it, send an email to nimble-users...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/nimble-users/8ff30da2-a016-427d-8581-83630d0c7ad0%40googlegroups.com.

Tom Hull

unread,
Dec 2, 2019, 7:25:09 PM12/2/19
to nimble-users
haha, yes, that would be it. Thank you.

Tom Hull

unread,
Dec 4, 2019, 6:01:52 PM12/4/19
to nimble-users
Is it possible to have more than one latent state space variable?

I get the error "all latent nodes must come from same variable" when trying to specify a model that as a time varying "c".

Perry de Valpine

unread,
Dec 4, 2019, 7:21:30 PM12/4/19
to Tom Hull, nimble-users
It is possible, but they must be set up in the same variable.

For example, x[t, 1:nStates] can be a vector of the states at time t, and x[t+1, 1:nStates] depends on x[t, 1:nStates] however you need.  You can still declare x[t, 1], x[t, 2], etc. separately, so it should be simply a matter of how latent states are named and organized.



To unsubscribe from this group and stop receiving emails from it, send an email to nimble-users...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/nimble-users/fdc25f83-7f74-43b8-a48d-29652769fffa%40googlegroups.com.

Tom Hull

unread,
Jan 15, 2020, 9:24:03 AM1/15/20
to nimble-users
Hi again

I'm having some trouble understanding how to scale my state space model for variable time intervals correctly.
basically the below model works, but only if dT = 1, if I change it to something closer to my real problem say 600. many of the samples for x[,1] end up being very very small, with low logProb of data node errors. Not quite sure what is wrong.

Example code:
require(nimble)
require(ggplot2)
require(data.table)

timeModelCode <- nimbleCode({
    # inital states
  x[1,1] ~ dnorm(y[1], sd = obs_error) # C
  x[1,2] ~ dnorm(0, sd = 5) # J
  x[1,3] ~ dnorm(kw[1] * (Csat - x[1,1]), sd = 0.0001) # gas exchange

    # state space evolution
  for(i in 2:t){
    x[i,1] ~ dnorm(x[i-1,1] + (x[i-1,2] * dT) + (x[i-1,3] * dT), sd = pro_error) # change in C
    x[i,2] ~ dnorm(x[i-1,2], sd = var_J) # J, (random walk)
    x[i,3] ~ dnorm(kw[i] * (Csat - x[i,1]), sd = 0.0001) # G
    y[i] ~ dnorm(x[i,1], sd = obs_error) # observation
  }

  obs_error ~ dnorm(0, sd=5)
  pro_error ~ dnorm(0, sd=1)
  var_J ~ dnorm(0, 1) # how much J varies between steps

})

# simulate some data

t <- 48
pro_error <- 0.1
obs_error <- 0.5
x <- rep(300, t) # inital state + populate vector
y <- rnorm(1, x, obs_error) # first observation
dT = 600 # time between observations (seconds)

J = (8 * sin(0.03 * pi * 1:t) + 0.5) / dT # our unobserved latent state of interest
kw = rep(3e-05, t)
Csat = 250
G = rep(0, t)
G[1] <- kw[1] * (Csat - x[1])

for(i in 2:t){
  x[i] <- x[i-1] + (J[i-1] * dT) + (G[i-1] * dT)
  G[i] <- rnorm(1, kw[i] * (Csat - x[i]), 0.0001)

  y[i] <- rnorm(1, x[i], obs_error)
}

truth = data.table(x, J, G, t = 1:t)
truth = melt(truth, id.vars = "t")
truth[, var := as.numeric(as.factor(variable))]


# y[10:20] = NA # for testing missing observations


plot(x, type="l")
points(y, col="red")

plot(J, type="l")
lines(G, col="blue")
plot(J + G, type="l", col="red")

# build the model
rTimeModel <- nimbleModel(timeModelCode, constants = list(t = t, dT = dT),
                          data <- list(y = y, kw = kw, Csat = Csat), check = T,
                          inits = list(var_J = 1, obs_error = 0.1, pro_error = 0.1)) # should x be init?
rTimeModel$initializeInfo()

cTimeModel <- compileNimble(rTimeModel)

timeConf <- configureMCMC(rTimeModel)

# Add random walk PMCMC sampler with particle number optimization.

timeConf$addSampler(target = c("var_J", "obs_error", "pro_error"), type = "RW_PF_block",
                    control <- list(pfType = "auxiliary", pfControl = list(saveAll = T, smoothing = T),
                                    latents = c("x")))

timeConf$addMonitors("x")
RtimeConf <- buildMCMC(timeConf)
CtimeConf <- compileNimble(RtimeConf, project=rTimeModel, resetFunctions=T)
CtimeConf$run(5000)

  # check top level parameters against truth
hist(as.matrix(CtimeConf$mvSamples, "obs_error"), breaks=100)
abline(v=obs_error, col="red")

hist(as.matrix(CtimeConf$mvSamples, "pro_error"), breaks=100)
abline(v=pro_error, col="red")

hist(as.matrix(CtimeConf$mvSamples, "var_J"), breaks=100)
abline(v=sd(diff(J)), col="red")


  # extract time series latent state samples, if there is more elegant way of doing this I don't know it.
out = as.matrix(CtimeConf$mvSamples, "x")
out = melt(data.table(out))
out = out[,.(mean = mean(value, na.rm = F), min= quantile(value, 0.05, na.rm=F), max=quantile(value, 0.95, na.rm=F)), by=variable]
out[, t := stringr::str_extract_all(variable, "\\d+", simplify = T)[,1]]
out[, var := stringr::str_extract_all(variable, "\\d+", simplify = T)[,2]]
out[, t := as.numeric(t)]
#
ggplot(out) +
  geom_ribbon(aes(t, ymin=min, ymax=max), alpha=0.5) +
  geom_point(data=truth, aes(t, value)) +
  facet_grid(var ~ ., scales = "free_y") +
  geom_line(aes(t, mean))



Perry de Valpine

unread,
Jan 15, 2020, 8:37:33 PM1/15/20
to Tom Hull, nimble-users
Hi Tom,

It's not clear if this is a nimble issue or an issue with the statistical modeling problem.  It sounds like the latter, unless you are sure about the results you should be getting and believe the PMCMC in nimble is not working correctly.  Is that the case?

-Perry


To unsubscribe from this group and stop receiving emails from it, send an email to nimble-users...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/nimble-users/8eba743f-4914-4426-b4fc-40ad14474354%40googlegroups.com.

Tom Hull

unread,
Jan 16, 2020, 9:36:59 AM1/16/20
to nimble-users
It's possibly a modelling problem, but if so I don't understand the nimble error messages.

Looking closely the samples for x[40,1] for instance are oddly skewed. most are normally distributed around the true value, but some (the first samples) are very very high (should be ~ 400, but closer to 70,000).

Not enough burn-in for the sampler?
Bad initiations for the x variable?

Perry de Valpine

unread,
Jan 16, 2020, 12:20:30 PM1/16/20
to Tom Hull, nimble-users
I haven't tried to run your code yet.  Could you say what the error messages are?  I'm not seeing them in your most recent emails.

The x values should be drawn by the particle filter, starting from the prior for x[1, 1:3].  So I wonder if there are bad initial values for the top-level parameters (obs_error, pro_error, var_J)?  It might help to set the prior for var_J to be only positive.

BTW, nimble 0.9-0 is released and has the updates to particle filtering methods mentioned by Chris earlier in this email thread.

-Perry


To unsubscribe from this group and stop receiving emails from it, send an email to nimble-users...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/nimble-users/40aa090b-44fc-4c20-ba86-491e013e8bba%40googlegroups.com.

Tom Hull

unread,
Jan 17, 2020, 5:24:27 AM1/17/20
to nimble-users
Hi Perry,

As the run starts I get:
warning: logProb of data node y[36]: logProb less than -1e12.
and the first few samples for the latent states at these times (t=36) are wildly wrong.

I'm using nimble 0.9

setting the prior on var_J to dunif(0, 2) didn't help. and there are no samples with var_J negative.
Reply all
Reply to author
Forward
0 new messages