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.

[Quresh Latif]

Hi Daniel. I think the functions you wrote will serve my purposes beautifully. Thanks so much!

I'm not sure how to explain the ghost of previous functions in my code if they never existed. I did publish the results from that analysis, and the code I shared is the code that was archived for that study. Additionally, I have a vague memory of my co-worker discovering those functions and pointing them out to me. So, I'm not sure how all this is possible if those functions never existed. I am seeing Nimble version 0.12.1 listed in the supplementary materials for that paper.

I got a bit of a discrepancy between values calculated from your function and those calculated manually from inside the Nimble model. Here is the relevant code:

### Your function ###

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

### GOF code in my model ###

## For this code, 'Traffic[j]' are the data values representing the dependent variable for this model.

 ...

    Traffic_DOY_mn[j] ~ dgamma(shape.Traffic_DOY_mn, rate.Traffic_DOY_mn[j])

    log(pred.Traffic_DOY_mn[j]) <- BETA0.Traffic_DOY_mn +
      BETA.TrailTotm.Traffic_DOY_mn * X.beta[ind.DOY_Speed[j], ind.TrailTotm] +
      BETA.RoadTotm.Traffic_DOY_mn * X.beta[ind.DOY_Speed[j], ind.RoadTotm] +
      BETA.Prp_MotRestricted.Traffic_DOY_mn * X.beta[ind.DOY_Speed[j], ind.Prp_MotRestricted] +
      BETA.Prp_HorseRestricted.Traffic_DOY_mn * X.beta[ind.DOY_Speed[j], ind.Prp_HorseRestricted]
    rate.Traffic_DOY_mn[j] <- shape.Traffic_DOY_mn / pred.Traffic_DOY_mn[j]

    #_____ GOF _____#
    LLobs.Traffic_function[j] <- logdens.gamma(Traffic[j], shape.Traffic, rate.Traffic[j]) # Log-likelihood using your function
    LLobs.Traffic_manual[j] <- log((pow(rate.Traffic[j], shape.Traffic) * # My manual calculation of log-likelihood
      Traffic[j] * exp(-1 * rate.Traffic[j] * Traffic[j]))) -
      loggam(shape.Traffic)
    X.sim.Traffic[j] ~ dgamma(shape.Traffic, rate.Traffic[j])
    LLsim.Traffic_function[j] <- logdens.gamma(X.sim.Traffic[j], shape.Traffic, rate.Traffic[j]) # Log-likelihood for simulated data with your function
    LLsim.Traffic_manual[j] <- log((pow(rate.Traffic[j], shape.Traffic) *   # Log-likelihood for simulated data calculated manually by me
      X.sim.Traffic[j] * exp(-1 * rate.Traffic[j] * X.sim.Traffic[j]))) -
      loggam(shape.Traffic)
    #_______________#

  dev.obs.Traffic_function <- -2 * sum(LLobs.Traffic_function[1:ngrdyrs.hpresent])
  dev.sim.Traffic_function <- -2 * sum(LLsim.Traffic_function[1:ngrdyrs.hpresent])
  test.Traffic_function <- step(dev.sim.Traffic_function - dev.obs.Traffic_function) # GOF test with your function
 
  dev.obs.Traffic_manual <- -2 * sum(LLobs.Traffic_manual[1:ngrdyrs.hpresent])
  dev.sim.Traffic_manual <- -2 * sum(LLsim.Traffic_manual[1:ngrdyrs.hpresent])
  test.Traffic_manual <- step(dev.sim.Traffic_manual - dev.obs.Traffic_manual) # GOF test with my manual calculations

...

If I save "test.Traffic_function" and "test.Traffic_manual", they don't come out quite the same. I'm wondering if you (or anyone) is able to spot any errors in my manual calculations of the log-likelihood. If not, I can try to dig deeper (e.g., save the actual values from both log-likelihood calculations and compare).

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

https://link.edgepilot.com/s/4d97bdc9/yT7STHq8wEuQkzCwdox02A?u=http://www.birdconservancy.org/

To view this discussion on the web visit https://link.edgepilot.com/s/c3d107f4/JbSGWf0McUeYZGUDwf57ZQ?u=https://groups.google.com/d/msgid/nimble-users/CACDDDpACxKC0v6Dt4X2%253Dq5P%252Bop9GgS9PB9BZh0NtMfxEyQMXzA%2540mail.gmail.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/81d27591/HWW95AKriE6NYkdTV-_O2Q?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/88488fb1/NUrxuBaksEm_yJBXgaiguQ?u=https://groups.google.com/d/msgid/nimble-users/CAKbe0hpw1V-uyUbvBKs6OBGMQ2a4BrCxywmfmi1-xYCn7xb7AA%2540mail.gmail.com.

