question about the functioning of buildLaplace

[frederic....@gmail.com]

Dear all,

in an effort to check what is included in the marginalized log-Likelihood of a model buildt using builLaplace, I have changed one of the examples given with buildLaplace to have the logLikelihood constant. I expected the output not to depend on model parameters. Yet it does. Here is the code:

library(nimble)

pumpCode <- nimbleCode({
  for (i in 1:N){
    theta[i] ~ dgamma(exp(alpha), exp(beta))

    lambda[i] <- theta[i] * t[i]

    #x[i] ~ dpois(4.0) # does not work

x[i] ~ dpois(0*lambda[i]+4.0)

  }

  #alpha ~ dexp(1.0)

  #beta ~ dgamma(0.1, 1.0)

alpha ~ dnorm(0.0,1.0)

  beta ~ dnorm(0.0, 1.0)

 

})

 

pumpConsts <- list(N = 10, t = c(94.3, 15.7, 62.9, 126, 5.24, 31.4, 1.05, 1.05, 2.1, 10.5))

 

pumpData <- list(x = c(5, 1, 5, 14, 3, 19, 1, 1, 4, 22))

 

pumpInits <- list(alpha = 0.1, beta = 0.1, theta = rep(0.1, pumpConsts$N))

 

pump <- nimbleModel(code = pumpCode, name = "pump", constants = pumpConsts,
                    data = pumpData, inits = pumpInits, buildDerivs = TRUE)

 

pumpLaplace <- buildLaplace(pump,c('alpha','beta'))

 

CpumpLaplace <- compileNimble(pumpLaplace, project = pump)

 

CpumpLaplace$calcLaplace(c(1,1))

CpumpLaplace$calcLaplace(c(0,0))

CpumpLaplace$calcLaplace(c(0,1))

CpumpLaplace$calcLaplace(c(1,0))

 

Results equal -64.38562 with c(0,0) and c(0,1) and -63.88024 with c(1,1) and c(1,0). Do these results make sense ? I am aware that I probably push the system at its margins since the logLikelhood no longer depends on the parameters...

Thank in advance for enlightening me on what happens here.

All the best,

Frédéric

[Perry de Valpine]

Dear Frédéric,

Thanks for the question. It's interesting to see your exploration.

Here is what I believe is happening. buildLaplace, by default, uses the model graph (what depends on what) to decide what should be treated as a random effect to integrate over. When you use x[i] ~ dpois(0*lambda[i]+4.0)there is a graph dependency that x[i] depends on lambda[i] (even though we humans can see it is always multiplied by 0), so theta[i] (on which lambda[i] depends deterministically) can be determined to be a random effect. When you use x[i] ~ dpois(4.0)nimble believes (correctly) that there is no dependency from theta-->lambda-->x, and so theta is not a random effect.

You can see the relevant sets of (default) nodes by

setupMargNodes(pump, c('alpha','beta'))

There are other arguments to setupMargNodes and/or buildLaplace for fine-grained control, all the way up to specifying your choices of parameters, random effects, and calculation nodes.

There is a first attempt to explain this and guide you to setupMargNodes in help(buildLaplace), but I'm sure it could be improved and we'd welcome your suggestions.

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 on the web visit https://groups.google.com/d/msgid/nimble-users/6cd3df58-7898-4d4e-bee5-aa74ef6aba80n%40googlegroups.com.

[frederic....@gmail.com]

Thank you Perry for the explanation that is very clear and I think this is what I had in mind (and why I introduced the 0* in the code).

Yet my question was also about the marginal logLik values themselves that are depend strangely from the hyperparameter values. I would have expected them to be the same: the integration over the random effects should be 0 in logLik (since they are deconnected from observations), whatever the parameter values: but maybe I'm wrong on this.

Sincerely,

Frédéric

[Wei Zhang]

Hi Frédéric,

For Laplace approximation, we need to maximize wrt random effects the complete data likelihood p(x, theta | alpha, beta) = p(theta | alpha, beta) p(x | theta, alpha, beta) for any pair of values for alpha and beta. In your case, p(x | theta, alpha, beta) is always a constant given the setup, but p(theta | alpha, beta) depends on alpha and beta, so it makes sense that the result depends on the value of alpha and beta. I guess the four values are not equal to each other at all but some of them are close so after rounding they appear to be equal. Hope this helps. 

