Hi.
I am trying to analyze with BRugs the Box-Tiao variance components example
in WinBUGS. The output from BRugs,
mean sd MC_error val2.5pc median val97.5pc start sample
sigma2.btw 681.9 1161 10.89 0.7016 253.8 4232 25001 100000
sigma2.with 4266.0 1246 4.92 2480.0000 4057.0 7262 25001 100000
doesn't match the output from WinBUGS,
node mean sd MC error 2.5% median
97.5% start
sample
sigma2.btw 2929.0 2191.0 12.0 263.9 2306.0 8658.0 10001
100000
sigma2.with 2725.0 880.3 4.112 1505.0 2562.0 4897.0 10001
100000
Can anyone offer some insight, or--better yet--a solution to fix this?
Thanks in advance,
Ethan
Here is my code:
############################################
############################################
> rm(list=ls(all=T))
> library(R2WinBUGS)
> library(BRugs)
> library(coda)
> #################################
> dats = list(batches = 6, samples = 5, y = structure(
+ .Data = c(1545, 1440, 1440, 1520, 1580,
+ 1540, 1555, 1490, 1560, 1495,
+ 1595, 1550, 1605, 1510, 1560,
+ 1445, 1440, 1595, 1465, 1545,
+ 1595, 1630, 1515, 1635, 1625,
+ 1520, 1455, 1450, 1480, 1445)
+ , .Dim = c(6,5))
+ )
> ##########################################
> model = function()
+ {
+ for( i in 1 : batches )
+ {
+ mu[i] ~ dnorm(theta, tau.btw)
+ for( j in 1 : samples )
+ {
+ y[i , j] ~ dnorm(mu[i], tau.with)
+ }
+ }
+ theta ~ dnorm(0.0, 1.0E-10)
+ # prior for within-variation
+ sigma2.with <- 1 / tau.with
+ tau.with ~ dgamma(0.005, 0.005)
+
+ # Choice of priors for between-variation
+ # Prior 1: uniform on SD
+ sigma.btw~ dunif(0,100)
+ sigma2.btw<-sigma.btw*sigma.btw
+ tau.btw<-1/sigma2.btw
+ }
> ##########################################
> inits = function()
+ list(theta=1500, tau.with=1, sigma.btw=1)
> ##########################################
> setwd('C:/Users/ethan/Desktop/Dissertation/R/Dyes')
>
> print(dats)
$batches
[1] 6
$samples
[1] 5
$y
[,1] [,2] [,3] [,4] [,5]
[1,] 1545 1555 1605 1465 1625
[2,] 1440 1490 1510 1545 1520
[3,] 1440 1560 1560 1595 1455
[4,] 1520 1495 1445 1630 1450
[5,] 1580 1595 1440 1515 1480
[6,] 1540 1550 1595 1635 1445
> bugs.data(dats)
[1] "data.txt"
> model.file.name = 'C:/Users/ethan/Desktop/Dissertation/R/Dyes/Model.txt'
> model.data.name = 'C:/Users/ethan/Desktop/Dissertation/R/Dyes/data.txt'
> n.chains = 1
> write.model(model, model.file.name)
> modelCheck(model.file.name)
model is syntactically correct
> modelData(model.data.name)
data loaded
> modelCompile(numChains=n.chains)
model compiled
> model.init.name = bugsInits(inits)
> modelInits(rep(model.init.name,n.chains))
Initializing chain 1: initial values loaded but chain contain uninitialized
variables
> modelGenInits()
initial values generated, model initialized
> modelUpdate(25000)
25000 updates took 1 s
> samplesSet(c("sigma2.with","sigma2.btw"))
monitor set for variable 'sigma2.with'
monitor set for variable 'sigma2.btw'
> #dicSet()
> modelUpdate(100000)
100000 updates took 5 s
> samplesStats("*")
mean sd MC_error val2.5pc median val97.5pc start sample
sigma2.btw 681.9 1161 10.89 0.7016 253.8 4232 25001 100000
sigma2.with 4266.0 1246 4.92 2480.0000 4057.0 7262 25001 100000
>
############################################
############################################
--
View this message in context:
http://r.789695.n4.nabble.com/Output-from-BRugs-Doesn-t-Match-That-from-OpenBUGS-tp3932530p3932530.html
Sent from the R help mailing list archive at Nabble.com.
______________________________________________
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
