"dynamic index out of bounds" in HMM

[tenan....@gmail.com]

Hi all,

I have got a series of warnings first from nimbleModel() function and then from runMCMC() function, with a Hidden Markov model.

I wonder whether warnings are due to the way I specify priors for latent states ("z").

Here below the code and the errors. I can provide the data out of the list if needed.

Many thanks for your help,

Simone

# Read in the data:
CH_f <- read.table(file='data.txt',header=T,sep=" ")
head(CH_f)
dim(CH_f)
nind <- dim(CH_f)[1]
nyears <- dim(CH_f)[2]

# Compute the date of first capture for each individual:
first <- NULL
for (i in 1:nind){
temp <- 1:nyears
first <- c(first,min(temp[CH_f[i,]>=1]))}

# Model:
code <- nimbleCode({
  # ELEMENTARY MATRICIES

  # STATE PROCESS:
  # step 1:
  ps1[1,1:8] <- c(phi,0,0,0,0,0,0,1-phi)
  ps1[2,1:8] <- c(0,phi,0,0,0,0,0,1-phi)
  ps1[3,1:8] <- c(0,0,phi,0,0,0,0,1-phi)
  ps1[4,1:8] <- c(0,0,0,phi,0,0,0,1-phi)
  ps1[5,1:8] <- c(0,0,0,0,phi,0,0,1-phi)
  ps1[6,1:8] <- c(0,0,0,0,0,phi,0,1-phi)
  ps1[7,1:8] <- c(0,0,0,0,0,0,phi,1-phi)
  ps1[8,1:8] <- c(0,0,0,0,0,0,0,1)

  # step 2:
  ps2[1,1:8] <-c(1-psi,psi,0,0,0,0,0,0)
  ps2[2,1:8] <-c(0,1-psi,psi,0,0,0,0,0)
  ps2[3,1:8] <-c(0,0,1,0,0,0,0,0)
  ps2[4,1:8] <-c(0,0,0,1,0,0,0,0)
  ps2[5,1:8] <-c(0,0,0,0,1,0,0,0)
  ps2[6,1:8] <-c(0,0,0,0,0,1,0,0)
  ps2[7,1:8] <-c(0,0,0,0,0,0,1,0)
  ps2[8,1:8] <-c(0,0,0,0,0,0,0,1)

  # step 3:
  ps3[1,1:11] <- c(1,0,0,0,0,0,0,0,0,0,0)
  ps3[2,1:11] <- c(0,1,0,0,0,0,0,0,0,0,0)
  ps3[3,1:11] <- c(0,0,phi_c,1-phi_c,0,0,0,0,0,0,0)
  ps3[4,1:11] <- c(0,0,0,0,phi_c,1-phi_c,0,0,0,0,0)
  ps3[5,1:11] <- c(0,0,0,0,0,0,phi_c,1-phi_c,0,0,0)
  ps3[6,1:11] <- c(0,0,0,0,0,0,0,0,1,0,0)
  ps3[7,1:11] <- c(0,0,0,0,0,0,0,0,0,1,0)
  ps3[8,1:11] <- c(0,0,0,0,0,0,0,0,0,0,1)

  # step 4:
  ps4[1,1:11] <- c(1,0,0,0,0,0,0,0,0,0,0)
  ps4[2,1:11] <- c(0,1,0,0,0,0,0,0,0,0,0)
  ps4[3,1:11] <- c(0,0,0,1,0,0,0,0,0,0,0)
  ps4[4,1:11] <- c(0,0,1,0,0,0,0,0,0,0,0)
  ps4[5,1:11] <- c(0,0,0,0,0,1,0,0,0,0,0)
  ps4[6,1:11] <- c(0,0,0,0,1,0,0,0,0,0,0)
  ps4[7,1:11] <- c(0,0,0,0,0,0,0,1,0,0,0)
  ps4[8,1:11] <- c(0,0,0,0,0,0,1,0,0,0,0)
  ps4[9,1:11] <- c(0,0,0,0,0,0,0,0,1,0,0)
  ps4[10,1:11] <- c(0,0,0,0,0,0,0,0,0,1,0)
  ps4[11,1:11] <- c(0,0,0,0,0,0,0,0,0,0,1)

  # step 5:
  ps5[1,1:8] <- c(1,0,0,0,0,0,0,0)
  ps5[2,1:8] <- c(0,1,0,0,0,0,0,0)
  ps5[3,1:8] <- c(0,0,0,0,0,0,1,0)
  ps5[4,1:8] <- c(0,0,0,1,0,0,0,0)
  ps5[5,1:8] <- c(0,0,0,0,0,0,1,0)
  ps5[6,1:8] <- c(0,0,0,0,1,0,0,0)
  ps5[7,1:8] <- c(0,0,0,0,0,0,1,0)
  ps5[8,1:8] <- c(0,0,0,0,0,1,0,0)
  ps5[9,1:8] <- c(0,0,0,0,0,0,1,0)
  ps5[10,1:8] <- c(0,0,gamma,0,0,0,1-gamma,0)
  ps5[11,1:8] <- c(0,0,0,0,0,0,0,1)

  # State process: Matrix product
  ps[1:8,1:8] <- ps1[1:8,1:8] %*% ps2[1:8,1:8] %*% ps3[1:8,1:11] %*% ps4[1:11,1:11] %*% ps5[1:11,1:8]

  # OBSERVATION PROCESS:
  # step 1:
  po1[1,1:8] <- c(1-p,p,0,0,0,0,0,0)
  po1[2,1:8] <- c(1-p,0,p,0,0,0,0,0)
  po1[3,1:8] <- c(1-p,0,0,p,0,0,0,0)
  po1[4,1:8] <- c(1-p,0,0,0,p,0,0,0)
  po1[5,1:8] <- c(1-p,0,0,0,0,p,0,0)
  po1[6,1:8] <- c(1-p,0,0,0,0,0,p,0)
  po1[7,1:8] <- c(1-p,0,0,0,0,0,0,p)
  po1[8,1:8] <- c(1,0,0,0,0,0,0,0)

  # step 2:
  po2[1,1:7] <- c(1,0,0,0,0,0,0)
  po2[2,1:7] <- c(0,1,0,0,0,0,0)
  po2[3,1:7] <- c(0,0,deltaPB,0,0,0,1-deltaPB)
  po2[4,1:7] <- c(0,0,0,1,0,0,0)
  po2[5,1:7] <- c(0,0,0,0,1,0,0)
  po2[6,1:7] <- c(0,0,0,0,1,0,0)
  po2[7,1:7] <- c(0,0,0,0,1,0,0)
  po2[8,1:7] <- c(0,0,0,0,0,deltaNB,1-deltaNB)

  # Observation process: matrix product
  po[1:8,1:7] <- po1[1:8,1:8] %*% po2[1:8,1:7]

  #
  for (i in 1:nind){  
    z[i,first[i]] <- zFirst[i]

    for (j in (first[i]+1):nyears){  # loop over time
      ## STATE EQUATIONS ##
      z[i,j] ~ dcat(ps[z[i,j-1],1:8])
      ## OBSERVATION EQUATIONS ##
      y[i,j] ~ dcat(po[z[i,j],1:7])
    }
  }
 
  # PRIORS
  phi ~ dunif(0, 1)
  psi ~ dunif(0, 1)
  phi_c ~ dunif(0, 1)
  gamma ~ dunif(0, 1)
  p ~ dunif(0, 1)
  deltaPB ~ dunif(0, 1)
  deltaNB ~ dunif(0, 1)
}) # end of model



