Как преобразовать код winbugs в R

Поскольку я новичок в 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. Кто-нибудь знает, что не так и как это исправить. Любая помощь будет принята с благодарностью.

Спасибо заранее

Гюналь


person Günal    schedule 05.12.2012    source источник
comment
Попробуйте data <- list("n","time","cases","distance")   -  person BenBarnes    schedule 05.12.2012
comment
...и удалить n из inits   -  person BenBarnes    schedule 05.12.2012
comment
Большое спасибо, Бен Барнс. Очень признателен   -  person Günal    schedule 05.12.2012


Ответы (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/")

должно сработать. Проблема была в бета-версии. Ваши бета-версии в коде ОШИБОК варьируются от 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/")

некоторые другие начальные значения (для узлов time[i]) по-прежнему будут генерироваться WinBUGS.

person guyabel    schedule 05.12.2012
comment
Большое спасибо gjabel. Я мог бы быть более благодарен. Еще одну вещь, которую я хотел бы спросить: бета-версии в выводах winbugs и R немного отличаются. Например, бета в R составляет 2,2, 1,7 и 0,0, а в winbugs — 2,32, 1,61 и 0,01458. Есть идеи, почему? Большое спасибо Гюналь - person Günal; 05.12.2012
comment
@Günal Вы отказались от периода выгорания? - person Roman Luštrik; 31.01.2013