MCMC sampler samples just a subset of the model's variable without any warnings/errors

[luchinoprince]

Hello to everyone.

I am a new user of Nimble, but with some experience in other sampling software(Stan for example). I have been trying to sample from a model for the last days, but somehow the model just samples(or reports samples) for just a subset of the variables. The sampler and model do not give any warning/error message for why it is doing so, which makes it very disappointing from a software point of view as I truly have no idea where to turn to solve this problem. The code to generate the model with the toy data is the following:

library(nimble)

extractpi <- nimbleFunction(
  run=function(pis=double(1), idx_ij=integer(0), K=integer(0)){
    returnType(double(1))
    lower_ij <- (K*idx_ij)+1
    upper_ij <- (idx_ij+1)*K
    return(pis[lower_ij:upper_ij])
  }
)

extractpi2 <- nimbleFunction(
  run=function(pi_0=double(1), pi_1=double(1), idx_ij=integer(0)){
    returnType(double(1))
    if (idx_ij != 0) return(pi_1)
    else return(pi_0)
  }
)

dweigthbm <- nimbleFunction(
  run=function(x=double(0), theta_ij=double(0), A_ij=integer(0), log=integer(0, default=0)){
    returnType(double(0))
    if (A_ij != 0){
      val <- dnorm(x, mean=theta_ij, sd=1)
    }
    else{
      if (x != `0){`
        val <- 0.0
      }
      else{
        val <- 1.0
      }
    }
    if (log) return(log(val))
    else return(log)
  }
)

rweigthbm <- nimbleFunction(
  run=function(n=integer(0,default=1), theta_ij=double(0), A_ij=integer(0)){
    returnType(double(0))
    if (A_ij != 0){
    return (rnorm(n=1, mean=theta_ij, sd=1))
    }
    else{
      return(0)
    }
  }
)

model_code2 <- nimbleCode({
  pi_0[1:K] ~ ddirch(alpha[1:K])
  pi_1[1:K] ~ ddirch(beta[1:K])

  for (i in 1:N){
    Z[i] ~ dcat(r[1:nB])
  }
  for (i in 1:N){
    for (j in 1:N){
      idx[i,j] <- M[Z[i], Z[j]]
      #pi <- pis[1:K]
      pi[i,j,1:K] <- extractpi2(pi_0[1:K], pi_1[1:K], idx[i,j])
      s[i,j] ~ dcat(pi[i,j,1:K])
    }
  }

  for (i in 1:N){
    for(j in 1:N){
      A[i,j] ~ dbinom(size=1, prob=B[Z[i], Z[j]])
      W[i,j] ~ dweigthbm(theta[s[i,j]], A[i,j])
    }
  }
})

###################### GENERATION TOY DATA ####################################
K <- 2
nB <- 2
n <- 10
## Percentages of vertices in each block
r = array(1/nB, dim=c(nB))

alpha <- array(c(1.0,2.0))
beta <- array(c(2.0,1.0))

pi_0 <- rdirch(alpha=alpha)
pi_1 <- rdirch(alpha=beta)

theta <- array(c(1,3))
B = matrix(0,nrow=nB, ncol=nB)
B[1,1] = 0.5
B[2,2] = 0.8
B[1,2] = 0.2
B[2,1] = B[1,2]

Z_true = array(0, n)
for (i in 1:n){
  Z_true[i] = rcat(prob=r)
}

A = matrix(0, nrow=n, ncol=n)
W = matrix(0, nrow=n, ncol=n)
s_true = matrix(nrow=n, ncol=n)

for (i in 1:n){
  for (j in 1:i){
    if (Z_true[i]==Z_true[j]){
      s_true[i,j] = rcat(prob=pi_0)
      s_true[j,i] = s_true[i,j]
    }
    else{
      s_true[i,j] = rcat(prob=pi_1)
      s_true[j,i] = s_true[i,j]
    }
    A[i,j] = rbinom(n=1, size=1, prob=B[Z_true[i], Z_true[j]])
    A[j,i] = A[i,j]
    if (A[i,j] == 1){
      W[i,j] = rnorm(n=1, mean=theta[s_true[i,j]], sd=1)
      W[j,i] = W[i,j]
    }
  }
}