# Data and constants
constants <- list(nind = nind,
                  nyears = nyears,
                  first = first
                 )

# state at first encounter
zz <- as.matrix(read.table(file='z.txt',header=T,sep=" "))
zz[zz==0] <- NA
zFirst <- NULL
for(i in 1:nind){
    zFirst[i] <- zz[i,first[i]]
}

data = list(y=as.matrix(CH_f+1),zFirst=zFirst)

# initial values for latent states
z_inits <- as.matrix(read.table(file='z_inits.txt',header=T,sep=" "))

# List of initial values:
inits <- list(
  phi=runif(1,0.92,0.99),
  psi=runif(1,0.30,0.37),
  phi_c=runif(1,0.5,0.6),
  gamma=runif(1,0.5,0.6),
  p=runif(1,0.3,0.4),
  deltaPB=runif(1,0.30,0.31),
  deltaNB=runif(1,0.80,0.85),
  z = z_inits)

head(data$y)
     C2004 C2005 C2006 C2007 C2008 C2009 C2010 C2011 C2012 C2013 C2014 C2015
[1,]     1     1     1     1     1     1     1     1     7     1     5     1
[2,]     1     1     7     1     1     1     1     1     1     1     1     1
[3,]     1     7     1     1     1     1     1     1     1     1     5     1
[4,]     1     1     7     4     1     1     1     1     1     1     1     1
[5,]     2     1     1     1     1     4     1     1     1     5     1     1
[6,]     2     2     2     1     1     1     1     1     1     7     7     1
     C2016 C2017 C2018 C2019
