1

由于我是 winbugs 新手，我需要在 R 中进行贝叶斯分析。为此，我需要将下面的 winbugs 代码转换为 R 代码：

model{
# model’s likelihood
for (i in 1:n){
time[i] ~ dnorm( mu[i], tau ) # stochastic componenent
# link and linear predictor
mu[i] <- beta0 + beta1 * cases[i] + beta2 * distance[i]
}
# prior distributions
tau ~ dgamma( 0.01, 0.01 )
beta0 ~ dnorm( 0.0, 1.0E-4)
beta1 ~ dnorm( 0.0, 1.0E-4)
beta2 ~ dnorm( 0.0, 1.0E-4)
# definition of sigma
s2<-1/tau
s <-sqrt(s2)
# calculation of the sample variance
for (i in 1:n){ c.time[i]<-time[i]-mean(time[]) }
sy2 <- inprod( c.time[], c.time[] )/(n-1)
# calculation of Bayesian version R squared
R2B <- 1 - s2/sy2
# Expected y for a typical delivery time
typical.y <- beta0 + beta1 * mean(cases[]) + beta2 * mean(distance[])
}
INITS
list( tau=1, beta0=1, beta1=0, beta2=0 )
DATA (LIST)
list( n=25,
time = c(16.68, 11.5, 12.03, 14.88, 13.75, 18.11, 8, 17.83,
79.24, 21.5, 40.33, 21, 13.5, 19.75, 24, 29, 15.35,
19, 9.5, 35.1, 17.9, 52.32, 18.75, 19.83, 10.75),
distance = c(560, 220, 340, 80, 150, 330, 110, 210, 1460,
605, 688, 215, 255, 462, 448, 776, 200, 132,
36, 770, 140, 810, 450, 635, 150),
cases = c( 7, 3, 3, 4, 6, 7, 2, 7, 30, 5, 16, 10, 4, 6, 9,
10, 6, 7, 3, 17, 10, 26, 9, 8, 4) )

这是我使用的 R 代码

library(R2WinBUGS)

time <-      c(16.68,11.5,12.03,14.88,13.75,18.11,8,17.83,79.24,21.5,40.33,21,13.5,19.75,24,29,15.35,19,9.5,35.1,17.9,52.32,18.75,19.83,10.75)
cases <- c(7,3,3,4,6,7,2,7,30,5,16,10,4,6,9,10,6,7,3,17,10,26,9,8,4)
distance <-   c(560,220,340,80,150,330,110,210,1260,605,688,215,255,462,448,776,200,132,36,770,140,810,450,635,150)
n <- length(time)

data <- list("time","cases","distance")
inits <- function(){
list(tau=1,n=25,beta1=rep(0,25),beta2=rep(0.25),beta3=rep(0,25))
}

sim <- bugs(data, inits, model.file = "C:/Users/Gunal/Desktop/dummy/reg.txt",
parameters = c("beta1", "beta2","beta3"),
n.chains = 3, n.iter = 1000,codaPkg = FALSE,
bugs.directory = "D:/PROGRAMLAR/WinBUGS14/",debug=TRUE)

print(sim)

这是上面使用的“reg.txt”文件

model{
# model’s likelihood
for (i in 1:n){
time[i] ~ dnorm( mu[i], tau ) # stochastic componenent
# link and linear predictor
mu[i] <- beta0 + beta1 * cases[i] + beta2 * distance[i]
}
# prior distributions
tau ~ dgamma( 0.01, 0.01 )
beta0 ~ dnorm( 0.0, 1.0E-4)
beta1 ~ dnorm( 0.0, 1.0E-4)
beta2 ~ dnorm( 0.0, 1.0E-4)
# definition of sigma
s2<-1/tau
s <-sqrt(s2)
# calculation of the sample variance
for (i in 1:n){ c.time[i]<-time[i]-mean(time[]) }
sy2 <- inprod( c.time[], c.time[] )/(n-1)
# calculation of Bayesian version R squared
R2B <- 1 - s2/sy2
# Expected y for a typical delivery time
typical.y <- beta0 + beta1 * mean(cases[]) + beta2 * mean(distance[])
}

