Error in mapCopy. Sizes don't match

[Matthijs Hollanders]

Hi all,

I'm trying to implement some simple vectorised distributions, specifically a vectorised normal distribution and a vectorised dLogExp (from {nimbleNoBounds}). dnorm_vec() and dLogExp_vec() are specified as follows (note that mean/sd/rate are not vectors to make them easy to use as prior distributions):

# nimbleFunction for vectorised dnorm
dnorm_vec <- nimbleFunction(
  run = function(x = double(1),
                 mean = double(0, default = 0),
                 sd = double(0, default = 1),
                 log = logical(0, default = 0)) {
    returnType(double(0))
    log_dens <- sum(dnorm(x, mean, sd, log = TRUE))
    if (log)
      return(log_dens)
    else
      return(exp(log_dens))
  }
)

# nimbleFunction for vectorised dLogExp

dLogExp_vec <- nimbleFunction (
  run = function(x = double(1),
                 rate = double(0, default = 1),
                 log = integer(0, default = 0)) {
    returnType(double(0))
    y <- exp(x)
    log_dens <- sum(log(rate) - (rate * y) + x)
    if (log)
      return(log_dens)
    else
      return(exp(log_dens))
  }
)

When I put this distributions to use in some parts of the model, it works fine. For instance, for two log hazard rates, there's no problem at all:

# mortality hazard rates (monthly, juvenile and adult)
  log_h_alpha[1:2] ~ dLogExp_vec(3)
  h_alpha[1:2] <- exp(log_h_alpha[1:2])

For a different part of the model, nimble spits out the error: Error in mapCopy.  Sizes don't match: 2 != 1

This is the part of the model where it shows up (note that the vectorised distributions are commented out—the model works fine without them, but when I apply the vectorised distributions nimble throws hundreds of those errors):

# random individual effects on mortality and capture (bivariate, noncentered)
#  log_epsilon_sigma[1:2] ~ dLogExp_vec(2)
  for (s in 1:2) {
    log_epsilon_sigma[s] ~ dLogExp(2)
  } # s
  epsilon_sigma[1:2] <- exp(log_epsilon_sigma[1:2])
  epsilon_chol[1:2, 1:2] ~ dlkj_corr_cholesky(2, 2)
  epsilon_cor[1:2, 1:2] <- t(epsilon_chol[1:2, 1:2]) %*% epsilon_chol[1:2, 1:2]
 
  # open individual loop
  for (i in 1:n_ind) {
   
    # individual effects cont.
#    epsilon_z[1:2, i] ~ dnorm_vec(0, 1)
    for (s in 1:2) {
      epsilon_z[s, i] ~ dnorm(0, 1)
    }
    epsilon[1:2, i] <- (diag(epsilon_sigma[1:2])
                        %*% t(epsilon_chol[1:2, 1:2])
                        %*% epsilon_z[1:2, i])

What does this error mean and how can I avoid it? Am I making a mistake somewhere?

Cheers,

Matt

[Chris Paciorek]

Hi Matt,

Can you provide a (as simple as possible) reproducible example?

thanks,

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/3fdfdad0-48ab-4cf2-9c38-b1a998da2bbdn%40googlegroups.com.

[Matthijs Hollanders]

Hey Chris,

Sure thing, see attached. Once the z-scores get vectorised the error shows up.

Cheers,

Matt

[Chris Paciorek]

Hi Matt,

It looks like this is because of a subtle issue in your `rnorm_vec` definition. The error is coming from the MCMC initialization.

For a user-defined distribution, `n` in the "r" is the number of random variables to draw. In a multivariate case, this is the number of random vectors (or matrices, as the case might be) to draw, not the number of random elements. Nimble is drawing "n=1" values from your `rnorm_vec` but this is only simulating a single scalar value, not a vector of values. So the number of simulated values that your `rnorm_vec` produces doesn't match what nimble expects (i.e., two values) based on the declaration `z[1:2, i] ~ dnorm_vec(0, 1)`.

You'll need to add an argument to both `dnorm_vec` and `rnorm_vec` that is the length of the vector, since that can't be determined from any of the other arguments. Something like:

rnorm_vec <- nimbleFunction(
  run = function(n = integer(0, default = 1),

                 mean = double(0, default = 0),
                 sd = double(0, default = 1),

                 p = double(0)) {

    returnType(double(1))
    x <- rnorm(p, mean, sd)  # note 'p' not 'n'
    return(x)
  }
)

You'll also need `p` as an arg to `dnorm_vec`, even if you don't use it in the code.

Note that while we require `n` as the first arg of the "r" function, in most cases the "r" functions just assume `n=1`, as will be the case here.

If you're curious, this same issue came up as we wrote nimble's internal `rlkj_corr_cholesky` (seen in the `distributions_implementations.R` file in the nimble source code).

-chris

To view this discussion on the web visit https://groups.google.com/d/msgid/nimble-users/5af45f53-79c8-4f0f-a8b3-e74c15cfab9bn%40googlegroups.com.