nimble-dev / nimble

The base NIMBLE package for R
http://R-nimble.org
BSD 3-Clause "New" or "Revised" License
156 stars 23 forks source link

Error in MCMC Parameter Estimation using Nimble: 'object 'i' not found'" #1367

Closed ghost closed 9 months ago

ghost commented 9 months ago

Hello everyone,

I have been struggling with an issue while running an MCMC parameter estimation using Nimble. I am encountering an error message that says "object 'i' not found," even though I have provided the 'I' constant in the code. It's worth mentioning that similar programs have run successfully before.

Here is the code snippet I am working with: `library(nimble)

data

Y = sample(1:5, size = 1000, replace = T) Y = matrix(Y, nrow = 50)

logT = rnorm(1000, 1, 1) logT = matrix(logT, nrow = 50)

data.list = list(Y = Y, logT = logT)

constants

N = nrow(Y) I = ncol(Y)

Q_m1 = rep(1, 20) Q_m2 = c(0, 0, 1, 1, 0, 1, 1, 0, 1, 1, 1, 0, 0, 1, 0, 1, 1, 0, 0, 0)

Q = cbind(Q_m1, Q_m2) sign = Q_m2

K = as.numeric(apply(Y, 2, max)) M = 2 mu_theta = rep(0, M) R = diag(M)

constants <- list(N = N, I = I, K = K, M = M, Q = Q, sign = sign, mu_theta = mu_theta, R = R)

parameters

monitors <- c("theta", "tao","a", "b", "delta", "beta1", "beta2", "sigma")

initial values

sim_theta <- function(N, Phi) { mu <- rep(0, nrow(Phi)) mvtnorm::rmvnorm(N, mu, Phi) }

inits <- list(theta = sim_theta(N, Phi = diag(M)), tao = rnorm(N, 0, 1), a = matrix(abs(rnorm(IM)), nrow = I, ncol = M) Q, b = matrix((rnorm(I*5)), nrow = I), delta = rep(0, I), beta1 = rep(1, I), beta2 = rep(1, I), sigma = rep(0.5, I))

model code

model_code <- nimbleCode({ for (n in 1:N) { for (i in 1:I){ Y[n, i] ~ dcat(prob[n, i, 1:K[i]]) logT[n, i] ~ dnorm(mu_rt[n, i], sigma[i]) } theta[n, 1:M] ~ dmnorm(mu_theta[1:M], R[1:M,1:M]) tao[n] ~ dnorm(0, 1)

# Item category probabilities
for (i in 1:I){
  for (k in 1:K[i]){
    eta[n, i, k] <- sum(a[i, ] * Q[i, ] * (theta[n, ] - b[i, k]))
    psum[n, i, k] <- sum(eta[n, i, 1:k])
    exp.psum[n, i, k] <- exp(psum[n, i, k])
    prob[n, i, k] <- exp.psum[n, i, k] / sum(exp.psum[n, i, 1:K[i]])
  }
  gini[n, i] <- 1 - sum((prob[n, i, 1:K[i]])^2)
}

for (i in 1:I) {
  if (sign[i] == 1){
    we[i] <- a[i, 1]/a[i, 2]
  }
  if (sign[i] == 0){
    we[i] <- 0
  }
  mu_rt[n, i] <- tao[n] + delta[i] + beta1[i] * gini[n, i] + beta2[i] * we[i]
}

}

Priors

for (i in 1:I){ for (m in 1:M) { a[i, m] ~ T(dnorm(1, sd = 1), 0, 3) } }

for (i in 1:I){ b[i, 1] <- 0 for (k in 2:K[i]){ b[i, k] ~ dnorm(0, sd = 1) } }

for (i in 1:I){ delta[i] ~ dnorm(0, sd = 1) beta1[i] ~ dnorm(0, sd = 1) beta2[i] ~ dnorm(0, sd = 1) sigma[i] ~ dinvgamma(1, 1) } })

run mcmc

system.time( fit <- nimbleMCMC(code = model_code, constants = constants, data = data.list, inits = inits, monitors = monitors, niter = 2000, nburnin = 1000, nchains = 2, summary = TRUE)) ` I would greatly appreciate it if someone could help me resolve this issue. I have been stuck on it for quite some time now.

Thank you in advance for your assistance!

paciorek commented 9 months ago

Please post your problem to our nimble-users list (https://groups.google.com/g/nimble-users).