user-defined distribution based on CDF of multivariate Gaussian

[Henry Scharf]

Many thanks to the developers and countless contributors for all the great work building this useful software.

I am attempting to fit a model for jointly dependent binary response vectors (e.g., presence/absence of multiple species at a site). The model is based on a latent multivariate Gaussian distribution of the same dimension as the binary data vector, where realizations of the binary process arise as element-wise truncations above (success) or below (failure) zero. I have been able to develop an MCMC algorithm for this model using NIMBLE for the case when the latent multivariate Gaussian variables are updated with the rest of the parameters. The approach is analogous to a multivariate Albert and Chib (1993) probit model. However, the MCMC algorithm is, unsurprisingly, inefficient. I would like to try integrating out the latent multivariate Gaussian variables. This involves computing CDFs of multivariate Gaussian distributions with general covariance matrices, which I don't believe is currently supported with base NIMBLE functions. Thus, I am trying to write my own custom distributions for the binary response vectors, relying on existing R packages (tmvtnorm or TruncatedNormal) to do the CDF evaluations. I am using nimbleRcall() to try and make the function pmvnorm() from tmvtnorm available to use in my custom probability distribution. Unfortunately, I get an error when I try to build the nimble model, although the custom density function and sampler I defined appear to work fine when I test them out. The error appears to come from a check within pmvnorm(). 

A toy example of the hierarchical model is shown below. This should reproduce the error I get (Error in checkmvArgs(lower = lower, upper = upper, mean = mean, corr = corr, : ‘lower’ not specified or contains NA). However, as I mentioned above, I think the error is probably coming from the pmvnorm() function, not NIMBLE. I suspect that my issue comes from the way I am trying to use pmvnorm() via nimbleRcall + nimbleFunction. Any insight into this issue would be super helpful. 

Thanks!

--Henry

library(nimble); library(tmvtnorm)

## define multivariate probit function ----

mvprobit <- function(z, mu, sigma, log = 0){

  lb <- z

  ub <- z

  lb[z == 0] <- -Inf

  lb[z == 1] <- 0

  ub[z == 0] <- 0

  ub[z == 1] <- Inf

  as.numeric(pmvnorm(mean = mu, sigma = sigma, lower = lb, upper = ub))

}

## nimble-ize mvprobit ----

nimble_mvprobit <- 

  nimbleRcall(prototype = function(z = double(1), mu = double(1), 

                                   sigma = double(2), log = integer(0, default = 0)){}, 

              Rfun = 'mvprobit', 

              returnType = double(0))

## define multivariate bernoulli based on multivariate Gaussian ----

dmbern <- nimbleFunction(

  run = function(x = double(1), mu = double(1), sigma = double(2), 

                 log = integer(0, default = 0)) {

    returnType(double(0))

    Prob <- nimble_mvprobit(z = x, mu = mu, sigma = sigma, log = log)

    if(log) return(log(Prob))

    else return(Prob)

  })

rmbern <- nimbleFunction(

  run = function(n = integer(0, default = 1), mu = double(1), sigma = double(2)) {

    returnType(double(1))

    chol_mat <- chol(sigma)

    latent <- rmnorm_chol(n = 1, mean = mu, cholesky = chol_mat, prec_param = F) 

    return(latent > 0)

  })

## define the model ----

code <- nimbleCode({

  ## mixture probabilities + response distribution ----

  mean[1:n_species] ~ dmnorm(mean = prior_mean[1:n_species], 

                             prec = prior_prec[1:n_species, 1:n_species])

  for(i in 1:n_observations){

    z[i, 1:n_species] ~ dmbern(mu = means[1:n_species], sigma = cov[1:n_species, 1:n_species])

  }

})

## constants + inits ----

n_observations <- 1e2

n_species <- 2

constants <- list(

  n_observations = n_observations,

  n_species = n_species,

  prior_mean = rep(0, n_species),

  prior_prec = diag(n_species),

  cov = diag(n_species)

)

inits <- list(

  mean = rep(0, n_species)

)

## nimbleModel ----

model <- nimbleModel(code = code, constants = constants, inits = inits)

[Perry de Valpine]

Hi Henry,

I got a slightly different error, but close enough that I think it's related.  (Have you updated nimble recently?)

The output I get is:

> model <- nimbleModel(code = code, constants = constants, inits = inits)

defining model...
Registering the following user-provided distributions: dmbern .
NIMBLE has registered dmbern as a distribution based on its use in BUGS code. Note that if you make changes to the nimbleFunctions for the distribution, you must call 'deregisterDistributions' before using the distribution in BUGS code for those changes to take effect.
Adding prior_mean, prior_prec, cov as data for building model.
building model...
setting data and initial values...
running calculate on model (any error reports that follow may simply reflect missing values in model variables) ... Error in checkmvArgs(lower = lower, upper = upper, mean = mean, corr = corr,  :
  ‘mean’ is not a numeric vector

checking model sizes and dimensions... This model is not fully initialized. This is not an error. To see which variables are not initialized, use model$initializeInfo(). For more information on model initialization, see help(modelInitialization).
model building finished.

The "any error reports that follow may simply reflect missing values in model variables" is important.  It means you don't necessarily need to worry about the reported error message if it simply results from uninitialized variables.  The function did exit successfully and a model was created.  I following the instructions to do model$initializeInfo() and got the following:

> model$initializeInfo()
Missing values (NAs) or non-finite values were found in model variables: z, means. This is not an error, but some or all variables may need to be initialized for certain algorithms to operate properly. For more information on model initialization, see help(modelInitialization).

So it seems that z and means were not initialized (did you intend to have one variable called "mean" and another called "means"?).

So far, no problems, just incomplete initialization.

I did check to see if your model would compile and note the following:

In rmbern, it won't successfully cast a logical vector to a double vector return type, so please use:

   result <- numeric(length = length(latent), value = latent > 0)
   return(result)

   F is not recognized for FALSE, so please use "prec_param = FALSE"

With those changes I was able to compile the model.  I stopped there.

-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 post to this group, send email to nimble...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/nimble-users/44cecd71-92dd-4717-8cf7-ee8fc35deed5%40googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

[Henry Scharf]

I'm not sure how I managed to put this typo in there, but the variable 'means' in the dmbern() distribution for z should be 'mean'. With the typo, I got the same error you did. With that fixed, I get the same error I mentioned in the first post. That said, your point that I 'don't necessarily need to worry about the reported error message' is well taken. In fact, everything seems just fine. I am able to simulate the data z, and after doing so model$calculate() works without the error.

For what it's worth, I am at least using nimble 0.8.0, even if I can't read NIMBLE output very well.

Many thanks for taking a look!

Henry

To unsubscribe from this group and stop receiving emails from it, send an email to nimble...@googlegroups.com.

[Perry de Valpine]

Great - good luck with it!

To unsubscribe from this group and stop receiving emails from it, send an email to nimble-users...@googlegroups.com.

To post to this group, send email to nimble...@googlegroups.com.

To view this discussion on the web visit https://groups.google.com/d/msgid/nimble-users/6c78ad06-ddee-4f0a-af68-1ffa7c3641c5%40googlegroups.com.