[1,]     1     1     1     1
[2,]     1     1     1     1
[3,]     1     1     6     1
[4,]     1     1     1     1
[5,]     1     5     1     6
[6,]     1     1     1     1

head(z_inits)
     C2004 C2005 C2006 C2007 C2008 C2009 C2010 C2011 C2012 C2013 C2014 C2015
[1,]     0     0     0     0     0     0     0     0    NA     3    NA     5
[2,]     0     0    NA     7     7     7     7     7     7     7     7     7
[3,]     0    NA     3     4     7     7     7     7     7     3     4     5
[4,]     0     0    NA    NA     4     5     7     7     7     7     7     7
[5,]    NA     1     2     2     2    NA     7     7     3    NA     5     7
[6,]    NA    NA    NA     2     2     2     3     4     5    NA    NA     7
     C2016 C2017 C2018 C2019
[1,]     7     7     7     7
[2,]     7     7     7     7
[3,]     7     7    NA     7
[4,]     7     7     7     7
[5,]     3    NA     5    NA
[6,]     7     7     7     3

head(zFirst)
[1] 2 7 2 2 1 1

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

Warning: dynamic index out of bounds: dcat(prob = ps[z[i, j_minus_1], 1:8])
Warning: dynamic index out of bounds: dcat(prob = po[z[i, j], 1:7])
Warning: dynamic index out of bounds: dcat(prob = ps[z[i, j_minus_1], 1:8])

Warning: dynamic index out of bounds: dcat(prob = ps[z[i, j_minus_1], 1:8])
Warning: dynamic index out of bounds: dcat(prob = po[z[i, j], 1:7])
(...)

Rmodel$initializeInfo()
Missing values (NAs) or non-finite values were found in model variables: z. 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).

C_model <- compileNimble(Rmodel)

# configure an MCMC:
conf <- configureMCMC(Rmodel,monitors=parameters, print = TRUE)

# build the MCMC:
mcmc <- buildMCMC(conf, enableWAIC = FALSE)

# compile the MCMC:
Cmcmc <- compileNimble(mcmc, project = Rmodel)

# run the model
samples <- runMCMC(Cmcmc, niter=2000, nburnin=1000, thin = 1, nchains=1,
                   samplesAsCodaMCMC = TRUE, summary = TRUE, WAIC = FALSE)

