使用R2winBUGS从R获取WinBUGS Leuk示例

Question

我正在尝试从R学习使用winBUGS，我有使用R的经验。我已经成功地从R运行了一个简单的例子，没有任何问题。我一直在尝试运行Leuk:从winBUGS示例中存活第一卷。我已经成功地从winBUGS图形用户界面运行了这一卷，没有出现任何问题。我的问题是尽我所能(我已经尝试和搜索了几天)，我不能让它使用R2winBUGS运行。我相信这是很简单的事情。

如果我尝试在脚本中设置inits，我得到的错误消息是

Error in bugs(data = L, inits = inits, 
      parameters.to.save = params, model.file "model.txt",  : 
      Number of initialized chains (length(inits)) != n.chains

我知道这意味着我还没有初始化一些链，但我粘贴了winbugs示例手册中的初始化代码，所有其他设置对我来说都与在winBUGS图形用户界面上运行时相同。

如果我尝试使用inits=NULL，则会收到另一条错误消息

    display(log)
check(C:/BUGS/model.txt)
model is syntactically correct
data(C:/BUGS/data.txt)
data loaded
compile(1)
model compiled
gen.inits()
shape parameter (r) of gamma dL0[1] too small -- cannot sample
thin.updater(1)
update(500)
command #Bugs:update cannot be executed (is greyed out)
set(beta)

这向我表明，在解决了第一个问题之后，我仍然会有问题！我要放弃使用winBUGS了，有没有人能救救我？我知道我可能会看起来很愚蠢，但每个人都必须学习：-)

我在Windows XP 2002 SP3上使用winBUGS 1.4.3

我的R代码在下面，非常感谢你至少读到这里。

   rm(list = ls())

L<-list(N = 42, T = 17, eps = 1.0E-10,
        obs.t = c(1, 1, 2, 2, 3, 4, 4, 5, 5, 8, 8, 8, 8, 11, 11, 12, 12, 15, 17, 22, 23, 6, 
        6, 6, 6, 7, 9, 10, 10, 11, 13, 16, 17, 19, 20, 22, 23, 25, 32, 32, 34, 35),
        fail = c(1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 1, 0, 1, 0, 0, 1, 1, 0, 0, 0, 1, 1, 0, 0, 0, 0, 0),
        Z = c(0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5,
               -0.5, -0.5, -0.5, -0.5, -0.5, -0.5, -0.5, -0.5, -0.5, -0.5, -0.5, -0.5, -0.5, -0.5, -0.5, -0.5, -0.5, -0.5, -0.5, -0.5, -0.5),
        t = c(1, 2, 3, 4, 5, 6, 7, 8, 10, 11, 12, 13, 15, 16, 17, 22, 23, 35))

### 5.4. Analysis using WinBUGS
library(R2WinBUGS)      # Load the R2WinBUGS library CHOOSE to use WinBUGS
    #library(R2OpenBUGS)        # Load the R2OpenBUGS library CHOOSE to use OpenBUGS
setwd("C://BUGS") 

# Save BUGS description of the model to working directory
sink("model.txt")
cat("

model
{
# Set up data
for(i in 1:N) {
for(j in 1:T) {

# risk set = 1 if obs.t >= t
Y[i,j] <- step(obs.t[i] - t[j] + eps)
# counting process jump = 1 if obs.t in [ t[j], t[j+1] )
# i.e. if t[j] <= obs.t < t[j+1]
dN[i, j] <- Y[i, j] * step(t[j + 1] - obs.t[i] - eps) * fail[i]
}
}
# Model
for(j in 1:T) {
for(i in 1:N) {
dN[i, j] ~ dpois(Idt[i, j]) # Likelihood
Idt[i, j] <- Y[i, j] * exp(beta * Z[i]) * dL0[j] # Intensity
}
dL0[j] ~ dgamma(mu[j], c)
mu[j] <- dL0.star[j] * c # prior mean hazard
# Survivor function = exp(-Integral{l0(u)du})^exp(beta*z)
S.treat[j] <- pow(exp(-sum(dL0[1 : j])), exp(beta * -0.5));
S.placebo[j] <- pow(exp(-sum(dL0[1 : j])), exp(beta * 0.5));
}
c <- 0.001
r <- 0.1
for (j in 1 : T) {
dL0.star[j] <- r * (t[j + 1] - t[j])
}
beta ~ dnorm(0.0,0.000001)
}


",fill=TRUE)
sink()

params<- c("beta","S.placebo","S.treat")

inits<-list( beta = 0.0, 
     dL0 = c(1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,
             1.0,1.0,1.0,1.0,1.0,1.0, 1.0,1.0))

# MCMC settings
nc <-1                  # Number of chains
ni <- 1000              # Number of draws from posterior (for each chain)
ns<-1000 #Number of sims (n.sims)
nb <- floor(ni/2)                   # Number of draws to discard as burn-in
nt <- max(1, floor(nc * (ni - nb) / ns))# Thinning rate
Lout<-list()



# Start Gibbs sampler: Run model in WinBUGS and save results in object called out
out <- bugs(data = L, inits =inits, parameters.to.save = params, model.file = "model.txt", 
n.thin = nt, n.chains = nc, n.burnin = nb, n.iter = ni, debug = T, DIC = TRUE,digits=5,
codaPkg=FALSE, working.directory = getwd())

当我设置nc<-2时，我认为它只是设置每个样本节点的MCMC链的数量，我得到了以下结果

display(log)
check(C:/BUGS/model.txt)
model is syntactically correct
data(C:/BUGS/data.txt)
data loaded
compile(2)
model compiled
inits(1,C:/BUGS/inits1.txt)
this chain contains uninitialized variables
inits(2,C:/BUGS/inits2.txt)
this chain contains uninitialized variables
gen.inits()
shape parameter (r) of gamma dL0[1] too small -- cannot sample
thin.updater(1)
update(500)
command #Bugs:update cannot be executed (is greyed out)

Answer 1

好的，我自己解决了这个问题

我改变的第一件事是，我不喜欢很多变量的名字和r函数't‘相同，所以我把所有东西都改成了唯一的变量名。

其次，我没有为每个链提供唯一的起始值列表

所以inits被改成

inits <- function() list (list (beta1 = rnorm(1)，dL0 = runif(17，0，1))，list (beta1 = rnorm(1)，dL0 = runif(17，0，1)

现在它起作用了。