如何WinBUGS軟件代碼轉換至R

Question

由於我是新來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)) 

這就是我用

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) 

和這裏將R代碼是「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[]) 
} 

aused最後這是埃羅r由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。有誰知道什麼是錯的，以及如何解決它。任何幫助將不勝感激。

提前

Günal

Answer 1

喜歡的東西非常感謝......

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代碼中的測試版從beta1運行到beta3。

如果你也想設置初始值（你的R代碼中的初始值是長度爲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]節點）。