log density functions

[Quresh Latif]

In a previous analysis implemented a couple years ago (might be able to pull up the Nimble version if needed), I was able to calculate the log densities of continuous probability functions (e.g., lines 169, 199, and 212 in here). In the current version, I get an error:

Error in getSymbolicParentNodesRecurse(code, constNames, indexNames, nimbleFunctionNames, : R function 'logdensity.gamma' in the code 'logdensity.gamma(LogTraffic[j], shape.LogTraffic, shape.LogTraffic/pred.LogTraffic[j])' does not exist.

Are these types of functions still available under a different name, or are they now deprecated? I don't find any mention of them in the latest documentation or anywhere on the forum. If these do still exist, are they described in documentation somewhere?

[Quresh Latif]

I have not written any nimble functions before. I tried writing a function to calculate the log density of a gamma distribution:

nimLogDgamma <- nimbleFunction(
  run = function(x = double(0), shape = double(0), rate = double(0)) {
    ans <- log(dgamma(x, shape = shape, rate = rate))
    return(ans)
    returnType(double(0))
  })
CnimLogDgamma <- compileNimble(nimLogDgamma)

When I try to use this in my model, I get an error:

Error in getSymbolicParentNodesRecurse(code, constNames, indexNames, nimbleFunctionNames, : R function 'CnimLogDgamma' has arguments that cannot be evaluated in the code 'CnimLogDgamma(LogTraffic[j], shape.LogTraffic, shape.LogTraffic/pred.LogTraffic[j])'. Either the function must be a nimbleFunction or values for the following inputs must be specified as constants in the model: LogTraffic[j],shape.LogTraffic,c("shape.LogTraffic", "pred.LogTraffic[j]").

Still working through the manual to try and figure out how to specify this function so that it can take in stochastic nodes, but any help in the meantime would be much appreciated.

[Chris Paciorek]

Yes, please see Section 12.2 of the manual. 

A couple brief pointers:

- your set of arguments to your user-defined distribution isn't quite right. The manual discusses this.

- you shouldn't compile the nimbleFunction yourself - nimble will do this automatically if you use it in a model. You just call the original, uncompiled nimbleFunction (`nimLogDgamma`) in the model.

If you still run into problems, please make sure to show us your nimble model code as well.

-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/08de2faf-5000-4b1d-b4b9-af9364ecc3b6n%40googlegroups.com.

[Quresh Latif]

OK, thanks. I am attaching my model code. I'd like to check that my approach to what I am trying to do makes the most of Nimble's functionality. Ultimately, I want to run posterior predictive checks based on deviance for my 5 regression models (one bernoulli and the other four gamma). To get deviance, I need to get the log-likelihoods, and I figure it makes sense to write a function to calculate the log-likelihoods. In the attached version, I've worked out (hopefully correctly) how to get the log-likelihoods as part of the model fit, but seems like it would be more elegant (i.e., at least simpler code) to call a user-defined function to calculate the log-likelihoods. Do you agree, or can you think of a better approach that better leverages all of what Nimble has to offer? If my approach sounds good to you, let me know if you still think the chapter on user-defined distributions has the info I need to move the likelihood calculation into a Nimble function.

Quresh S. Latif
Biometrician

Bird Conservancy of the Rockies

230 Cherry St., Ste. 150, Fort Collins, CO 80521

970-482-1707 (ext. 15)

Connecting people, birds and land

www.birdconservancy.org

To view this discussion on the web visit https://link.edgepilot.com/s/e7bb7523/mXRUl_xKrEKzr8yFq9Kkhw?u=https://groups.google.com/d/msgid/nimble-users/08de2faf-5000-4b1d-b4b9-af9364ecc3b6n%2540googlegroups.com.

--
You received this message because you are subscribed to a topic in the Google Groups "nimble-users" group.
To unsubscribe from this topic, visit https://link.edgepilot.com/s/e320a83d/5aDwrvo9b0ilCj0JP6P6cQ?u=https://groups.google.com/d/topic/nimble-users/dZRvJeKUVdI/unsubscribe.
To unsubscribe from this group and all its topics, send an email to nimble-users...@googlegroups.com.
To view this discussion on the web visit https://link.edgepilot.com/s/53853f55/pyZUjIvDrkePEhpqjSxqXg?u=https://groups.google.com/d/msgid/nimble-users/CA%252B7Ts_pri0bWxZOnhsEM%253DLrBo0G%252BQJrYpOYaZkrGOeAD0r826w%2540mail.gmail.com.

[Daniel Turek]

Quresh, thanks for writing, and your continued use of nimble.  Let me give a few comments, which might help get us thinking about the right path.

I looked at the earlier code you linked.  To the best of my knowledge, nimble has never provided (out-of-the-box) the deterministic model functions that you seem to be using here, for example the logdensity.beta, logdensity.gamma, and logdensity.negbin functions.  Are you able to confirm what version of nimble you were using for that analysis?  Or, did you manually define these functions? (which I don't think is the case, since you also said you have not written any nimble functions before).  In any case, this surprised me, since I don't believe these functions were ever provided, and therefore were never deprecated.

If you want to get that approach running again, I've included a short example below, which shows implementations of logdens.beta and logdens.gamma functions, uses these in a model, and also verifies that the log-densities are correct.  I also added another approach in the code below, which directly monitors internal variables of the model (logProb_x and logProby).  This approach might seem the most straightforward, but it might also become less clear when either 'x' or 'y' are non-scalars (for example, they represent arrays or matrices).  Let me know if you have questions about understanding or adapting this code.  This is one approach you could take, which most closely mimics the original code you referenced.

If you don't want to do this node-by-node in the model (which may be a lot of coding), it would actually be easier to do this for the entire model at once (to record the (sum) model log-density (the sum of the log-density of all stochastic nodes in the model).  This approach could also be adapted to only sum the log-densities of whatever subset of the (stochastic) model nodes you'd like to.  Adapting this code, to only record the (sum) log-density of the predictive data nodes in your posterior predictive checks might be the easiest way to go.  Here is an example of doing that:

Record model sum-posterior-log-density

Let me know when you've had time to look at these, and if you give it a try.  In any case, I'm sure we can find a path to get this working for you.  And the example code for calculating (and saving) model log-densities on a node-by-node basis if directly below.

Cheers,

Daniel

library(nimble)

code <- nimbleCode({
    a ~ dunif(0, 10)
    b ~ dunif(0, 10)
    x ~ dbeta(a, b)
    xdens <- logdens.beta(x, a, b)
    y ~ dgamma(a, b)
    ydens <- logdens.gamma(y, a, b)
})

constants <- list()
data <- list(x = 0.3, y = 0.5)
inits <- list(a = 1, b = 1)

logdens.beta <- nimbleFunction(
    run = function(x = double(), a = double(), b = double()) {
        lp <- dbeta(x, a, b, log = TRUE)
        returnType(double())
        return(lp)
    }
)

logdens.gamma <- nimbleFunction(
    run = function(x = double(), a = double(), b = double()) {
        lp <- dgamma(x, a, b, log = TRUE)
        returnType(double())
        return(lp)
    }
)

Rmodel <- nimbleModel(code, constants, data, inits)

conf <- configureMCMC(Rmodel)
conf$addMonitors(c('xdens', 'ydens'))
conf$addMonitors(c('logProb_x', 'logProb_y'))

Rmcmc <- buildMCMC(conf)

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

set.seed(0)
samples <- runMCMC(Cmcmc, 100)

dbeta (data$x, samples[,'a'], samples[,'b'], log = TRUE) - samples[, 'xdens']
samples[, 'logProb_x'] - samples[, 'xdens']

dgamma(data$y, samples[,'a'], samples[,'b'], log = TRUE) - samples[, 'ydens']
samples[, 'logProb_y'] - samples[, 'ydens']

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