Best wishes,

Wei

To view this discussion on the web visit https://groups.google.com/d/msgid/nimble-users/6a40a57e-e9fe-4bf2-a56c-3e84623aca7cn%40googlegroups.com.

[Wei Zhang]

PS: for a quick check you can run the code below, where the constant log-likelihood component from x is separated out  

pumpLaplace <- buildLaplace(pump, paramNodes = c('alpha','beta'), randomEffectsNodes = "theta", calcNodes = "theta",
                            control = list(check = FALSE))

CpumpLaplace <- compileNimble(pumpLaplace, project = pump)

CpumpLaplace$calcLaplace(c(1,1))
CpumpLaplace$calcLaplace(c(0,0))
CpumpLaplace$calcLaplace(c(0,1))
CpumpLaplace$calcLaplace(c(1,0))

CpumpLaplace$calcLaplace(c(1,1)) + pump$calculate("x")
CpumpLaplace$calcLaplace(c(0,0)) + pump$calculate("x")
CpumpLaplace$calcLaplace(c(0,1)) + pump$calculate("x")
CpumpLaplace$calcLaplace(c(1,0)) + pump$calculate("x")

[frederic....@gmail.com]

Thank you Wei for the answers. I tried your code - using project= compileNimble(pump)  and found what you told. Yet I have difficulties understanding this: for me, as it integrates over theta it should be close to 0.

Actually, if I change the distribution of theta to a Gaussian one this is what I obtain - approximately: (results all less than 1e-11) (see code below). Would there be problems of Laplace marginalization with some distributions, esp. asymetric ones?

library(nimble)

pumpCode <- nimbleCode({
  for (i in 1:N){

    theta[i] ~ dnorm(alpha,exp(beta))

    lambda[i] <- theta[i] * t[i]

    #x[i] ~ dpois(4.0) # does not work

x[i] ~ dpois(0*lambda[i]+4.0)

  }

  #alpha ~ dexp(1.0)

  #beta ~ dgamma(0.1, 1.0)

alpha ~ dnorm(0.0,1.0)

  beta ~ dnorm(0.0, 1.0)

 

})

 

pumpConsts <- list(N = 10, t = c(94.3, 15.7, 62.9, 126, 5.24, 31.4, 1.05, 1.05, 2.1, 10.5))

 

pumpData <- list(x = c(5, 1, 5, 14, 3, 19, 1, 1, 4, 22))

 

pumpInits <- list(alpha = 0.1, beta = 0.1, theta = rep(0.1, pumpConsts$N))

 

pump <- nimbleModel(code = pumpCode, name = "pump", constants = pumpConsts,
                    data = pumpData, inits = pumpInits, buildDerivs = TRUE)

pumpLaplace <- buildLaplace(pump, paramNodes = c('alpha','beta'), randomEffectsNodes = "theta", calcNodes = "theta",

                            control = list(check = FALSE))

 

 

CpumpLaplace <- compileNimble(pumpLaplace, project = compileNimble(pump))

 

CpumpLaplace$calcLaplace(c(1,1))

CpumpLaplace$calcLaplace(c(0,0))

CpumpLaplace$calcLaplace(c(0,1))

CpumpLaplace$calcLaplace(c(1,0))

[frederic....@gmail.com]

Dear all,

I am reopening this discussion as I actually meet a series of problems using buildLalpace. My (three) problems are the following ones - with the version 1.0.1 of Nimble in Windows and 1.1.0 on Linux:

     (i) when part of the data Likelihood is not related to the random effects being integrated out, they seem not to be included in the calculation of the likelihood in $calcLaplace: here are two codes that should give the same results but that do not, due to this "problem". I am unsure this is a problem. At least it should be mentioned in the Nimble help what is included in calcLaplace. Note that I have to run these two codes in different R sessions due to the problem reported in point (ii) below. :

 

library(nimble)

