另外，你是否知道如何计算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)

 马上下载

testmodel.txt
大小:(189 Bytes)

 马上下载

奇怪了，我把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)

我试试，另外beg究竟该取多少合适，比如iteration是10000次，从history上也看不出从哪里开始burn-in，是不是beg取1000或2000，这个对最后计算rd的后验信息有无实质影响呢？

如果你在winbugs是如此设置
1.sample monitor tool
  beg:1001  (就是 burn in)
2.update tool
  updates:11000

按"stats" button输出的结果与
底下的R code是相同的. start 1001,  sample 10000

modelUpdate(1000)         # burn in
modelUpdate(10000)       # 10000 more iterations ...
samplesStats("*")            # the summarized results

> 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

若模型简单,要估的参数少,burn in设1000,2000就够了.
详细可参阅 manual14.pdf page-53 Checking convergence

也可以由winbugs直接输出文件,
供"CODA"继续分析HPDinterval.
这里我以 p node说明:
sample monitor tool 按"coda"button
会跳出两个窗口:(1)CODA index (2)CODA for chain 1
请分别存成pCODAindex.txt, pCODAchain1.txt
放在R可读到的working directory.
若不确定R目前的working directory
可键入getwd()
知道R目前的working directory后,
将pCODAindex.txt, pCODAchain1.txt放入.

执行下列程序:
library(coda)
## Get the coda chains from WinBUGS output
coda.out1 <- read.coda("pCODAchain1.txt", "pCODAindex.txt")
HPDinterval(coda.out1,prob=0.95)

另补充:
有了这些数据,HPDinterval也可自行算出:
library(coda)
coda.out1 <- read.coda("pCODAchain1.txt", "pCODAindex.txt")
M = 10000
prob = 0.95
chain1 = matrix(scan("pCODAchain1.txt"),4*M,2,byrow=T)
chain1 = matrix(chain1[,2],M,4)
chain1=apply(chain1,2,sort)
nsamp <- nrow(chain1)
npar <- ncol(chain1 )
gap <- max(1, min(nsamp - 1, round(nsamp * prob)))
init <- 1:(nsamp - gap)
inds <- apply(chain1[init + gap, ,drop=FALSE] - chain1[init, ,drop=FALSE],

2, which.min)
ans <- cbind(chain1[cbind(inds, 1:npar)],

chain1 [cbind(inds + gap, 1:npar)])
dimnames(ans) <- list(colnames(ans), c("lower", "upper"))
attr(ans, "Probability") <- gap/nsamp
ans
> ans

      lower upper
[1,] 0.2450 0.2620
[2,] 0.2436 0.2609
[3,] 0.2403 0.2573
[4,] 0.2368 0.2537
attr(,"Probability")
[1] 0.95

够义气，谢谢！研究中

哦  天啦！

8# epoh
你好，我也是十分想求教winbugS这个软件的使用放个，我在做毕业论文，能告诉我你的qq号吗？感谢啊~~~~~~~···

我收藏了，努力学习啊！！