eigen decomposition and type="double" when combined with dnorm vs dunif

[fbud...@gmail.com]

Hi all,

I've been trying to work an eigen decomposition into a nimble model, but am running into an issue. When I try to manipulate the extracted eigenvalue by adding or multiplying it by another parameter, I get the error "cannot coerce type 'closure' to vector of type 'double'." 

But it only seems to happen when the other parameter is drawn from a normal distribution, and if I instead draw it from a uniform distribution, it works. I can't understand why those would produce different behavior. Simplified example below.

Thanks,

Franny

pumpCode <- nimbleCode({
  S ~ dunif(0,2)
  V ~ dunif(0,2)
 
  for (t in 1:T){
    D[1,1,t]<-0
    D[1,2,t]<-V
    D[1,3,t]<-V
   
    D[2,1,t]<-S
    D[2,2,t]<-0
    D[2,3,t]<-0
   
    D[3,1,t]<-0
    D[3,2,t]<-S
    D[3,3,t]<-S
   
    lambda_d[t]<-eigen(D[1:3,1:3,t])$values[1]

    C[t]~dpois(phi[t])
    log(phi[t])~dnorm(a+log(lambda_d[t]),sd=0.01)
  }
 
  #a~dunif(-1,1)
  a~dnorm(0,sd=0.0001)
})

S<-0.95
V<-0.25
T<-50
D<-matrix(c(0,V,V,S,0,0,0,S,S),byrow=T,nrow=3)
lambda_d<-eigen(D)$values[1]

C<-rep(0,T)
for (t in 2:T){
  log.lambda_d<-0.1+log(lambda_d)
  C[t]<-rpois(1,exp(log.lambda_d))
}

pumpInits <- list(S = 1, V = 1, a=0)
pumpConsts <- list(T=T)
pumpData <- list(C=C)

out<-nimbleMCMC(code = pumpCode, constants = pumpConsts,
                data = pumpData,inits=pumpInits,
                nchains = 1, niter = 60000,
                monitors = c("a"))

[Daniel Turek]

Franny, this is a new one for sure.  The error is occurring at the stage of MCMC configuration, when the model object is inspected, and sampler choices are made for updating the dimensions of the model.  Apparently, the use of eigen(D)$values inside the model trips up the system for identifying conjugate dimensions of the model, thus giving the incomprehensible error that you're seeing.  The reason you don't see this when using a uniform prior, is the use of a uniform prior distribution precludes the possibility of possible conjugacy, so the conjugacy checking routine exits sooner.  However, the normal prior distribution allows the possibility of a conjugate relationship, so the conjugacy checking routine proceeds further - far enough to get itself into trouble.

This is something we'll have to look into, and fix.  In the meantime, you can use either prior distribution if you use the (slightly longer) workflow shown below.  I've updated the steps to be consistent with your original call to nimbleMCMC().  In the code below, you'll see the addition of "useConjugacy = FALSE", which skips all conjugacy checking entirely, so this problem doesn't arise.

Please let me know if this works for you.

Daniel

.......

.... earlier code here ....

.......

pumpInits <- list(S = 1, V = 1, a=0)
pumpConsts <- list(T=T)
pumpData <- list(C=C)

Rmodel <- nimbleModel(pumpCode, pumpConsts, pumpData, pumpInits)

conf <- configureMCMC(Rmodel, monitors = c("a"), useConjugacy = FALSE)

Rmcmc <- buildMCMC(conf)

Cmodel <- compileNimble(Rmodel)
Cmcmc <- compileNimble(Rmcmc, project = Rmodel)

out <- runMCMC(Cmcmc, 60000)

--
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/b63ad1cb-8ebd-4b77-a728-e1178e38f1e4n%40googlegroups.com.

[fbud...@gmail.com]

That works Daniel, thank you!