pumpCode <- nimbleCode({
  for (i in 1:N){

    theta[i] ~ dnorm(alpha,exp(beta))

    lambda[i] <- theta[i] * t[i]

    #x[i] ~ dpois(4.0) # does not work

x[i] ~ dpois(0*lambda[i]+4.0)

x2[i] ~ dpois(4.0)

 

  }

  #alpha ~ dexp(1.0)

  #beta ~ dgamma(0.1, 1.0)

alpha ~ dnorm(0.0,1.0)

  beta ~ dnorm(0.0, 1.0)

 

})

 

pumpConsts <- list(N = 10, t = c(94.3, 15.7, 62.9, 126, 5.24, 31.4, 1.05, 1.05, 2.1, 10.5))

 

pumpData <- list(x = c(5, 1, 5, 14, 3, 19, 1, 1, 4, 22), x2 = c(5, 1, 5, 14, 3, 19, 1, 1, 4, 22))

 

pumpInits <- list(alpha = 0.1, beta = 0.1, theta = rep(0.1, pumpConsts$N))

 

pump <- nimbleModel(code = pumpCode, name = "pump", constants = pumpConsts,
                    data = pumpData, inits = pumpInits, buildDerivs = TRUE)

 

pumpLaplace <- buildLaplace(pump, randomEffectsNodes = "theta")

 

 

CpumpLaplace <- compileNimble(pumpLaplace, project = compileNimble(pump))

 

CpumpLaplace$calcLaplace(c(1,1))

# -63.57501 (not OK)

 

 

 

 

 

library(nimble)

pumpCode <- nimbleCode({
  for (i in 1:N){

    theta[i] ~ dnorm(alpha,exp(beta))

    lambda[i] <- theta[i] * t[i]

    #x[i] ~ dpois(4.0) # does not work

x[i] ~ dpois(0*lambda[i]+4.0)

x2[i] ~ dpois(0*lambda[i]+4.0)

 

  }

  #alpha ~ dexp(1.0)

  #beta ~ dgamma(0.1, 1.0)

alpha ~ dnorm(0.0,1.0)

  beta ~ dnorm(0.0, 1.0)

 

})

 

pumpConsts <- list(N = 10, t = c(94.3, 15.7, 62.9, 126, 5.24, 31.4, 1.05, 1.05, 2.1, 10.5))

 

pumpData <- list(x = c(5, 1, 5, 14, 3, 19, 1, 1, 4, 22), x2 = c(5, 1, 5, 14, 3, 19, 1, 1, 4, 22))

 

pumpInits <- list(alpha = 0.1, beta = 0.1, theta = rep(0.1, pumpConsts$N))

 

pump <- nimbleModel(code = pumpCode, name = "pump", constants = pumpConsts,
                    data = pumpData, inits = pumpInits, buildDerivs = TRUE)

 

pumpLaplace <- buildLaplace(pump, randomEffectsNodes = "theta")

 

 

CpumpLaplace <- compileNimble(pumpLaplace, project = compileNimble(pump))

 

CpumpLaplace$calcLaplace(c(1,1))

#-127.15 (OK)

 

(ii) my second problem - which either is specific to Windows are to version 1.0.1 of Nimble - is that when I modify the modelCode in the same R session and run different buildLaplace calls I get either problems with the compileNimble phase or R crashes. For example with the code above ran in only one R session, R crashes (maybe not systematically). If I do the above code without the last $calcLaplace, R does not crash but then the CompileNimble call returns these following strange messages:

Compiling

  [Note] This may take a minute.
  [Note] Use 'showCompilerOutput = TRUE' to see C++ compilation details.
Error, could not find name
Name = lifted_d0_times_lambda_oBi_cB_plus_4_L5
Available Name 1 = alpha
Error, could not find name
Name = lifted_d0_times_lambda_oBi_cB_plus_4_L5
Available Name 1 = alpha
Compiling
  [Note] This may take a minute.
  [Note] Use 'showCompilerOutput = TRUE' to see C++ compilation details.
