[R] Output from BRugs Doesn't Match That from OpenBUGS

[ethanarenson]

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.