running chain 1...
warning: problem initializing stochastic node z[1, 10]: logProb is -Inf.
warning: problem initializing stochastic node z[3, 3]: logProb is -Inf.
warning: problem initializing stochastic node z[7, 2]: logProb is -Inf.
warning: problem initializing stochastic node z[4, 5]: logProb is -Inf.
warning: logProb of data node y[4, 4]: logProb is -Inf.
warning: logProb of data node y[6, 2]: logProb is -Inf.
warning: problem initializing stochastic node z[8, 4]: logProb is -Inf.
warning: logProb of data node y[8, 3]: logProb is -Inf.
warning: logProb of data node y[9, 3]: logProb is -Inf.
warning: logProb of data node y[10, 3]: logProb is -Inf.
warning: logProb of data node y[16, 3]: logProb is -Inf.
warning: problem initializing stochastic node z[29, 12]: logProb is -Inf.
warning: logProb of data node y[29, 11]: logProb is -Inf.
warning: logProb of data node y[42, 2]: logProb is -Inf.
warning: logProb of data node y[43, 2]: logProb is -Inf.
warning: logProb of data node y[6, 3]: logProb is -Inf.
warning: problem initializing stochastic node z[7, 4]: logProb is -Inf.
warning: problem initializing stochastic node z[10, 5]: logProb is -Inf.
warning: logProb of data node y[10, 4]: logProb is -Inf.
warning: logProb of data node y[14, 4]: logProb is -Inf.
warning: logProb of data node y[35, 14]: logProb is -Inf.
warning: logProb of data node y[42, 3]: logProb is -Inf.
warning: problem initializing stochastic node z[43, 4]: logProb is -Inf.
warning: logProb of data node y[43, 3]: logProb is -Inf.
warning: logProb of data node y[13, 5]: logProb is -Inf.
warning: problem initializing stochastic node z[14, 6]: logProb is -Inf.
warning: logProb of data node y[14, 5]: logProb is -Inf.
warning: logProb of data node y[28, 13]: logProb is -Inf.
warning: logProb of data node y[35, 15]: logProb is -Inf.
warning: problem initializing stochastic node z[10, 7]: logProb is -Inf.
warning: problem initializing stochastic node z[28, 15]: logProb is -Inf.
warning: logProb of data node y[28, 14]: logProb is -Inf.
warning: logProb of data node y[30, 14]: logProb is -Inf.
warning: logProb of data node y[33, 15]: logProb is -Inf.
warning: logProb of data node y[34, 15]: logProb is -Inf.
warning: logProb of data node y[40, 16]: logProb is -Inf.
warning: logProb of data node y[42, 5]: logProb is -Inf.
warning: logProb of data node y[5, 6]: logProb is -Inf.
warning: problem initializing stochastic node z[6, 7]: logProb is -Inf.
warning: logProb of data node y[8, 7]: logProb is -Inf.
warning: problem initializing stochastic node z[29, 16]: logProb is -Inf.
warning: logProb of data node y[34, 16]: logProb is -Inf.
warning: problem initializing stochastic node z[43, 7]: logProb is -Inf.
warning: problem initializing stochastic node z[11, 9]: logProb is -Inf.
warning: problem initializing stochastic node z[18, 10]: logProb is -Inf.
warning: logProb of data node y[18, 9]: logProb is -Inf.
warning: logProb of data node y[26, 15]: logProb is -Inf.
warning: problem initializing stochastic node z[19, 11]: logProb is -Inf.
warning: logProb of data node y[19, 10]: logProb is -Inf.
warning: problem initializing stochastic node z[21, 11]: logProb is -Inf.
warning: logProb of data node y[26, 16]: logProb is -Inf.
warning: problem initializing stochastic node z[43, 9]: logProb is -Inf.
warning: logProb of data node y[43, 8]: logProb is -Inf.
warning: problem initializing stochastic node z[5, 11]: logProb is -Inf.
warning: logProb of data node y[5, 10]: logProb is -Inf.
warning: logProb of data node y[6, 10]: logProb is -Inf.
warning: problem initializing stochastic node z[9, 12]: logProb is -Inf.
warning: logProb of data node y[9, 11]: logProb is -Inf.
warning: problem initializing stochastic node z[13, 12]: logProb is -Inf.
warning: logProb of data node y[13, 11]: logProb is -Inf.
warning: logProb of data node y[18, 12]: logProb is -Inf.
warning: problem initializing stochastic node z[20, 13]: logProb is -Inf.
warning: logProb of data node y[20, 12]: logProb is -Inf.
warning: problem initializing stochastic node z[21, 13]: logProb is -Inf.
warning: logProb of data node y[21, 12]: logProb is -Inf.
warning: problem initializing stochastic node z[42, 12]: logProb is -Inf.
warning: logProb of data node y[42, 11]: logProb is -Inf.
warning: logProb of data node y[43, 11]: logProb is -Inf.
warning: logProb of data node y[8, 13]: logProb is -Inf.
warning: problem initializing stochastic node z[14, 14]: logProb is -Inf.
warning: problem initializing stochastic node z[15, 14]: logProb is -Inf.
warning: logProb of data node y[15, 13]: logProb is -Inf.
warning: problem initializing stochastic node z[8, 15]: logProb is -Inf.
warning: logProb of data node y[8, 14]: logProb is -Inf.
warning: logProb of data node y[3, 15]: logProb is -Inf.
warning: problem initializing stochastic node z[5, 15]: logProb is -Inf.
warning: logProb of data node y[5, 14]: logProb is -Inf.
warning: logProb of data node y[9, 16]: logProb is -Inf.
warning: logProb of data node y[14, 16]: logProb is -Inf.
warning: logProb of data node y[5, 16]: logProb is -Inf.
warning: logProb of data node y[7, 16]: logProb is -Inf.
|-------------|-------------|-------------|-------------|
|-------------------------------------------------------|
Error in quantile.default(x, 0.025) :
  missing values and NaN's not allowed if 'na.rm' is FALSE