Error, could not find name
Name = lifted_d0_times_lambda_oBi_cB_plus_4_L5
Available Name 1 = alpha
Error, could not find name
Name = lifted_d0_times_lambda_oBi_cB_plus_4_L5
Available Name 1 = alpha
Error, could not find name
Name = lifted_d0_times_lambda_oBi_cB_plus_4_L5
Available Name 1 = alpha
Error, could not find name
Name = lifted_d0_times_lambda_oBi_cB_plus_4_L5
Available Name 1 = alpha
Error, could not find name
Name = lifted_d0_times_lambda_oBi_cB_plus_4_L5
Available Name 1 = alpha
Error, could not find name
Name = lifted_d0_times_lambda_oBi_cB_plus_4_L5
Available Name 1 = alpha
Error, could not find name
Name = lifted_d0_times_lambda_oBi_cB_plus_4_L5
Available Name 1 = alpha
Error, could not find name
Name = lifted_d0_times_lambda_oBi_cB_plus_4_L5
Available Name 1 = alpha
Error, could not find name
Name = lifted_d0_times_lambda_oBi_cB_plus_4_L5
Available Name 1 = alpha
Error, could not find name
Name = lifted_d0_times_lambda_oBi_cB_plus_4_L5
Available Name 1 = alpha
Error, could not find name
Name = lifted_d0_times_lambda_oBi_cB_plus_4_L5
Available Name 1 = alpha
Error, could not find name
Name = lifted_d0_times_lambda_oBi_cB_plus_4_L5
Available Name 1 = alpha
Error, could not find name
Name = lifted_d0_times_lambda_oBi_cB_plus_4_L5
Available Name 1 = alpha
Error, could not find name
Name = lifted_d0_times_lambda_oBi_cB_plus_4_L5
Available Name 1 = alpha
Error, could not find name
Name = lifted_d0_times_lambda_oBi_cB_plus_4_L5
Available Name 1 = alpha
Error, could not find name
Name = lifted_d0_times_lambda_oBi_cB_plus_4_L5
Available Name 1 = alpha
Error, could not find name
Name = lifted_d0_times_lambda_oBi_cB_plus_4_L5
Available Name 1 = alpha
Error, could not find name
Name = lifted_d0_times_lambda_oBi_cB_plus_4_L5
Available Name 1 = alpha
Error, could not find name
Name = lifted_d0_times_lambda_oBi_cB_plus_4_L5
Available Name 1 = alpha
Error, could not find name
Name = lifted_d0_times_lambda_oBi_cB_plus_4_L5
Available Name 1 = alpha
>  

    (iii) my third problem is the following one: using variations around the above models (parameters of x and/or x2), depending on whether x and x2 share parameters, have distinct parameters... I have found very different behaviors of builLaplace. Here are four examples ran in four different R sessions due to problem (ii) above:

 

library(nimble)

pumpCode <- nimbleCode({
  for (i in 1:N){

    theta[i] ~ dnorm(alpha,exp(beta))

    lambda[i] <- theta[i] * t[i]*exp(meanplot[1])

    #x[i] ~ dpois(4.0) # does not work

x[i] ~ dpois(0*lambda[i]+4.0)

x2[i] ~ dpois(4.0*exp(gamma))

 

  }

  #alpha ~ dexp(1.0)

  #beta ~ dgamma(0.1, 1.0)

alpha ~ dnorm(0.0,1.0)

  beta ~ dnorm(0.0, 1.0)

gamma ~ dnorm(0.0, 1.0)

 

for (i in 1:N)

{meanplot[i]~dnorm(0,1)

}

 

})

 

pumpConsts <- list(N = 10, t = c(94.3, 15.7, 62.9, 126, 5.24, 31.4, 1.05, 1.05, 2.1, 10.5))

 

pumpData <- list(x = c(5, 1, 5, 14, 3, 19, 1, 1, 4, 22), x2 = c(5, 1, 5, 14, 3, 19, 1, 1, 4, 22))

 

pumpInits <- list(alpha = 0.1, beta = 0.1, theta = rep(0.1, pumpConsts$N))

 

pump <- nimbleModel(code = pumpCode, name = "pump", constants = pumpConsts,
                    data = pumpData, inits = pumpInits, buildDerivs = TRUE)

 

pumpLaplace <- buildLaplace(pump, randomEffectsNodes = "theta")

 

 

CpumpLaplace <- compileNimble(pumpLaplace, project = compileNimble(pump))

 

CpumpLaplace$getNodeNamesVec(returnParams = TRUE)

# [1] "alpha"       "beta"        "meanplot[1]" (OK)

 

