Some parameters don't update in parametric stochastic block model

Hello to everyone!

I created a relatively simple code implementing a parametric weighted stochastic block model to kind of play around and get familiar with Nimble. The model code is the following

model_code_param_ext <- nimbleCode({

  ## Prior for probability for vertices to be in blocks
  r[1:nB] ~ ddirch(pc[1:nB])

  ## Prior for connectivity within blocks + edges parameter distr.
  for (a in 1:nB){
    for (b in 1:nB){
      B[a,b] ~ dbeta(beta1, beta2)
      ## We assume same hyperparameters for shapes and rates.
      shapes[a,b] ~ dgamma(nu, beta)
      rates[a,b] ~ dgamma(nu, beta)
    }
  }

  ###### Vertex block membership
  for (i in 1:N){
    Z[i] ~ dcat(r[1:nB])
  }

  ##### Now we go to the likelihood of the adiecency matrix and of the weights.
  for (i in 1:N){
    for(j in 1:N){
    #for(j in 1:i-1){
      A[i,j] ~ dbinom(size=1, prob=B[Z[i], Z[j]])
      W[i,j] ~ dweigthbm_gamma(shapes[Z[i], Z[j]], rates[Z[i], Z[j]], A[i,j])
    }
  }
})

with dweigthbm_gamma

dweigthbm_gamma <- nimbleFunction(
  run=function(x=double(0), shape_ij=double(0), rate_ij=double(0), A_ij=integer(0), log=integer(0, default=0)){
    returnType(double(0))
    if (A_ij != 0){
      val <- dgamma(x, shape=shape_ij, rate=rate_ij)
    }
    else{
      if (x != 0){
        val <- 0.0
      }
      else{
        val <- 1.0
      }
    }
    if (log) return(log(val))
    else return(log)
  }
)

rweigthbm_gamma <- nimbleFunction(
  run=function(n=integer(0,default=1), shape_ij=double(0), rate_ij=double(0), A_ij=integer(0)){
    returnType(double(0))
    if (A_ij != 0){
      return (rgamma(n=1, shape=shape_ij, rate=rate_ij))
    }
    else{
      return(0)
    }
  }
)

I then create a slightly misspecified dataset to test the sample in the following way:

K <- 2
## Number of blocks
nB <- 2
## Number of vertices in the graph
n <- 50

B = matrix(0,nrow=nB, ncol=nB)
B[1,1] = 0.5
B[2,2] = 0.7
B[1,2] = 0.2
B[2,1] = B[1,2]
shapes <- array(c(100.0,50.0))
rates <- array(c(100.0, 10.0))

pi_0 <- array(c(0.75, 0.25))
pi_1 <- array(c(0.25, 0.75))

## Generate the data according to the model
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 2:n){
  for (j in 1:(i-1)){
    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] = rgamma(n=1, shape=shapes[s_true[i,j]], rate=rates[s_true[i,j]])
      W[j,i] = W[i,j]
    }
  }
}

where the misspecification is in the weights. I then compile the model and sample

pc <- 3 * array(1, dim=c(nB))

## Seem to give diffuse enough priors
nu <- 10
beta <- 4
beta1 <- 1.5
beta2 <- 1.5

data_pext = list(A=A, W=W)
constants_pext = list(N=n, nB=nB, pc=pc , nu=nu, beta=beta, beta1=beta1, beta2=beta2)
Rmodel_param <- nimbleModel(code=model_code_param_ext, constants = constants_pext, data=data_pext)
conf_param <- configureMCMC(Rmodel_param, monitors=c("Z", "B", "shapes", "rates", "r"))
Rmcmc_param <- buildMCMC(conf_param)
Cmodel_param <- compileNimble(Rmodel_param)
Cmcmc_param <- compileNimble(Rmcmc_param, project = Rmodel_param)

mcmc.out_param <- runMCMC(Cmcmc_param, niter=2000, nburnin=0, nchains=1, summary=TRUE)

During the sampling I get no error messages/warnings, and when I look at the output I get:

                  Mean    Median    St.Dev. 95%CI_low 95%CI_upp
B[1, 1]      0.3757850 0.3725258 0.01902426 0.3456662 0.4180413
B[2, 1]      0.4485430 0.4367798 0.11018092 0.2109650 0.7865345
B[1, 2]      0.4477253 0.4370324 0.11126312 0.2121075 0.8060270
B[2, 2]      0.6186998 0.6378723 0.14178728 0.1740158 0.8263060
Z[1]         1.8010000 2.0000000 0.39934789 1.0000000 2.0000000
Z[2]         1.0000000 1.0000000 0.00000000 1.0000000 1.0000000
Z[3]         1.0000000 1.0000000 0.00000000 1.0000000 1.0000000
Z[4]         1.0000000 1.0000000 0.00000000 1.0000000 1.0000000
Z[5]         1.0000000 1.0000000 0.00000000 1.0000000 1.0000000
Z[6]         1.0000000 1.0000000 0.00000000 1.0000000 1.0000000
Z[7]         1.0000000 1.0000000 0.00000000 1.0000000 1.0000000
Z[8]         1.0000000 1.0000000 0.00000000 1.0000000 1.0000000
Z[9]         1.8050000 2.0000000 0.39629979 1.0000000 2.0000000
Z[10]        1.8115000 2.0000000 0.39120873 1.0000000 2.0000000
Z[11]        1.0000000 1.0000000 0.00000000 1.0000000 1.0000000
Z[12]        1.0000000 1.0000000 0.00000000 1.0000000 1.0000000
Z[13]        1.0000000 1.0000000 0.00000000 1.0000000 1.0000000
......
......

Hence most of the latent variables have not moved during the sampling, which seems very very odd. Does anyone have an idea of what could cause this behaviour?

nimble-dev / nimble

Some parameters don't update in parametric stochastic block model #1353