Laplace approximation within MCMC

[frederic....@gmail.com]

Dear all,

I would like to test the use of the Laplace approximation to marginalize the log posterior density within a MCMC - as I have seen during ISEC 2022 conference. I have tried to do it based on the example in the help file for the function buildLaplace. But here weems to be an error just at the end.

Here is the temptative code - stopped at the place where the Error occurs:

library(nimble)

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

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

    x[i] ~ dpois(lambda[i])

  }

  alpha ~ dexp(1.0)

  beta ~ dgamma(0.1, 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)

 

Cpump<-          compileNimble(pump)

 

Confpump<-configureMCMC(pump, thin =10, monitors = c("alpha","beta"))         
# Build Laplace approximation

pumpLaplace <- buildLaplace(pump)

 

Confpump$removeSamplers(c("alpha","beta"))

Confpump$addSampler(target = "alpha", type = "RW_llFunction",

control = list(llFunction = pumpLaplace, includesTarget = FALSE))

Confpump$addSampler(target = "beta", type = "RW_llFunction",

control = list(llFunction = pumpLaplace, includesTarget = FALSE))

 

ModelMCMCpump <- nimble::buildMCMC(Confpump)

CModelMCMCpump <- nimble::compileNimble(ModelMCMCpump, project =Cpump)

The error message is:

Error: In sizeAssignAfterRecursing: LHS is not in typeEnv or symTab and cannot be added now. This occurred for: modelLP0 <- nfMethod(cppPointerDereference(llFunction),"run")() This was part of the call: { modelLP0 <- nfMethod(cppPointerDereference(llFunction),"run")() if(!(includesTarget)) { modelLP0 <- modelLP0 + getLogProb(nodeFxnVector=model_target_nodeFxnVector_includeData_TRUE_SU_) } propValue <- rnorm(n=1,mean=model_oB_oBtarget_cB_cB[model_oB_oBtarget_cB_cB_flatIndex],sd=scale) nfMethod(my_setAndCalculateOne,"run")(propValue) modelLP1 <- nfMethod(llFunction,"run")() if(!(includesTarget)) { modelLP1 <- modelLP1 + getLogProb(nodeFxnVector=model_target_nodeFxnVector_includeData_TRUE_SU_) } jump <- nfMethod(my_decideAndJump,"run")(modelLP1, modelLP0, 0, 0) if(adaptive) { nfMethod(targetRWSamplerFunction,"adaptiveProcedure")(jump) scale <<- nfVar(targetRWSamplerFunction,"scale") } }

Do you see where my code went wrong? Is there an issue within the package?

I thank you very much in advance,

Frédéric 

[Daniel Turek]

Frédéric, great that you're trying out the new functionality.  Give it a try using the modified code, below.  The main change is that you needed to define a custom nimbleFunction to serve the role of the llFunction for the RW_llFunction_block sampler.  I defined the nimbleFunction as llFun_Laplace, and the specialized instance of it (specialized to the model, and the {alpha, beta} parameters) is called llFun.  This custom function stores new values of alpha and beta into the model, executes the Laplace approximation (calcLaplace method), then returns the (approximate) marginal log-likelihood of alpha and beta.

The code below also removes the default (conjugate) samplers which are assigned to the theta[] elements, since those will be marginalized over in the Laplace approximation.

The code also uses the RW_llFunction_block sampler (rather than the RW_llFunction sampler, as you had), to jointly update both {alpha, beta} together.  I believe that is the necessary approach for the Laplace in MCMC for this model, since {alpha, beta} together, jointly, define the prior distribution for the theta[] elements which are integrated out.

But, hopefully / likely others will have better ideas and suggestions for you, as well.

Cheers,

Daniel

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

Cpump <- compileNimble(pump)


Confpump<-configureMCMC(pump, thin =10, monitors = c("alpha","beta"))

Confpump$removeSamplers(c("alpha","beta","theta"))

llFun_Laplace <- nimbleFunction(
    setup = function(model, paramNodes) {
        laplace <- buildLaplace(model, paramNodes)
    },
    run = function() {
        paramValues <- values(model, paramNodes)
        ll <- laplace$calcLaplace(paramValues)
        returnType(double())
        return(ll)
    }
)

