[Warning] Dynamic index out of bounds_Help

[ARJUN BANIK]

Hi All,

I am getting a series of warnings while running the nimbleMCMC() function, but nimbleModel is working well, and $calculate is returning me a value. Can you please help me understand why I am getting this warning "[Warning] Dynamic index out of bounds: dcatt(prob = pt[g[i, t_minus_1], 1:4, i, t])"?

I have provided my code below. Have also attached a small sample datasets which I am using for testing my code. 

Many thanks for your help,

Arjun

--

library(nimble)
library(MCMCvis)

####################################################################
## ------------------Custom Distribution------------------------- ##
####################################################################

dcatt <- nimbleFunction(
  run = function(x = integer(0),
                 prob = double(1),
                 log = logical(0)) {
    returnType(double(0))
   
    # if (x < 1 || x > length(prob)) {
    #   if (log) return(-Inf)
    #   else
    #     return(0)
    # }
   
    if (prob[x] == 0) {
      if (log)
        return(0)
      else
        return(prob[x])
    }
    if (prob[x] != 0) {
      if (log) {
        return(log(prob[x]))
      } else {
        return(prob[x])
      }
    }
  }
)

rcatt <- nimbleFunction(
  run = function(n = integer(0), prob = double(1)) {
    returnType(integer(0))   # scalar integer
    if(n != 1) stop("rcatt only allows n = 1")
   
    # normalize
    prob <- prob / sum(prob)
    K <- length(prob)
   
    # build cumulative probabilities manually
    cum_prob <- numeric(K)
    cum_prob[1] <- prob[1]
    for(k in 2:K) {
      cum_prob[k] <- cum_prob[k-1] + prob[k]
    }
   
    # sample
    u <- runif(1, 0, 1)
    cat <- 1
    for(k in 1:K) {
      if(u > cum_prob[k]) cat <- cat + 1
    }
   
    return(as.integer(cat))  # scalar integer
  }
)

registerDistributions(list(
 
  dcatt = list(
    BUGSdist = "dcatt(prob)",
    Rdist = "rcatt(prob)",
    types = c("value = integer(0)", "prob = double(1)")
  )
))

####################################################################
## ---------------------------CODE------------------------------- ##
####################################################################

get_first_non_zero_index <- function(row) {
  non_zero_indices <- which(row != 0)
  if (length(non_zero_indices) > 0) {
    return(non_zero_indices[1])
  } else {
    return(NA)
  }
}

cap <- read.csv(file.choose(),header = T)
alive <- read.csv(file.choose(), header= T) # Partially Latent
tag <- read.csv(file.choose(), header = T) # Partially Latent

retag_1 <- read.csv(file.choose(), header = T) # When retagging happen for tag 1 (left tag)
retag_2 <- read.csv(file.choose(), header = T) # When retagging happen for tag 2 (right tag)

code_m1 <- nimbleCode({
  # Priors
  phi ~ dbeta(1,1) # Survival Probability
  p ~ dbeta(1,1) # Capture Probability
 
  lam_1 ~ dbeta(1,1) # Tag 1 Retention Probability
  lam_2 ~ dbeta(1,1) # Tag 2 Retention Probability
 
  # Latent state process probabilities
  px[1, 1] <- phi
  px[1, 2] <- 1 - phi
 
  px[2, 1] <- 0
  px[2, 2] <- 1
 
  # Observation process probabilities
  po[1, 1] <- 1 - p
  po[1, 2] <- p
 
  po[2, 1] <- 1
  po[2, 2] <- 0
 
  # Tag state process probabilities
  # probabilities of state z(t+1) given z(t)
  for (i in 1:N) {
    for (t in f_1[i]:nyears) {
       
      pt[1,1,i,t] <- lam_1*lam_2
      pt[1,2,i,t] <- lam_1*(1-lam_2)
      pt[1,3,i,t] <- (1-lam_1)*lam_2
      pt[1,4,i,t] <- (1-lam_1)*(1-lam_2)
     
      pt[2,1,i,t] <- lam_1*r_2[i,t]*lam_2
      pt[2,2,i,t] <- lam_1*((1-r_2[i,t])+(r_2[i,t]*(1-lam_2)))
      pt[2,3,i,t] <- (1-lam_1)*r_2[i,t]*lam_2
      pt[2,4,i,t] <- (1-lam_1)*((1-r_2[i,t])+(r_2[i,t]*(1-lam_2)))
     
      pt[3,1,i,t] <- r_1[i,t]*lam_1*lam_2
      pt[3,2,i,t] <- r_1[i,t]*lam_1*(1-lam_2)
      pt[3,3,i,t] <- ((1-r_1[i,t])+(r_1[i,t]*(1-lam_1)))*lam_2
      pt[3,4,i,t] <- ((1-r_1[i,t])+(r_1[i,t]*(1-lam_1)))*(1-lam_2)
     
      pt[4,1,i,t] <- r_1[i,t]*lam_1*r_2[i,t]*lam_2
      pt[4,2,i,t] <- r_1[i,t]*lam_1*((1-r_2[i,t])+(r_2[i,t]*(1-lam_2)))
      pt[4,3,i,t] <- ((1-r_1[i,t])+(r_1[i,t]*(1-lam_1)))*r_2[i,t]*lam_2
      pt[4,4,i,t] <- ((1-r_1[i,t])+(r_1[i,t]*(1-lam_1)))*((1-r_2[i,t])+(r_2[i,t]*(1-lam_2)))
    }
  }
 
  #Likelihood 2 (from 2nd year of capture)
  for (i in 1:N){
    a[i, f_1[i]] ~ dcat(inprob_a[1:2])
    g[i, f_1[i]] ~ dcatt(inprob_g[1:4])
   
    for (t in seq(f_1[i] + 1, nyears,
                  length.out = max(0, nyears - f_1[i]))) {
     
      a[i,t] ~ dcat(px[a[i, t - 1], 1:2])
      g[i,t] ~ dcatt(pt[g[i, t - 1], 1:4,i,t])
      h[i,t] ~ dcat(po[a[i, t], 1:2])
     
    }
  }
})

