为什么rjags在这里给出y错误子集而导致维度不匹配？

编程问答 2022-07-06

问题描述

我已经编写了此模型，但是rjags给出了尺寸不匹配错误；发生什么事了？

jags.model（textConnection（model1），data = jags_data，n.chains = n_chains，中的错误： 运行时错误： 第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]],median))

解决方法

正如@ user20650所建议的那样，问题是您正在将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]],median)) 中稍作更改即可尝试以下代码：

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

输出：

y <- y[,drop=T]

您也可以在同一功能中尝试使用as.vector()代替TextFormField

jagsrrrjags