llFun <- llFun_Laplace(pump, paramNodes = c("alpha", "beta"))

Confpump$addSampler(target = c("alpha", "beta"),
                    type = 'RW_llFunction_block',
                    control = list(llFunction = llFun, includesTarget = FALSE))

ModelMCMCpump <- nimble::buildMCMC(Confpump)

CModelMCMCpump <- nimble::compileNimble(ModelMCMCpump, project =Cpump)

samples <- runMCMC(CModelMCMCpump, 1000)

--
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/fae0f1be-3389-42e3-8036-772da8d1cc46n%40googlegroups.com.

[frederic....@gmail.com]

Thank you very much Daniel for both the solution (it works!) and the explanation. Will be useful for others. 

Sincerely,

Frédéric

[frederic....@gmail.com]

Dear Daniel,

just for sharing (and checking): I propose the following function when what we specify are the random effects on which to marginalize: is this OK?

llFun_Laplace <- nimbleFunction(
    setup = function(model, RENodes) {
        laplace <- buildLaplace(model, randomEffectsNodes=RENodes)
paramNodes<-laplace$getNodeNamesVec(returnParams = TRUE)

    },
    run = function() {
        paramValues <- values(model, paramNodes)
        ll <- laplace$calcLaplace(paramValues)
        returnType(double())
        return(ll)
    }
)

Sincerely,

Frédéric

[Daniel Turek]

Frédéric, great that you're making progress on this.  Yes, I think that code will work as you intend, I'll only add that when you provide the randomEffectsNodes argument, this will (in some sense) affect what the parameter nodes are, in that the parameter nodes will now be defined as the parent nodes (parameters) of the randomEffectsNodes - in contrast to the otherwise default set of parameter nodes, being the model's top-level stochastic nodes.  In most cases this designation should logically be what you want, but I wanted to highlight the difference.

Also I should note that Wei is the guru on the Laplace functionality, and he will be more qualified to give a complete answer.

Cheers!

Daniel

To view this discussion on the web visit https://groups.google.com/d/msgid/nimble-users/aea26fb3-9692-4dc9-bbe3-b9823d0aa694n%40googlegroups.com.

[Wei Zhang]

Daniel has already explained this well. One more point I want to add here is that since you are using this within MCMC then you should be careful about setting up the RW_llFunction/RW_llFunction_block samplers. More specifically, when you specify particular random effects nodes to integrate out, their stochastic parent nodes (i.e. paramNodes within Laplace) need to be assigned with a RW_llFunction/RW_llFunction_block sampler for MCMC. However, sometimes the llFun_Laplace function you write here does not include all the log-likelihood needed for the paramNodes, i.e. they may have other stochastic dependencies that are not inside the random effects nodes you specified. Then you need to add the log-likelihood component from those nodes into the custom log-likelihood calculation like this: 

ll <- laplace$calcLaplace(paramValues) + model$calculate(extraDepsofParams)

Actually, we have used the pump example to demonstrate this kind of scenario. Please take a look at the code given here: https://github.com/nimble-dev/nimble/blob/devel/packages/Laplace-examples/pump/pumpNimble.R And by the way in the same place, we have some other Laplace examples. 

Hope this helps, and we are happy to discuss any questions you may have. 

Cheers,

Wei

To view this discussion on the web visit https://groups.google.com/d/msgid/nimble-users/CAKbe0hqiAmrsqLWtqbRf_zkYjqqY39CJX7EEsOS%3DN2s2RjcOAA%40mail.gmail.com.

[frederic....@gmail.com]

Thank you very much Wei and Daniel. Actually I was to write that the code I propose did not allow convergence. I think this is for this reason and due to the way I wrote my model (in which some of the data on which Likelhood should be calculated) do not depend on the random effects I am integarting on. Actually seems to explain some of my problems in a post I made in another discussion: https://groups.google.com/g/nimble-users/c/69byAolqNAo/m/qoejQEbrAAAJ

I will try bot to rewrite my model with the 0*random effect trick (so that all the data should be included) to see if this work. But this is a context specific solution.

Nad will try tyo modify the sampler as well (but unsure at this stage I can have a general solution as I discover the structure behind).

Sincerely,

Frédéric