最后这是winbugs发起的错误：

display(log)
check(C:/Users/Gunal/Desktop/dummy/reg.txt)
model is syntactically correct
data(C:/Users/Gunal/AppData/Local/Temp/RtmpUbKAAJ/data.txt)
data loaded
compile(3)
variable n is not defined
inits(1,C:/Users/Gunal/AppData/Local/Temp/RtmpUbKAAJ/inits1.txt)
command #Bugs:inits cannot be executed (is greyed out)
inits(2,C:/Users/Gunal/AppData/Local/Temp/RtmpUbKAAJ/inits2.txt)
command #Bugs:inits cannot be executed (is greyed out)
inits(3,C:/Users/Gunal/AppData/Local/Temp/RtmpUbKAAJ/inits3.txt)
command #Bugs:inits cannot be executed (is greyed out)
gen.inits()
command #Bugs:gen.inits cannot be executed (is greyed out)
thin.updater(1)
update(500)
command #Bugs:update cannot be executed (is greyed out)
set(beta1)
command #Bugs:set cannot be executed (is greyed out)
set(beta2)
command #Bugs:set cannot be executed (is greyed out)
set(beta3)
command #Bugs:set cannot be executed (is greyed out)
set(deviance)
command #Bugs:set cannot be executed (is greyed out)
dic.set()
command #Bugs:dic.set cannot be executed (is greyed out)
update(500)
command #Bugs:update cannot be executed (is greyed out)
coda(*,C:/Users/Gunal/AppData/Local/Temp/RtmpUbKAAJ/coda)
command #Bugs:coda cannot be executed (is greyed out)
stats(*)
command #Bugs:stats cannot be executed (is greyed out)
dic.stats()

DIC
history(*,C:/Users/Gunal/AppData/Local/Temp/RtmpUbKAAJ/history.odc)
command #Bugs:history cannot be executed (is greyed out)
save(C:/Users/Gunal/AppData/Local/Temp/RtmpUbKAAJ/log.odc)
save(C:/Users/Gunal/AppData/Local/Temp/RtmpUbKAAJ/log.txt)

显然，我没有定义变量 n。有谁知道出了什么问题以及如何解决。任何帮助将不胜感激。

提前谢谢了

居纳尔

4

1 回答 1

5

就像是...

time <-c(16.68,11.5,12.03,14.88,13.75,18.11,8,17.83,79.24,21.5,40.33,21,13.5,19.75,24,29,15.35,19,9.5,35.1,17.9,52.32,18.75,19.83,10.75)
cases <- c(7,3,3,4,6,7,2,7,30,5,16,10,4,6,9,10,6,7,3,17,10,26,9,8,4)
distance <-c(560,220,340,80,150,330,110,210,1260,605,688,215,255,462,448,776,200,132,36,770,140,810,450,635,150)
n <- length(time)
data <- list(n=n,time=time,cases=cases,distance=distance)

sim <- bugs(data, inits=NULL, model.file = "C:/Users/Gunal/Desktop/dummy/reg.txt",
        parameters = c("beta0","beta1", "beta2"),
        n.chains = 3, n.iter = 1000, codaPkg = FALSE,
        bugs.directory = "D:/PROGRAMLAR/WinBUGS14/")

应该管用。测试版中存在问题。BUGS 代码中的测试版从 beta0 运行到 beta2。R 代码中的 beta 版本从 beta1 运行到 beta3。

如果您也想设置初始值（R 代码中用于 beta 的初始值是长度为 25 的向量，而它们应该只是一个数字，如在您的 BUGS 文件中），那么这应该可以工作......

inits <- function(){
  list(tau=1,beta0=0,beta1=0,beta2=0)
}
sim <- bugs(data, inits=inits, model.file = "C:/Users/Gunal/Desktop/dummy/reg.txt",
        parameters = c("beta0","beta1", "beta2"),
        n.chains = 3, n.iter = 1000, codaPkg = FALSE,
        bugs.directory = "D:/PROGRAMLAR/WinBUGS14/")

WinBUGS 仍会生成一些其他初始值（用于 time[i] 节点）。

于 2012-12-05T11:45:53.493 回答