M <- 1 - diag(nB)

################## MODEL CREATION AND SAMPLING #######################
data = list(A=A, W=W, B=B, theta=theta, alpha=alpha, beta=beta, r=r, M=M)
inits = list(Z=Z_true, s=s_true, pi_0=as.double(pi_0), pi_1=as.double(pi_1))
constants = list(N=n, K=K, nB=nB)
model <- nimbleModel(code=model_code2, constants = constants, data=data, inits=inits)
conf <- configureMCMC(model)
conf$printSamplers()
#samplerlist <- conf$getSamplers()
mcmc.out <- nimbleMCMC(model, data=data, inits=inits, constants=constants, nchains = 2, niter = 500, summary = TRUE)

The above code copied and pasted in R should run. I get no warning, nor error messages, but when I then look at mcmc.out$samples I have only pi_0, pi_1, Z and no trace of the matrix s, as can be seen from here:

> mcmc.out[["samples"]][["chain1"]]
       Z[1] Z[2] Z[3] Z[4] Z[5] Z[6] Z[7] Z[8] Z[9] Z[10]     pi_0[1]   pi_0[2]   pi_1[1]    pi_1[2]
  [1,]    2    2    1    2    2    2    1    1    1     1 0.149406620 0.8505934 0.4066943 0.59330565
  [2,]    2    2    1    2    2    2    1    1    1     1 0.087870575 0.9121294 0.4639563 0.53604366
  [3,]    2    2    1    2    2    2    1    1    1     1 0.102816744 0.8971833 0.5137899 0.48621006
  [4,]    2    2    1    2    2    2    1    1    1     1 0.041543735 0.9584563 0.3910091 0.60899085

If I check the MCMC configuration I get

> conf <- configureMCMC(model)
===== Monitors =====
thin = 1: pi_0, pi_1, Z
===== Samplers =====
RW_dirichlet sampler (2)
  - pi_0[1:2] 
  - pi_1[1:2] 
categorical sampler (110)
  - Z[]  (10 elements)
  - s[]  (100 elements)

So the model correctly sees the variables, and also their type, but then just avoids sampling or does not report it. Should s also be in monitor? I doubt that is the problem though.

Any idea on how to solve this? I really don't know how to solve it since I have no error messages form Nimble.

Also as I side note for Nimble developers, which I hope you take as constructive criticism. I really want to stress that I find it very bad that a code runs, gives no warning nor error, when clearly there is. So I would try to improve the error/warning messaging.

[danielturek]

@luchinoprince There's no error here. The output displays that categorical samplers are assigned to all 100 elements of s[], it is indeed being updated by the MCMC. The Montiors section of the output also shows you that monitors are assigned (to record the values of) pi_0, pi_1, and Z, but not s, since s is latent. You can add a monitor to s if you want to record those samples also, using conf$addMonitors("s") before calling buildMCMC.

[luchinoprince]

Ok this was faster than expected :). I wasted so much time on it that it is a little embarassing ahahaha. Thanks! Maybe add this in the documentation, as reading trough it did not jump out to me that this was going to happen.

Thanks again!

[perrydv]

Hi @luchinoprince, in case you have not found your way there yet, I'll point you to r-nimble.org and the user manual there. The best place for support requests is the nimble-users list. A good strategy if you're getting started is to reduce the size of your example to something small that still displays the problem or question and post that on nimble-users. In this case, it looks like everything is working fine. Not every dimension is monitored by default because for large models doing so could be huge and undesirable. As @danielturek pointed out, the output displayed did tell you what was being monitored and what not, and did show that samplers were assigned to s. Thanks and good luck.