[Quresh Latif]

Hi all. Found my error. I get the same log-likelihood values when I calculate them manually with this:

...

LLobs.Traffic_manual[j] <- log((pow(rate.Traffic[j], shape.Traffic) *

      pow(Traffic[j], shape.Traffic - 1) * exp(-1 * rate.Traffic[j] * Traffic[j]))) -
      loggam(shape.Traffic)

...

Avoiding such errors is why I wanted such a function in the first place, so glad I asked. Thanks so much Daniel!

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

[Daniel Turek]

Wonderful, thanks for following up Qureh, I'm glad you were able to track down the problem.

And yes, your previous use of those functions seems to be an unsolved mystery.  But for now, I'm glad your analysis is working ok.

Cheers,

Daniel

[Quresh Latif]

Well, the mystery may be mostly solved. It seems that JAGS provides those functions, and I did start out that analysis in JAGS, and later moved to Nimble. I'm not sure why they weren't deleted when I moved to Nimble, but at the time, I don't think I was coordinating tracking of my scripts in git as well as I am now, so I suspect I disabled or removed those functions when actually running the analysis but failed to commit those changes to the master repository. Oh well. Hopefully there isn't anyone out there using my scripts as an example and getting confused. :D

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

On Fri, Jul 5, 2024 at 10:39 AM Daniel Turek <danie...@gmail.com> wrote:

Wonderful, thanks for following up Qureh, I'm glad you were able to track down the problem.

And yes, your previous use of those functions seems to be an unsolved mystery.  But for now, I'm glad your analysis is working ok.

Cheers,

Daniel

On Thu, Jul 4, 2024 at 8:07 PM Quresh Latif <quresh...@birdconservancy.org> wrote:

Hi all. Found my error. I get the same log-likelihood values when I calculate them manually with this:

...

LLobs.Traffic_manual[j] <- log((pow(rate.Traffic[j], shape.Traffic) *

      pow(Traffic[j], shape.Traffic - 1) * exp(-1 * rate.Traffic[j] * Traffic[j]))) -
      loggam(shape.Traffic)

...

Avoiding such errors is why I wanted such a function in the first place, so glad I asked. Thanks so much Daniel!

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

https://link.edgepilot.com/s/15209b6c/U77ZPJeOZUaVW1lOL2Gl1A?u=http://www.birdconservancy.org/

https://link.edgepilot.com/s/15209b6c/U77ZPJeOZUaVW1lOL2Gl1A?u=http://www.birdconservancy.org/

To unsubscribe from this topic, visit https://link.edgepilot.com/s/b72f6374/D4EfhaOjTU6wX5QPznXw9w?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/d0c19022/5elGa3NarUu4uP6Ow8wCbg?u=https://groups.google.com/d/msgid/nimble-users/CAKbe0hqXUZc8TdEem_oggLmsLe3iD8cKAMWEzf2f-6CzwB%253Ddiw%2540mail.gmail.com.

[Quresh Latif]

Hi Daniel. I've got a follow-up question. I tried to run a model that uses your function for calculating log densities with parallel processing using my RunNimbleParallel function and I get an error:

"Error in checkForRemoteErrors(lapply(cl, recvResult)) :
  3 nodes produced errors; first error: R function 'logdens.gamma' in the code 'logdens.gamma(Traffic[j], shape.Traffic, rate.Traffic[j])' does not exist."

I'm guessing the logdens.gamma function is not making it into the parallelized clusters. Do I just need to list the "logdens.gamma" in the list of objects to send to each cluster on lines 26-27 of my function ('"clusterExport(cl, c("model", "inits", "data.inp", "constants", "parameters", "ni", "nt"))')?

This is a low priority, so I am just going to turn off the GOF portion of my model for now, and I can report back when I get to trying this, but would be nice to get confirmation or alternative/additional suggestions in the meantime so I know I'm on the right track.

[Daniel Turek]

Quresh, great to hear you're still working on this.

Yes, I think you've diagnosed the problem correctly.  I'll add one caveat, that I'm not entirely certain how the exporting of logdens.gamma - since it's a nimbleFunction, rather than a standard R object or function - will work.  It may work fine, but I'm not certain of that.  In the event that simply including logdens.gamma in your clusterExport (lines 26-27) doesn't play nicely with the exporting mechanism, then a very easy workaround is just to define the logdens.gamma nimbleFunction right inside your parallelization function, for example right after loading the nimble package, on line 39 of your function, before you create the model object or anything else.  I think that would work just fine.

Cheers,

Daniel

To view this discussion on the web visit https://groups.google.com/d/msgid/nimble-users/be631026-34bd-4ca0-b398-74917114828cn%40googlegroups.com.