CpumpLaplace$calcLaplace(c(1,1,0))

#-63.57501

 

 

library(nimble)

pumpCode <- nimbleCode({
  for (i in 1:N){

    theta[i] ~ dnorm(alpha,exp(beta))

    lambda[i] <- theta[i] * t[i]*exp(meanplot[1])

    #x[i] ~ dpois(4.0) # does not work

x[i] ~ dpois(0*lambda[i]+4.0)

x2[i] ~ dpois(4.0*exp(gamma+meanplot[i]))

 

  }

  #alpha ~ dexp(1.0)

  #beta ~ dgamma(0.1, 1.0)

alpha ~ dnorm(0.0,1.0)

  beta ~ dnorm(0.0, 1.0)

gamma ~ dnorm(0.0, 1.0)

 

for (i in 1:N)

{meanplot[i]~dnorm(0,1)

}

 

})

 

pumpConsts <- list(N = 10, t = c(94.3, 15.7, 62.9, 126, 5.24, 31.4, 1.05, 1.05, 2.1, 10.5))

 

pumpData <- list(x = c(5, 1, 5, 14, 3, 19, 1, 1, 4, 22), x2 = c(5, 1, 5, 14, 3, 19, 1, 1, 4, 22))

 

pumpInits <- list(alpha = 0.1, beta = 0.1, theta = rep(0.1, pumpConsts$N))

 

pump <- nimbleModel(code = pumpCode, name = "pump", constants = pumpConsts,
                    data = pumpData, inits = pumpInits, buildDerivs = TRUE)

 

pumpLaplace <- buildLaplace(pump, randomEffectsNodes = "theta")

 

 

CpumpLaplace <- compileNimble(pumpLaplace, project = compileNimble(pump))

 

CpumpLaplace$getNodeNamesVec(returnParams = TRUE)

# [1] "alpha"       "beta"        "meanplot[1]" (OK)

 

CpumpLaplace$calcLaplace(c(1,1,0))

# -Inf (not OK)

 

 

 

library(nimble)

pumpCode <- nimbleCode({
  for (i in 1:N){

    theta[i] ~ dnorm(alpha,exp(beta))

    lambda[i] <- theta[i] * t[i]*exp(meanplot[1])

    #x[i] ~ dpois(4.0) # does not work

x[i] ~ dpois(0*lambda[i]+4.0)

x2[i] ~ dpois(4.0*exp(alpha))

 

  }

  #alpha ~ dexp(1.0)

  #beta ~ dgamma(0.1, 1.0)

alpha ~ dnorm(0.0,1.0)

  beta ~ dnorm(0.0, 1.0)

gamma ~ dnorm(0.0, 1.0)

 

for (i in 1:N)

{meanplot[i]~dnorm(0,1)

}

 

})

 

pumpConsts <- list(N = 10, t = c(94.3, 15.7, 62.9, 126, 5.24, 31.4, 1.05, 1.05, 2.1, 10.5))

 

pumpData <- list(x = c(5, 1, 5, 14, 3, 19, 1, 1, 4, 22), x2 = c(5, 1, 5, 14, 3, 19, 1, 1, 4, 22))

 

pumpInits <- list(alpha = 0.1, beta = 0.1, theta = rep(0.1, pumpConsts$N))

 

pump <- nimbleModel(code = pumpCode, name = "pump", constants = pumpConsts,
                    data = pumpData, inits = pumpInits, buildDerivs = TRUE)

 

pumpLaplace <- buildLaplace(pump, randomEffectsNodes = "theta")

 

 

CpumpLaplace <- compileNimble(pumpLaplace, project = compileNimble(pump))

 

CpumpLaplace$getNodeNamesVec(returnParams = TRUE)

# [1] "alpha"       "beta"        "meanplot[1]" (OK)

 

CpumpLaplace$calcLaplace(c(1,1,0))

# [1] -120.8813

#equal to sum(log(dpois(c(5, 1, 5, 14, 3, 19, 1, 1, 4, 22),4.0*exp(1))))+sum(log(dpois(c(5, 1, 5, 14, 3, 19, 1, 1, 4, 22),4)))

#so both the likelihood of x & x2 are here taken into account

 

 

 

library(nimble)

