Почему rjags дает здесь несоответствие размеров, взяв подмножество ошибки y?

Я написал эту модель, но rjags дает ошибку несоответствия размеров; Что происходит?

Ошибка в jags.model (textConnection (model1), data = jags_data, n.chains = n_chains,: RUNTIME ERROR: Ошибка компиляции в строке 8. Несоответствие размеров принимает подмножество y

library(rjags)
model1 <- "model {
        C <- 10000
        for (j in 1:nobs){
            zeros[j] ~ dpois(phi[j])
            
            phi[j] <- -log(L[j]) + C
            
            L[j] <- add[j]*(lambda[j]^y[j])*(1-lambda[j])^(1-y[j])
      
            add[j] = ifelse(lambda[j] == 0.5, 2, aux[j])
            aux[j] = 2*arctanh(1 - 2*lambda[j] + 10^(-323))/(1 - 2*lambda[j] + 10^(-323))
            
            logit(lambda[j]) <- inprod(X[j, ], beta)
        }
        beta[1] ~ dnorm(0,1)
        beta[2] ~ dgamma(1,1)
}"


n_chains = 1
n_adapt = 5000
n_iter = 10000
n_thin = 1
n_burnin = 5000

# generate data
n = 100

Ffun = plogis
design_mat = cbind(1, matrix(seq(0,1,by = 0.2), ncol=1))

gen_data = function(n, beta) {
X = design_mat[sample(nrow(design_mat), size = n, replace = T), ]
lambda = Ffun(X %*% beta)
y = rcbern(n,lambda)
idx = is.nan(y)
y[idx] = runif(length(idx))
list(X = X, y = y)
 }

rcbern = function(n,lam){
x = runif(n)
y = log((x*(2*lam-1) - (lam-1))/(1-lam))/log(lam/(1-lam))
return(y)
}

 beta = as.matrix(c(-3, 5))
jags_data = gen_data(n, beta)
jags_data$nobs = n
jg_model <- jags.model(textConnection(model1),
                   data = jags_data,
                   n.chains = n_chains,
                   n.adapt = n_adapt)
update(jg_model, n.iter = n_burnin)
result <- coda.samples(jg_model, 
                   variable.names = c("beta"),
                   n.iter = n_iter, 
                   thin = n_thin,
                   n.chains = n_chains)

beta_est = list(apply(result[[1]],2,median))

person ie86    schedule 13.09.2020    source источник
comment
y - это матрица, поэтому, возможно, попробуйте использовать индекс y[j,] или оставить его как есть и преобразовать y в вектор в соответствующей строке list(X = X, y = as.vector(y)) (непроверено)   -  person user20650    schedule 14.09.2020
comment
@ user20650 Спасибо, заработало. Мне было интересно, можете ли вы также взглянуть на мой другой пост?   -  person ie86    schedule 14.09.2020


Ответы (1)


Как предлагает @ user20650, проблема в том, что вы индексируете y как вектор, а ваши функции генерируются как матрица. Попробуйте этот код с небольшим изменением gen_data():

library(rjags)
model1 <- "model {
        C <- 10000
        for (j in 1:nobs){
            zeros[j] ~ dpois(phi[j])
            
            phi[j] <- -log(L[j]) + C
            
            L[j] <- add[j]*(lambda[j]^y[j])*(1-lambda[j])^(1-y[j])
      
            add[j] = ifelse(lambda[j] == 0.5, 2, aux[j])
            aux[j] = 2*arctanh(1 - 2*lambda[j] + 10^(-323))/(1 - 2*lambda[j] + 10^(-323))
            
            logit(lambda[j]) <- inprod(X[j, ], beta)
        }
        beta[1] ~ dnorm(0,1)
        beta[2] ~ dgamma(1,1)
}"


n_chains = 1
n_adapt = 5000
n_iter = 10000
n_thin = 1
n_burnin = 5000

# generate data
n = 100

Ffun = plogis
design_mat = cbind(1, matrix(seq(0,1,by = 0.2), ncol=1))

gen_data = function(n, beta) {
  X = design_mat[sample(nrow(design_mat), size = n, replace = T), ]
  lambda = Ffun(X %*% beta)
  y = rcbern(n,lambda)
  y <- as.vector(y)
  idx = is.nan(y)
  y[idx] = runif(length(idx))
  list(X = X, y = y)
}

rcbern = function(n,lam){
  x = runif(n)
  y = log((x*(2*lam-1) - (lam-1))/(1-lam))/log(lam/(1-lam))
  return(y)
}

beta = as.matrix(c(-3, 5))
jags_data = gen_data(n, beta)
jags_data$nobs = n
jg_model <- jags.model(textConnection(model1),
                       data = jags_data,
                       n.chains = n_chains,
                       n.adapt = n_adapt)
update(jg_model, n.iter = n_burnin)
result <- coda.samples(jg_model, 
                       variable.names = c("beta"),
                       n.iter = n_iter, 
                       thin = n_thin,
                       n.chains = n_chains)

beta_est = list(apply(result[[1]],2,median))

Вывод:

beta_est
[[1]]
     beta[1]      beta[2] 
-0.006031984  0.692007301 

Вы также можете попробовать y <- y[,1,drop=T] в той же функции вместо as.vector()

person Duck    schedule 14.09.2020
comment
Есть идеи, почему оценки бета-версии настолько ошибочны? - person ie86; 28.09.2020