my.data <- list(h = cap + 1, a = alive, g = tag,
                r_1 = retag_1, r_2 = retag_2)

N = nrow(cap)
nyears = ncol(cap)

a_inits <- alive
a_inits[is.na(alive)] <- 1
a_inits[!is.na(alive)] <- NA

g_inits <- tag
#g_inits[is.na(tag)] <- 1
g_inits[is.na(tag)] <-  sample(1:4, sum(is.na(tag)), replace = TRUE)
g_inits[!is.na(tag)] <- NA

f_1 <- apply(cap, 1, get_first_non_zero_index)

for (i in 1:N) {
  if (f_1[i] > 1) a_inits[i, 1:(f_1[i]-1)] <- NA
  if (f_1[i] > 1) g_inits[i, 1:(f_1[i]-1)] <- NA
 
}

my.constants <- list(N = N, nyears = nyears,
                     f_1 = f_1, inprob_a=c(1, 0),
                     inprob_g=c(0.25, 0.25, 0.25, 0.25))

initial.values <- list(phi = runif(1,0,1),
                       p = runif(1,0,1),
                       lam_1 = runif(1,0,1),
                       lam_2 = runif(1,0,1),
                       a = a_inits,
                       g = g_inits)

parameters.to.save <- c("phi", "p", "lam_1", "lam_2")

#parameters.to.save <- c("phi", "p", "lam_1", "lam_2", "g_1", "g_2", "h")

n.iter <- 10
n.burnin <- 2
n.chains <- 3

m <- nimbleModel(code = code_m1,
                 constants = my.constants,
                 data = my.data,
                 inits = initial.values)
m$calculate()

mcmc.multistate <- nimbleMCMC(code = code_m1,
                              constants = my.constants,
                              data = my.data,              
                              inits = initial.values,
                              monitors = parameters.to.save,
                              niter = n.iter,
                              nburnin = n.burnin,
                              nchains = n.chains,
                              summary = TRUE, WAIC = TRUE)

[Perry de Valpine]

Hi Arjun,

Thanks for providing a full example. I haven't run your code yet but instead will suggest some changes from reading the code. If you are still stuck, someone can run your code and look deeper.

The "dynamic index" in the warning is g[i, t_minus_1], so there is somehow a value of that <1 or >4.

All variables in models are doubles, so please change from integer(0) to double(0) at various places in dcatt and rcatt and remove as.integer.

if (prob[x] == 0) and log is TRUE, you should return -Inf, not 0. Returning log prob = 0 means the probability is 1 and a value of 0 being proposed by an MCMC sampler can be accepted. Returning log prob = -Inf means the probability is 0 and the value can never be accepted.

Some tricks for catching the problem in action:

If the warnings appear relatively quickly, then you could feasibly run the MCMC *uncompiled* (outside of nimbleMCMC, after buildMCMC) and insert debugging code into dcatt. For example:

if(x < 1 || x > 4) {

  print("We caught a problem case!")

  browser()

}

and then inspect the local variables inside dcatt as well as any variables in the (uncompiled) model. Or, you could even more simply do:

print(x) 

to see what values of x are being considered by the MCMC samplers (and hence passed as values to dcatt).

If it takes too long to catch an error running uncompiled, you can also catch it while running compiled. The first thing you can do is:

if(x < 1 || x > 4) {

 print("I found a problem with x = ", x)

}

which will allow you to see the problem values of x.

The next level would be to define a nimbleRcall, say "my_debugger", which calls an R function, say "Rmy_debugger", which contains a browser(). Then insert

if(x < 1 || x > 4) {

  print("We caught a problem case!")

  my_debugger()

}

When an invalid value of x if found, control will be passed back to R into my_debugger, which will then browser() to give you an R prompt, from which you can inspect model variables (in the compiled model, since this trick would be for a compiled case).

There are also two changes you can make in registerDistributions. Here is the input list for dcat from the nimble source code:

    dcat    = list(BUGSdist = 'dcat(prob)',
                   Rdist    = 'dcat(prob)',
                   altParams= c('k = length(prob)'),
                   types    = c('prob = double(1)'),
                   range    = c(1, Inf),
                   discrete = TRUE)

So you can give it a range from 1 to Inf, which at least tells it not to go below 1. In your case that might mask your bug so it may or may not be what you want to do to get to the bottom of things. The upper bound Inf is of course not right but nimble does not support ranges that are dynamic (depending on the data or parameters), so there is no upper bound provided for dcat (but you can still check based on length(prob) inside dcatt). You can also tell it discrete=TRUE, which in lieu of having actual integers in models serves to tell the MCMC system something about samplers to use. Without that, you may be getting non-integer proposals and your dcatt may not be working at all correctly in those cases.

So the three essential things are: always use doubles for model variables, always return -Inf for a log density representing log(0), and register the distribution as discrete.

HTH

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 view this discussion visit https://groups.google.com/d/msgid/nimble-users/10e9eccb-da1e-4338-8c44-5dae0d30ad4cn%40googlegroups.com.