1

我正在尝试学习 WINBUGS，并尝试构建一个小型模型，根据教科书中的示例进行调整（代码如下），假设隐藏的受感染携带者群体具有增长率（“R0”）和去除率（筛选和处理）随时间变化。但是，我倾向于收到一系列错误消息（“扫描的令牌无效或意外”、“无法执行初始化”等）。因此，有更多WINBUGS经验的人能否如此善意地关注我是否在我对WINBUGS的理解上犯了愚蠢的错误？我特别不确定是否可以在 WINBUGS 中完成人口的顺序更新（N.est[t+1] <- N.est[t] + newcases - obs）？提前谢谢了

set.seed(1)
n.years <- 10           # Number of years
N1 <- 30                # Initial population size
rnode <- 0.3
attendance <- 0.4

#generating some data. the model has a hidden population of infected carriers that grows, and is screened periodically (with x% attendance). screened individuals are removed from the population. the observed cases are used in the bayesian framework (BUGS) to infer the hidden parameters


obs <- N <- newcases <- numeric(n.years)
N[1] <- N1
for (t in 1:(n.years-1)){
  newcases[t] <- rbinom(n=1,size=N[t],prob=rnode) #
  obs[t] <- rbinom(1, N[t], attendance) #observed cases
  N[t+1] <- N[t] + newcases[t] - obs[t] }

# Specifying model in BUGS language
sink("ssm.bug")
cat("
    model { 
    # Priors and constraints
    N.est[1] ~ dbin(0.3, 100)    # Prior for initial population size
    rnode ~ dunif(0, 1)          # Prior for R0
    attendance ~ dunif(0, 1)     # Prior for attendance at screening

    # Likelihood
    # State process
    for (t in 1:(T-1)){
    newcases ~ dbin(rnode,N.est[t])
    obs ~ dbin(attendance,N.est[t])
    N.est[t+1] <- N.est[t] + newcases - obs    
    }
    ",fill = TRUE)
sink()

# Bundle data
bugs.data <- list(obs = obs, T = n.years)

# Initial values
inits <- function(){list(rnode=0.5,attendance=0.5, N.est = c(runif(1, 20, 40), rep(NA, (n.years-1))))} 

# Parameters monitored
parameters <- c("rnode", "attendance", "N.est")

# MCMC settings
ni <- 25000
nt <- 3
nb <- 10000
nc <- 3

# Call WinBUGS from R (BRT <1 min)
ssm <- bugs(bugs.data, inits, parameters, "ssm.bug", n.chains = nc, n.thin = nt, n.iter = ni, n.burnin = nb, debug = TRUE, bugs.directory = bugs.dir, working.directory = getwd())

# Define function to draw a graph to summarize results
graph.ssm <- function(ssm, N, y){
  fitted <- lower <- upper <- numeric()
  n.years <- length(y)
  for (i in 1:n.years){
    fitted[i] <- mean(ssm$sims.list$N.est[,i])
    lower[i] <- quantile(ssm$sims.list$N.est[,i], 0.025)
    upper[i] <- quantile(ssm$sims.list$N.est[,i], 0.975)}
  m1 <- min(c(y, fitted, N, lower))
  m2 <- max(c(y, fitted, N, upper))
  par(mar = c(4.5, 4, 1, 1), cex = 1.2)
  plot(0, 0, ylim = c(m1, m2), xlim = c(0.5, n.years), ylab = "Population size", xlab = "Year", las = 1, col = "black", type = "l", lwd = 2, frame = FALSE, axes = FALSE)
  axis(2, las = 1)
  axis(1, at = seq(0, n.years, 5), labels = seq(0, n.years, 5))
  axis(1, at = 0:n.years, labels = rep("", n.years + 1), tcl = -0.25)
  polygon(x = c(1:n.years, n.years:1), y = c(lower, upper[n.years:1]), col = "gray90", border = "gray90")
  points(N, type = "l", col = "red", lwd = 2)
  points(y, type = "l", col = "black", lwd = 2)
  points(fitted, type = "l", col = "blue", lwd = 2)
  legend(x = 1, y = m2, legend = c("True", "Observed", "Estimated"), lty = c(1, 1, 1), lwd = c(2, 2, 2), col = c("red", "black", "blue"), bty = "n", cex = 1)
}

# Execute function: Produce figure 
graph.ssm(ssm, N, y)
4

1 回答 1

1

一般来说，我建议在开始通过 R2WinBUGS（或 R2OpenBUGS）实现自动化之前先深入研究 WinBUGS 或 OpenBUGS 以调试您的模型。

在您上面的模型中，我立即发现您最后缺少一个右括号，即您需要

model { 
# Priors and constraints
N.est[1] ~ dbin(0.3, 100)    # Prior for initial population size
rnode ~ dunif(0, 1)          # Prior for R0
attendance ~ dunif(0, 1)     # Prior for attendance at screening

# Likelihood
# State process
for (t in 1:(T-1)){
  newcases ~ dbin(rnode,N.est[t])
  obs ~ dbin(attendance,N.est[t])
  N.est[t+1] <- N.est[t] + newcases - obs    
}
}

当您点击 Model | OpenBUGS 时，您可以很好地了解这些类型的错误。规格 -> 检查型号。如果无法编译模型，则光标会跳转到发生错误的位置（在模型代码的末尾）。

于 2014-07-23T09:39:47.377 回答