-
Notifications
You must be signed in to change notification settings - Fork 0
/
GMH_thin5.R
49 lines (32 loc) · 1.13 KB
/
GMH_thin5.R
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
Gibbs_MH_xMatrix=read.csv(file="Gibbs_MH_xMatrix.csv")
Gibbs_MH_yMatrix=read.csv(file="Gibbs_MH_yMatrix.csv")
dim(Gibbs_MH_yMatrix)
x=Gibbs_MH_xMatrix[,2:101]
y=Gibbs_MH_yMatrix[,2:101]
chain=cbind(x[,1],y[,1])
colnames(chain)=c("x","y")
par(mfrow=c(2,1))
plotAuto(chain)
x_thin5=x[seq(from=1,to=99000,by=5),]
y_thin5=y[seq(from=1,to=99000,by=5),]
dim(x_thin5)
corrs_GMH_thin5=rep(NA,100)
for(i in 1:100){
corrs_GMH_thin5[i]=cor(x_thin5[,i],y_thin5[,i])
}
mean(corrs_GMH_thin5)
sd(corrs_GMH_thin5)
LCI_GMH_thin5=mean(corrs_GMH_thin5)-qnorm(0.975,mean=0,sd=1)*sd(corrs_GMH_thin5)
UCI_GMH_thin5=mean(corrs_GMH_thin5)+qnorm(0.975,mean=0,sd=1)*sd(corrs_GMH_thin5)
par(mfrow=c(1,1))
plot(corrs_GMH_thin5,xlab="MEAN=-0.06070848 SD=0.005857476")
abline(h=mean(corrs_GMH_thin5),col="red")
abline(h=LCI_GMH_thin5,col="blue")
abline(h=UCI_GMH_thin5,col="blue")
chain_thin5=cbind(x[,1],y[,1])
colnames(chain_thin5)=c("x","y")
plotAuto(chain)
plotAuto(chain, thin=5)
plotAuto(chain_thin5)
write.csv(x_thin5,file="Gibbs_MH_thin5_xMatrix.csv")
write.csv(y_thin5,file="Gibbs_MH_thin5_yMatrix.csv")