Nimble/JAGS; Nimble can't recover parameters

[fbud...@gmail.com]

Hi all,

We have recently switched some models from JAGS to Nimble, mostly with success. However we have run into an issue with a change-point model, where JAGS can recover the simulated change-point and the means of the two distributions but Nimble can't.

There is also an autocorrelation process, but in the simplified examples below we don't estimate those parameters. We have also tried a block update for mu, but that didn't change the results. Below is our set of code for both JAGS (first) and Nimble (second).

We can't think of a reason why the two behave differently; any suggestions? 

Thank you,

Franny and Katie

## libraries

rm(list=ls())

library(mvtnorm)

library(jagsUI)

library(nimble)

##

## Simulated Data

##

set.seed(23)

T <- 116

mu<-matrix(c(5,10,8,12),byrow=T,nrow=2) ## attractant center

Delta<-matrix(c(0.5,0,0,0.5),nrow=2,byrow=2) ## variance

tau<-52 ## change-point

rho<-(.8)

M<-rho*diag(2)

z<-rep(0,T)

z.idx <- z[1]+1

y<-matrix(0,nrow=T,ncol=2)

y[1,]<-rmvnorm(1,mu[z.idx,],Delta)

for (i in 2:T){

  z[i]<-ifelse(i<tau,0,1)

  z.idx<-z[i]+1

  psi<-t(M%*%y[i-1,])+t((diag(2)-M)%*%(mu[z.idx,]))

  y[i,]<-rmvnorm(1,psi,Delta)

}

plot(y,col=z+1)

##

## JAGS

##

j.Model <- tempfile()

writeLines("

model{

  # Priors for mean (attractant)

  mu[1,1:2]~dmnorm.vcov(c(0,0),Sigma)

  mu[2,1:2]~dmnorm.vcov(c(0,0),Sigma)

  

  ## priors for variance in likelihood

  zeta[1,2]~dgamma(0.001,0.001)

  zeta[1,1]~dgamma(0.001,0.001)

  pi[1,1]<-1/zeta[1,1]

  pi[1,2]<-1/zeta[1,2]

  

  Delta[1,1]<-pi[1,1] 

  Delta[2,2]<-pi[1,2]

  Delta[2,1]<-0

  Delta[1,2]<-Delta[2,1]

  

  prob<-rep(1,T)

  tau~dcat(prob)

  z[1]<-1   

  for (i in 2:T){

    z[i]<-ifelse(i<tau, 1, 2) 

  }

  

  #Likelihood  

  locations[1,]~dmnorm.vcov(mu[z[1],],Delta)

  for (i in 2:T){

    locations[i,]~dmnorm.vcov(locations[i-1,]%*%M+mu[z[i],]%*%I.M,Delta)

  }   

}

",con=j.Model)

j.Data<-list(locations=y,T=T,M = diag(0.8,2,2),I.M = diag(0.2,2,2),Sigma = matrix(c(10,0,0,10),nrow=2,byrow=2))

j.Params<-c("mu","tau","Delta")

j.Inits <- function() list(mu = matrix(c(5,5,5,5),byrow=T,nrow=2),zeta = matrix(rep(5,2),ncol=2),tau = 5)

j.Results<-jags(data=j.Data, inits=j.Inits, parameters.to.save = j.Params, model.file = j.Model, n.chains=1, n.adapt=1000, n.iter = 5000, n.burnin = 2500)

plot(j.Results)

colMeans(j.Results$samples[[1]])

##

## NIMBLE

##

n.Code <- nimbleCode( {

  # Priors for mean (attractant)

  chi[1:2]<-c(0,0)

  mu[1,1:2]~dmnorm(chi[1:2],cov=Sigma[1:2,1:2])

  mu[2,1:2]~dmnorm(chi[1:2],cov=Sigma[1:2,1:2])

  

  Delta[1,1]~dinvgamma(1,1)

  Delta[2,2]~dinvgamma(1,1)

  Delta[2,1]<-0

  Delta[1,2]<-Delta[2,1]

  

  ## priors for tau

  prob[1:T]<-rep(1,T)

  tau[1,1]~dcat(prob[1:T])

  ## z indexing

  z[1,1]<-1

  for(j in 2:T){

    z[j,1] <- (j < tau[1,1]) + 1

  }

  ##likelihood

  y[1,1:2]~dmnorm(mu[z[1,1],1:2],cov=Delta[1:2,1:2])

  psi[1,1:2] <- mu[z[1,1],1:2]

  for(i in 2:T){

    psi[i,1:2]<-t(M[1:2,1:2]%*%y[i-1,1:2])+t(I.M[1:2,1:2]%*%mu[z[i,1],1:2])

    y[i,1:2]~dmnorm(psi[i,1:2],cov=Delta[1:2,1:2])

  }

})

n.Data <- list( y = y)

n.Const <- list(T = T, Sigma = diag(10,2,2), M = diag(0.8,2,2), I.M = diag(0.2,2,2))

n.Inits <- list(mu = matrix(c(5,5,5,5),byrow=T,nrow=2), Delta = diag(5,2,2),tau = 5)

n.Model <- nimbleModel(code=n.Code, data = n.Data, constants = n.Const, inits = n.Inits)

conf.model <- configureMCMC(n.Model)

CmyModel <- compileNimble(n.Model) 

myMCMC <- buildMCMC(conf.model)

CmyMCMC <- compileNimble(myMCMC,project=n.Model,resetFunctions = TRUE)

n.Results<-runMCMC(CmyMCMC,niter=5000,nburnin =2000)

par(mfrow = c(3, 2))

colMeans(n.Results)

plot(n.Results[,5],typ="l")

plot(n.Results[,6],typ="l")

plot(n.Results[,7],typ="l")

plot(n.Results[,8],typ="l")

plot(n.Results[,9],typ="l")

[Perry de Valpine]

Hi Franny and Katie,

This was a puzzle.  The problem seems to be in the translation of the ifelse in JAGS to a logical in NIMBLE.  In the JAGS code, time steps before the change point (tau) get index 1.  In the NIMBLE code, time steps before the change point get index 2 (from TRUE + 1).  (This meant that for tau of 1 or 2, two of the mu elements were unused and followed their priors, and related issues.) When I replaced the relevant line with

z[j,1] <- (j >= tau[1,1]) + 1

it appeared to work and at a quick glance to match the posterior from JAGS.

Some other notes: I believe the ifelse should have worked in NIMBLE but I see it gives an error, which perhaps you hit as well.  I will file an issue for us to look at that.  I think the dinvgamma equivalent to what you have in JAGS would be dinvgamma(0.001, 0.001).  You may have tweaked that in exploring the puzzle of what was happening.  I replaced that section with exactly the way you had it in the JAGS code as I checked on what was happening.  An alternative way to write the linear algebra piece is:

    psi[i,1:2]<- (y[i-1,1:2]%*%M[1:2,1:2]+mu[z[i,1],1:2]%*%I.M[1:2,1:2])[1,1:2]

At some point it looks like you put in an extra index in the NIMBLE code for z, and I'm not sure if that was from some type issue.

But the main thing looks to me like accidentally flipped logic on putting 1s vs 2s into elements of z.  Does that look correct and resolve the issue on your end?

Best,

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 on the web visit https://groups.google.com/d/msgid/nimble-users/8201e188-4137-4fe2-b3bc-10f372a45e82n%40googlegroups.com.

[fbud...@gmail.com]

That did resolve our issues, thank you!