求助统计学高手！ [推广有奖]

11楼

epoh 发表于 2011-10-1 14:39:57

我以另一个例子来解释:

以winbugs产生的物件schools.sim

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

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

请将其放进c:\Bugs\school

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

若有不同,请自行更改

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

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

school.bug , schools.dat

school.rar (432 Bytes)

# 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).