高手快来指点迷津啊

根据你给的文献jd empi.pdf

page 18/47,

作者是用matlab实证

burn=2000,iterations=3000,

为了要画page 19/47 Jump Probability图形,这也是恰巧合理.

因为Xiy matrix 3849 x 3000,就已近 300 Mb,

若iterations增加到10000,

一般软件无法读入及处理这么大的文件

R就无法处理,底下图形我是在matlab画的.

package "R2WinBUGS"顾名思义

就是R & Winbugs 的桥梁

通过R及package "coda"可以进一步检验模型及画图

如:gelman.diag() Gelman and Rubin＇s convergence diagnostic

library(coda)
out=as.mcmc.list(JD1.sim)
gelman.diag(out)

Potential scale reduction factors:

         Point est. Upper C.I.
deviance       1.04       1.12
lambday       1.04       1.12
mu              1.00       1.00
muy            1.00       1.00
sig              1.02       1.06
sigy            1.07       1.27
tau             1.02       1.05
tauq           1.06       1.23
########gelman.plot

gelman.plot(out)

########chain=2,burn in=2000,iter=5,000 result
########Node statistics
node      mean      sd         MC error  2.5%    median   97.5%  start sample
deviance8192.0   109.6        6.519    7971.0 8195.0    8401.0   2001 6000
lambday 0.05502 0.0084615.251E-4 0.03947 0.05464 0.07272 2001 6000
mu        0.07177 0.01326 2.825E-4  0.04573 0.07179 0.09771 2001 6000
muy     -0.4368   0.2444  0.01728  -0.9498  -0.4286   0.0098152001 6000
sig        0.7349   0.01331 6.43E-4   0.7098   0.7348    0.761     2001 6000
sigy      2.883     0.2375   0.01903   2.47      2.873      3.376     2001 6000
tau       1.854     0.06707 0.003248 1.727     1.852      1.985     2001 6000
tauq     0.06615 0.0099737.981E-4 0.04893  0.06544  0.08696 2001 6000

log.odc

epoh老师，您好！
在R2winbugs包打开的情况下，运行您的程序，为什么会出现以下错误？如何更正和运行？

> library(coda)
> out=as.mcmc.list(JD1.sim)
错误于as.mcmc.list(JD1.sim) : 找不到对象'JD1.sim'
> gelman.diag(out)
错误于as.mcmc.list(x) : 找不到对象'out'
> gelman.plot(out)

还有，您提供的：log.odc文件如何运行？
可以在winbugs14中打开log.odc文件，但是无法运行。
display(log)
check(c:/Bugs/JD1//JD1.bug.txt)
model is syntactically correct
data(c:/Bugs/JD1/data.txt)
data loaded
compile(2)
model compiled
inits(1,c:/Bugs/JD1/inits1.txt)
this chain contains uninitialized variables
inits(2,c:/Bugs/JD1/inits2.txt)
this chain contains uninitialized variables
gen.inits()
initial values generated, model initialized
thin.updater(1)
update(2000)
set(mu)
set(muy)
set(tau)
set(tauq)
set(sig)
set(sigy)
set(lambday)
set(deviance)
dic.set()
command #Bugs:dic.set cannot be executed (is greyed out)
update(3000)
coda(*,c:/Bugs/JD1/coda)
stats(*)

Node statistics
         node         mean         sd         MC error        2.5%        median        97.5%        start        sample
        deviance        8192.0        109.6        6.519        7971.0        8195.0        8401.0        2001        6000
        lambday        0.05502        0.008461        5.251E-4        0.03947        0.05464        0.07272        2001        6000
        mu        0.07177        0.01326        2.825E-4        0.04573        0.07179        0.09771        2001        6000
        muy        -0.4368        0.2444        0.01728        -0.9498        -0.4286        0.009815        2001        6000
        sig        0.7349        0.01331        6.43E-4        0.7098        0.7348        0.761        2001        6000
        sigy        2.883        0.2375        0.01903        2.47        2.873        3.376        2001        6000
        tau        1.854        0.06707        0.003248        1.727        1.852        1.985        2001        6000
        tauq        0.06615        0.009973        7.981E-4        0.04893        0.06544        0.08696        2001        6000
dic.stats()

DIC
history(*,c:/Bugs/JD1/history.odc)

我以另一个例子来解释:

以winbugs产生的物件schools.sim

可以经由library(coda),作进一步检验分析.

这个例子需要两个文件school.bug , schools.dat

请将其放进c:\Bugs\school

也请注意你的winbugs安装的位置,

若有不同,请自行更改

所有结果都存在c:\Bugs\school,

包含你说的log.odc,那就是执行结果的文件

school.bug , schools.dat

# in R command window

library(R2WinBUGS)
schools <- read.table ("c:/Bugs/school/schools.dat", header=TRUE)
J <- nrow(schools)
y <- schools$estimate
sigma.y <- schools$sd
data <- list ("J", "y", "sigma.y")
inits <- function() {list (theta=rnorm(J,0,100), mu.theta=rnorm(1,0,100), sigma.theta=runif(1,0,100))}
parameters <- c("theta", "mu.theta", "sigma.theta")

schools.sim <- bugs(data, inits, parameters.to.save=parameters,"school.bug",n.chains=2,
               n.thin=1,n.iter=10000,n.burnin=5000,debug=TRUE,DIC=TRUE,
               bugs.directory="d:/WinBUGS14/",working.directory = "c:/Bugs/school")
print(schools.sim)

library(coda)
out=as.mcmc.list(schools.sim)
gelman.plot(out)

#########################################################
Inference for Bugs model at "school.bug", fit using WinBUGS,
2 chains, each with 10000 iterations (first 5000 discarded)
n.sims = 10000 iterations saved
            mean  sd  2.5%  25%  50%  75% 97.5% Rhat n.eff
theta[1]    11.9 8.3  -1.8  6.4 10.9 15.8  31.8  1.0 10000
theta[2]     8.1 6.4  -5.2  4.2  8.1 12.3  20.7  1.0   800
theta[3]     6.3 7.9 -12.1  2.2  6.9 11.3  20.9  1.0  1600
theta[4]     7.9 6.6  -5.6  3.9  7.9 12.2  21.1  1.0   780
theta[5]     5.3 6.4  -8.9  1.4  5.7  9.7  16.2  1.0   420
theta[6]     6.3 6.8  -9.0  2.5  6.8 10.8  18.6  1.0  2700
theta[7]    11.1 6.7  -0.8  6.5 10.6 15.0  26.3  1.0  5400
theta[8]     8.8 8.0  -7.3  4.3  8.6 13.1  26.1  1.0 10000
mu.theta     8.2 5.1  -2.0  5.0  8.2 11.5  18.3  1.0  2100
sigma.theta  6.9 5.7   0.2  2.8  5.5  9.5  21.2  1.1    94
deviance    60.2 2.2  56.7 58.8 59.9 61.2  65.6  1.0  3800

For each parameter, n.eff is a crude measure of effective sample size,
and Rhat is the potential scale reduction factor (at convergence, Rhat=1).

DIC info (using the rule, pD = Dbar-Dhat)
pD = 2.9 and DIC = 63.1
DIC is an estimate of expected predictive error (lower deviance is better).

epoh老师，您好！
以上R估计的图片，我想问一下：theta[1]、theta[2]、theta[3]这三张图片
跑到哪儿去了？为什么在R中没有显示出来？
非常感谢！

点击图片添加到帖子内容中

好难啊