pumpCode <- nimbleCode({
  for (i in 1:N){

    theta[i] ~ dnorm(alpha,exp(beta))

    lambda[i] <- theta[i] * t[i]*exp(meanplot[1])

    #x[i] ~ dpois(4.0) # does not work

x[i] ~ dpois(0*lambda[i]+4.0)

x2[i] ~ dpois(4.0*exp(gamma+alpha))

 

  }

  #alpha ~ dexp(1.0)

  #beta ~ dgamma(0.1, 1.0)

alpha ~ dnorm(0.0,1.0)

  beta ~ dnorm(0.0, 1.0)

gamma ~ dnorm(0.0, 1.0)

 

for (i in 1:N)

{meanplot[i]~dnorm(0,1)

}

 

})

 

pumpConsts <- list(N = 10, t = c(94.3, 15.7, 62.9, 126, 5.24, 31.4, 1.05, 1.05, 2.1, 10.5))

 

pumpData <- list(x = c(5, 1, 5, 14, 3, 19, 1, 1, 4, 22), x2 = c(5, 1, 5, 14, 3, 19, 1, 1, 4, 22))

 

pumpInits <- list(alpha = 0.1, beta = 0.1, theta = rep(0.1, pumpConsts$N))

 

pump <- nimbleModel(code = pumpCode, name = "pump", constants = pumpConsts,
                    data = pumpData, inits = pumpInits, buildDerivs = TRUE)

 

pumpLaplace <- buildLaplace(pump, randomEffectsNodes = "theta")

 

 

CpumpLaplace <- compileNimble(pumpLaplace, project = compileNimble(pump))

 

CpumpLaplace$getNodeNamesVec(returnParams = TRUE)

# [1] "alpha"       "beta"        "meanplot[1]" (OK)

 

CpumpLaplace$calcLaplace(c(1,1,0))

# [1] -Inf

 

These three problems do not summarize all the issues I am around with buildLaplace. I am trying to find a replicable example for the other ones.

 

Sincerely,

Frédéric

[frederic....@gmail.com]

Dear all,

another issue I have with the buildLaplace function is when some of the random effects that are marginalized on correspond to no data. This may seem strange to you but it may be of use for example to calculate marginalized looIC (which is my case). Here is a replicable example: two models: in the first one observed data correspond to each of the levels of the random effect; in the second one, for theta[2] there is no related observation and the result of $calcLaplace is Infinity:

 

library(nimble)

pumpCode <- nimbleCode({

for (i in 1:2){ theta[i] ~ dnorm(alpha,exp(beta))}

for (j in 1:5) { thetab[j] ~ dnorm(0,1)}

 

  for (i in 1:2){

for (j in 1:5) {

   

    lambda[5*(i-1)+j] <- exp(theta[i]) * t[i]*exp(thetab[j])

    #x[i] ~ dpois(4.0) # does not work

x[5*(i-1)+j] ~ dpois(0*lambda[5*(i-1)+j]+4.0)

 

  }}

 

  #alpha ~ dexp(1.0)

  #beta ~ dgamma(0.1, 1.0)

alpha ~ dnorm(0.0,1.0)

  beta ~ dnorm(0.0, 1.0)

 

})

 

pumpConsts <- list(N = 10, t = c(94.3, 15.7, 62.9, 126, 5.24, 31.4, 1.05, 1.05, 2.1, 10.5))

 

pumpData <- list(x = c(5, 1, 5, 14, 3, 19, 1, 1, 4, 22))

 

pumpInits <- list(alpha = 0.1, beta = 0.1, theta = rep(0.1, pumpConsts$N))

 

pump <- nimbleModel(code = pumpCode, name = "pump", constants = pumpConsts,
                    data = pumpData, inits = pumpInits, buildDerivs = TRUE)

 

pumpLaplace <- buildLaplace(pump, randomEffectsNodes = "theta")

 

 

CpumpLaplace <- compileNimble(pumpLaplace, project = compileNimble(pump))

 

CpumpLaplace$calcLaplace(c(1,1,0,0,0,0,0))

# -63.57501 (OK)

 

 

  

library(nimble)