[Daniel Turek]

A few comments here.

I think the problems are coming from specifying zFirst as "data", then  using zFirst on the right-hand-side:

 z[i,first[i]] <- zFirst[i]

Generally, for initialization of latent z states such as this, what you want is:

A data structure z_data provided in the "data" list, which has 1's for each z[i,t] up to and including the time t of the final positive sighting of each individual, and NA's thereafter

Another data stricture, z_inits, provided in the "inits" list, which has NA's everywhere that z_data has 1's, and then valid initial values thereafter (perhaps, for example, all 1's thereafter.  Or, alternatively, all 0's.  Or, within each row, some number of 1's, followed exclusively by 0's).

That said, without your actual data and not building your model, it's hard to say for sure.

Another comment, with the large - and very sparse - matrices floating around, not to mention all the matrix multiplications taking place, then I have to feel that the model could be hugely sped up with a little more custom math going on.  Tayloring the model code more specifically to only do the calculations that actually use the non-zero matrix values, and not constructing the large matrices in the first place.  This would take a little careful thinking and coding.  But if efficiency ends up being a problem, or perhaps even prohibitive, then you certainly might think about it.

Daniel

--
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/11ea2e1e-d0f9-44e0-bb91-b7b9c4589eean%40googlegroups.com.

[tenan....@gmail.com]

Dear Daniel,

Thanks for you help. I actually found a bug in the matrices, but I would have two questions more about your comments:

- I think the problems are coming from specifying zFirst as "data", then  using zFirst on the right-hand-side:

 z[i,first[i]] <- zFirst[i]

"zFirst" are values for (known) individual state at first encounter. I have fond this in the NIMBLE User Manual: "Values for nodes that appear only on the right-hand side of BUGS declarations (e.g., covariates/predictors) can be provided as constants or as data or initial values. There is no real difference between providing as data or initial values and the values can be added after building a model via setInits or setData." Should I add them *after* nimbleModel()?

- Another comment, with the large - and very sparse - matrices floating around, not to mention all the matrix multiplications taking place, then I have to feel that the model could be hugely sped up with a little more custom math going on.  Tayloring the model code more specifically to only do the calculations that actually use the non-zero matrix values, and not constructing the large matrices in the first place.  This would take a little careful thinking and coding.  But if efficiency ends up being a problem, or perhaps even prohibitive, then you certainly might think about it.

Wonderful. Could you kindly suggest me some examples or material for improving the code?

Thanks,

Simone

[Daniel Turek]

This is true, everything that's stated in the User Manual.  Nodes that appear only on the RHS (zFirst, here) can be data, or inits, and it won't make a difference if it's data or initial values.  That said, I'm still not sure exactly what's happening with zFirst, and I do feel the observed values of z[,] could be provided as data to the model, again that would be done as:

dataList <- list(...., z = zData)

nimbleModel(...., data = dataList)

where zData is the same size & dimension as z, and contains NA's everywhere, except in the observed locations of z[,], where zData would contain the observed value.

Also, as for avoiding the large matrices, for example, I see that the po[] matrix is constructed as the matrix product of two other (somewhat) large matrics, po1 and po2.  These are fairly sparse.  I would instead suggest constructing the po matrix element-by-element, with most of the elements being zero, and the non-zero elements, defined as:

po[i, j] <- [some expression here involving p, deltaPB, and deltaNB]

This would avoid construction of po1 and po2, and would save a huge number of un-necessary calculations.

Similar could be done with ps[], which is the product of 5 large matrices, so it would really benefit the code to remove all thost multiplications, in a similar way to as I suggested for po[].  I hope this helps.  It will be a little bit of coding, but I believe would speed things up quite a lot.

Hope this helps,

Daniel

To view this discussion on the web visit https://groups.google.com/d/msgid/nimble-users/5d8efb6a-6f5f-4970-87bd-0a643cbd9a87n%40googlegroups.com.