跪求赐教！关于winbugs

这样一个问题的后验密度用mcmc的方法在winbugs中如何编程？研究了好几天了，不知从何处下手，请高人不吝赐教啊！急火攻心，食不下咽，请帮帮忙！ problem.doc (74 KB)

2010-2-25 10:13:32 上传

也许远水已经救不了近火,
不过还是提出来供参考.
由电脑产生一组n=10,000
的matched-pair table
        2536        2522
        2490        2452

########
model {
   ## Multinomial distribution
   x[1:4] ~ dmulti(p[],n)
   ## Dirichlet prior distributions
   p[1:4] ~ ddirch(alpha[])
   n <- sum(x[])
}
list(x=c( 2536, 2522,2490 ,2452),alpha = c(1,1,1,1))

在线等~~~~~

为什么？！！！！救命啊！

看来只能实施自救

谢谢你，朋友！不过问题没有完全解决，能不能请你再帮忙解答几个问题？！
1、由电脑产生一组n=10,000的matched-pair table，具体是怎样产生呢？可以详细说一下吗？而且是不是不这样产生数据的话是不行的，比如我任意写x=c（20，12，2，16）
2、你写的这段代码中，运行后可算出p1,p2,p3,p4的后验信息，如果我只想要rd=p[1]-p[3] 的后验信息，那么我加了
model {
   ## Multinomial distribution
   x[1:4] ~ dmulti(p[],n)
   ## Dirichlet prior distributions
   p[1:4] ~ ddirch(alpha[])
   n <- sum(x[])
   rd<-p[1]-p[3]
}
list(x=c( 2536, 2522,2490 ,2452),alpha = c(1,1,1,1))
为什么提示rd是multiple definition呢？
你好人做到底，再帮帮我，万分感谢！

另外，你是否知道如何计算HPD interval，因为winbugs计算出来的是5%，95%等等的置信区间，也就是所谓的等尾的置信区间，但没给出HPDinterval，如果要算这个，winbugs怎么实现呢？

n=10时,误差比较大.所以我取n=10,000.
matched-pair table,我是在matlab产生.
由两组vector相互比较而得.
看了下例你就明白

11

00

10

10

10

11

11

00

11

00


matched-pair=


43

03

HPDinterval可在 R package coda 取得
先安装package "BRugs","coda"
程序修改为testmodel.txt & testdata.txt
两个文件放在My Documents底下让R能读到
底下是R的执行程序及输出结果:
########
library(BRugs)
modelCheck("testmodel.txt")
# check model file
modelData("testdata.txt")
# read data file
modelCompile(numChains=1)
# compile model with 1 chains
modelGenInits()
modelUpdate(1000)
# burn in
samplesSet(c("p","rd"))

modelUpdate(10000)
# 10000 more iterations ....

samplesStats("*")
# the summarized results
## density plots
samplesDensity("*")
# plot the densities
## Write output to coda files
samplesCoda("p","test_p")
samplesCoda("rd","test_rd")
## process output with CODA
library(coda)
## Get the coda chains from WinBUGS output
coda.out1 <- read.coda("test_pCODAchain1.txt", "test_pCODAindex.txt")
HPDinterval(coda.out1,prob=0.95)

coda.out2 <- read.coda("test_rdCODAchain1.txt", "test_rdCODAindex.txt")
HPDinterval(coda.out2,prob=0.95)

#######
#output
> samplesStats("*")
# the summarized results

mean
sd
MC_error
val2.5pc
median val97.5pc start sample
p[1] 0.253600 0.004377 4.425e-05
0.245100 0.253600
0.26230
1001
10000
p[2] 0.252200 0.004338 4.379e-05
0.243700 0.252200
0.26070
1001
10000
p[3] 0.249000 0.004339 4.169e-05
0.240500 0.248900
0.25750
1001
10000
p[4] 0.245200 0.004324 4.284e-05
0.236700 0.245200
0.25360
1001
10000
rd
0.004679 0.007138 7.070e-05 -0.009277 0.004691
0.01851
1001
10000
> HPDinterval(coda.out1,prob=0.95)

      lower upper
p[1] 0.2453 0.2624
p[2] 0.2439 0.2607
p[3] 0.2405 0.2574
p[4] 0.2368 0.2536
attr(,"Probability")
[1] 0.95

> HPDinterval(coda.out2,prob=0.95)

      lower upper
rd -0.01095 0.01688
attr(,"Probability")
[1] 0.95

testdata.txt (52 Bytes)

2010-3-1 17:20:10 上传

testmodel.txt (189 Bytes)

2010-3-1 17:20:10 上传

奇怪了，我把testdata.txt和testmodel.txt两文件放在my documents下了，可R运行时却总提示这两个文件不存在

那就请把两个文件放在
C:\Program Files\R\R-2.9.2\library\BRugs\OpenBUGS\Examples
然后程序的最前面,加入这两行
## Now setting the working directory to the examples' one:
oldwd <- getwd()
setwd(system.file("OpenBUGS", "Examples", package="BRugs"))

程序的最后面,加入这一行
## switch back to the previous working directory:
setwd(oldwd)