pumpCode <- nimbleCode({

for (i in 1:2){ theta[i] ~ dnorm(alpha,exp(beta))}

for (j in 1:5) { thetab[j] ~ dnorm(0,1)}

 

  for (i in 1:1){

for (j in 1:5) {

   

    lambda[5*(i-1)+j] <- exp(theta[i]) * t[i]*exp(thetab[j])

    #x[i] ~ dpois(4.0) # does not work

x[5*(i-1)+j] ~ dpois(0*lambda[5*(i-1)+j]+4.0)

 

  }}

 

  #alpha ~ dexp(1.0)

  #beta ~ dgamma(0.1, 1.0)

alpha ~ dnorm(0.0,1.0)

  beta ~ dnorm(0.0, 1.0)

 

})

 

pumpConsts <- list(N = 10, t = c(94.3, 15.7, 62.9, 126, 5.24, 31.4, 1.05, 1.05, 2.1, 10.5))

 

pumpData <- list(x = c(5, 1, 5, 14, 3, 19, 1, 1, 4, 22))

 

pumpInits <- list(alpha = 0.1, beta = 0.1, theta = rep(0.1, pumpConsts$N))

 

pump <- nimbleModel(code = pumpCode, name = "pump", constants = pumpConsts,
                    data = pumpData, inits = pumpInits, buildDerivs = TRUE)

 

pumpLaplace <- buildLaplace(pump, randomEffectsNodes = "theta")

 

 

CpumpLaplace <- compileNimble(pumpLaplace, project = compileNimble(pump))

 

CpumpLaplace$calcLaplace(c(1,1,0,0,0,0,0))

# Inf

 

 Sincerely,

Frédéric

[frederic....@gmail.com]

Dear all,

my last problem at the moment with buildLaplace - but for which I cannot produce a replicable example - is that the names of active parameters produced with the  $getNodeNamesVec(returnParams = TRUE) of the compiled model does not correspond well (in terms of number of parametrs) to what is required by $calcLaplace. In my specific case I can correct for this but may indicate a problem somewhere in the code.

Sincerely,

Frédéric Gosselin

[Perry de Valpine]

Thanks for the questions and reports, Frédéric. I will try to address them in order.

1. Your question about why using

x[i] ~ dpois(4.0)

gives results you didn't expect, and what is documented.

I feel I am not following the question because it looks the same or similar to an earlier question we responded to. nimbleModel creates a graph of relationships among variables. This declaration for x[i] says it has no parent nodes, so the graph may not be what you think it would be. It is correct from the default graphical model perspective that there is nothing to integrate over here, as x[i] is structurally a top-level node or parameter and not a random effect or latent state. The documentation for buildLaplace says that you can call setupMargNodes directly to see how nodes will be handled based on your inputs. The documentation for setupMargNodes is complicated because there are a variety of complicated cases that arise. Again, I may have become lost by the multiple questions so please clarify if I am missing it. You may be asking about parts of the model that should be (or you want to be) included in the likelihood but simply are unrelated to what is being marginalized. Those are called "calcNodesOther" in Laplace, and you can see what setupMargNodes thinks they are, and you can provide them if the default isn't what you want. In your example I am not seeing what you think should be happening, as the x[i] are declared as a set of stand-alone (in effect, predictive) nodes.

2. Your report about crashes when creating and then recreating a model and Laplace algorithm in the same R session: There is a known issue here and it is something we'll try to fix. However I've normally seen that if you rebuild the model (nimbleModel) and recompile it, and then build the Laplace and compile it, I don't see crashes. One thing I see in your code that is unusual to me is:

CpumpLaplace <- compileNimble(pumpLaplace, project = compileNimble(pump)). It makes sense but I would comment that you can do:

CompiledBoth <- compileNimble(pumpLaplace, pump) and then you'll get a list of the two compiled interface objects.

3. Thanks for the report about (if I follow) the length of the node names returned not matching the length of the parameter input to calcLaplace. I will file an issue to look at that.

Perry

To view this discussion on the web visit https://groups.google.com/d/msgid/nimble-users/996fd121-b2c5-49bb-911d-b034a60eca5dn%40googlegroups.com.

[Perry de Valpine]

P.S. Something else to try on #2 would be to remove the name argument to nimbleModel. Please let us know if that seems to help.