1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
|
##################################################
# Section 6.8 Example of Output Analysis
##################################################
library(LearnBayes)
d=list(int.lo=c(-Inf,seq(66,74,by=2)),
int.hi=c(seq(66,74,by=2), Inf),
f=c(14,30,49,70,33,15))
library(coda)
library(lattice)
start=c(70,1)
fit=laplace(groupeddatapost,start,d)
start=c(65,1)
proposal=list(var=fit$var,scale=0.2)
bayesfit=rwmetrop(groupeddatapost,proposal,start,10000,d)
dimnames(bayesfit$par)[[2]]=c("mu","log sigma")
xyplot(mcmc(bayesfit$par[-c(1:2000),]),col="black")
S=readline(prompt="Type <Return> to continue : ")
windows()
par(mfrow=c(2,1))
autocorr.plot(mcmc(bayesfit$par[-c(1:2000),]),auto.layout=FALSE)
summary(mcmc(bayesfit$par[-c(1:2000),]))
batchSE(mcmc(bayesfit$par[-c(1:2000),]), batchSize=50)
S=readline(prompt="Type <Return> to continue : ")
start=c(70,1)
proposal=list(var=fit$var,scale=2.0)
bayesfit=rwmetrop(groupeddatapost,proposal,start,10000,d)
dimnames(bayesfit$par)[[2]]=c("mu","log sigma")
sim.parameters=mcmc(bayesfit$par[-c(1:2000),])
windows()
xyplot(mcmc(bayesfit$par[-c(1:2000),]),col="black")
s=readline(prompt="Type <Return> to continue : ")
windows()
par(mfrow=c(2,1))
autocorr.plot(sim.parameters,auto.layout=FALSE)
summary(sim.parameters)
batchSE(sim.parameters, batchSize=50)
|
