#### 第3章 ####
#### 3.1節 ####

# 時系列回帰モデルを作成
T <- 500
et <- rnorm(T,0,1)
ro <- 0.5
xt <- matrix(0,nrow=1,ncol=T)
xt[1,1] <- et[1]
for(i in 2:T){xt[,i]<-xt[,i-1]*ro+et[i]}
xt.mean <- matrix(0,nrow=1,ncol=T)
xt.mean[1]<-xt[1]
for(i in 2:T){xt.mean[,i]<-sum(xt[,1:i])/i}
par(mfrow=c(3,1))
plot(1:T, xt, type="l", xlab="T", ylab="Xt")
plot(xt[1:T-1], xt[2:T], xlab="Xt", ylab="Xt+1")
plot(1:T, xt.mean,type="l")

# 標本平均の収束過程
T <- 500
et <- matrix(0,nrow=5,ncol=T)
for(i in 1:5){et[1,]<-rnorm(T,0,1)}
ro <- 0.5
xt <- matrix(0,nrow=5,ncol=T)
xt[,1] <- runif(5,-10,10)
for(j in 1:5){
for(i in 2:T){xt[j,i]<-xt[j,i-1]*ro+et[i]}}
xt.mean <- matrix(0,nrow=5,ncol=T)
xt.mean[1]<-xt[1]
for(j in 1:5){
for(i in 2:T){xt.mean[j,i]<-sum(xt[j,1:i])/i}}
plot(xt.mean[1,],type="l", ylim=c(min(xt.mean),max(xt.mean)), xlim=c(1,T), lwd=2)
lines(xt.mean[2,],col="blue", lwd=2)
lines(xt.mean[3,],col="green", lwd=2)
lines(xt.mean[4,],col="red", lwd=2)
lines(xt.mean[5,],col="cyan", lwd=2)

#### 3.2節 ####
# 初期値(10, 10)の場合
T<-500; ro<-0.7; y0<-10; mu1<-mu2<-0; sig1<-sig2<-1
y1 <- rep(y0,T); y2<-rep(y0,T)
for(i in 2:T){
y2[i] <- rnorm(1, (mu2+ro*sig2/sig1*(y1[i-1]-mu1)), sqrt(sig2*(1-ro^2)))
y1[i] <- rnorm(1, (mu1+ro*sig1/sig2*(y2[i-1]-mu1)), sqrt(sig1*(1-ro^2)))}
# 初期値(-10, 10)の場合
T<-500; ro<-0.7; y00<--10; mu1<-mu2<-0; sig1<-sig2<-1
y3 <- rep(y0,n); y4<-rep(y00,n)
for(i in 2:T){
y4[i] <- rnorm(1, (mu2+ro*sig2/sig1*(y3[i-1]-mu1)), sqrt(sig2*(1-ro^2)))
y3[i] <- rnorm(1, (mu1+ro*sig1/sig2*(y4[i-1]-mu1)), sqrt(sig1*(1-ro^2)))}
# 図3.3
# par(mfrow=c(1,2))
# 図3.3 (a)
plot(y1,y2,type="o",xlim=c(min(y1)-2,max(y1)+2),ylim=c(min(y2)-2,max(y2)+2))
# 図3.3 (b)
plot(y3,y4,type="o",xlim=c(min(y3)-2,max(y3)+2),ylim=c(min(y4)-2,max(y4)+2))

# 初期値(10, 10)の場合の信頼区間
HPDUPy1<-mean(y1)+qt(0.975,T-2)*sd(y1)
HPDLWy1<-mean(y1)+qt(0.025,T-2)*sd(y1)
HPDUPy1
HPDLWy1
HPDUPy2<-mean(y2)+qt(0.975,T-2)*sd(y2)
HPDLWy2<-mean(y2)+qt(0.025,T-2)*sd(y2)
HPDUPy2
HPDLWy2

#### 3.5節 ####
# 線形回帰モデルへの適用
library(MCMCpack)
data(swiss)
n<-nrow(swiss)
k<-ncol(swiss)
# 通常最小二乗法の適用
sample.lm <- lm(Fertility~.,data=swiss)
summary(sample.lm)
# Gibbs sampling
# burnin=稼働検査期間，mcmc=MCMCの生成回数，b0とB0は事前分布の平均と分散
# c0とd0はτに関するガンマ関数の平均と分散に関するパラメータ
sample.gibbs.post1 <- 
MCMCregress(Fertility~.,data=swiss,
burnin=10000, mcmc=100000,
marginal.likelihood="Chib95", 
b0=0, B0=0.001,c0=0.001, d0=0.001)
summary(sample.gibbs.post1, digit=3)
plot(sample.gibbs.post1)

# Gelman-Rubin統計量
sample.gibbs.post2 <- 
MCMCregress(Fertility~.,data=swiss,burnin=10000, mcmc=100000,
marginal.likelihood="Chib95", b0=1, B0=0.05,c0=0.1, d0=0.1)
sample.gibbs.post.all<-mcmc.list(sample.gibbs.post1,sample.gibbs.post2)
gelman.diag(sample.gibbs.post.all)
# Gewekeの方法
geweke.diag(sample.gibbs.post1)
geweke.diag(sample.gibbs.post2)
# Raftery-Lewisの方法
raftery.diag(sample.gibbs.post1)

# boaパッケージを使う
#library(boa)
# codaパッケージを使う
library(coda)
codamenu()

# ベイズファクターの計算
# model1:最初に推定したモデル
sample.gibbs.post1 <- 
MCMCregress(Fertility~.,data=swiss,
marginal.likelihood="Chib95",
burnin=10000, mcmc=100000,
b0=0, B0=0.001,c0=0.001, d0=0.001)
summary(sample.gibbs.post1, digit=3)

# model3:AgricultureとExaminationを除いて推定
sample.gibbs.post3<-
MCMCregress(Fertility~Education+Catholic+Infant.Mortality,data=swiss,
##MCMCregress(Fertility~Education+Catholic+Examination,data=swiss,
marginal.likelihood="Chib95",
burnin=10000, mcmc=100000,
b0=0, B0=0.0015,c0=0.001, d0=0.001)

# ベイズファクター
BayesFactor(sample.gibbs.post1,sample.gibbs.post3)

# WinBUGSによる計算
# BUGSコード（model1.txt）# Rのディレクトリに保存
model{
for(i in 1:n){
y[i]~dnorm(mu[i],tau)
mu[i] <-beta0 + beta1*X[i,1]+ beta2*X[i,2]+ beta3*X[i,3]+ beta4*X[i,4]+ beta5*X[i,5]}
beta0 ~dnorm(0,0.0001)
beta1 ~dnorm(0,0.0001)
beta2 ~dnorm(0,0.0001)
beta3 ~dnorm(0,0.0001)
beta4 ~dnorm(0,0.0001)
beta5 ~dnorm(0,0.0001)
tau ~dgamma(0.001, 0.001)
}

# R上でWinBUGSを動かす
library(R2WinBUGS)
data(swiss)
y<-swiss[,1]
k<-ncol(swiss)
X<-as.matrix(swiss[,2:k])
n<-nrow(X)
#dim(X)<-c(n,k-1)
data <- list("n", "y", "X")
in1 = list(beta0=0, beta1=0, beta2=0, beta3=0, beta4=0, beta5=0, tau=1)
in2 = list(beta0=1, beta1=1, beta2=1, beta3=1, beta4=1, beta5=1, tau=1)
in3 = list(beta0=2, beta1=2, beta2=2, beta3=2, beta4=2, beta5=2, tau=1)
inits = list(in1, in2, in3)
# inits = function() {list(beta0 = 0, beta1 = 0, beta2 = 0, beta3 = 0, beta4 = 0, beta5 = 0, tau = 1)}
parameters <- c("beta0", "beta1","beta2","beta3","beta4","beta5","tau")

sample.wb <- bugs(data, inits, parameters,
model.file="C:/Program Files/R/R-2.7.2/model1.txt",
debug=FALSE,
n.chains=3, n.iter=100000, n.burnin=10000, n.thin=300,
codaPkg=TRUE,
bugs.directory="C:/Program Files/WinBUGS14",
program="WinBUGS",
working.directory="C:/Program Files/R/R-2.7.2")

print(sample.wb,digits=3)
plot(sample.wb)

# WinBUGSを使うと、うまく計算できないことが少なくありません。
# テキストでは　model.file="model1.txt"　としているが、R上で、
# 以下にエラー bugs(data, inits, parameters, model.file = "model1.txt", debug = FALSE,  : 
# model1.txt does not exist.
# と表示された場合には、上述のようにbugs()中で
# model.file="C:/Program Files/R/R-2.7.2/model1.txt" 
# とファイルのパスを指定するとよい。

# テキストでは　working.directory=NULL　としているが、
# WinBUGS上で、
# cannot open C:/ ..... /data.txt
# と表示されたら、
# working.directory="C:/Program Files/R/R-2.7.2"
# とディレクトリを指定するとよい。

# 初期値を指定せずに
# sample.wb <- bugs(data, inits=NULL, ...
# としてもよい

# WinBUGS上で
# shape parameter (r) of gamma tau too small -- cannot sample
# と表示されたら、model1.txtのtau ~dgamma(　)中の値を大きくするとよい。

# 上記の方法でうまくいかない場合、以下の方法も試してみてください。
y<- swiss[,1]
x1<-swiss[,2]
x2<-swiss[,3]
x3<-swiss[,4]
x4<-swiss[,5]
x5<-swiss[,6]
x <- matrix(c(x1,x2,x3,x4,x5),ncol=5,nrow=47)
n=47
data <- list ("n", "y", "x") 

in1 <- list(beta0=0, beta1=0, beta2=0, beta3=0, beta4=0, beta5=0, tau=1)
in2 <- list(beta0=1, beta1=1, beta2=1, beta3=1, beta4=1, beta5=1, tau=1)
in3 <- list(beta0=2, beta1=2, beta2=2, beta2=2, beta4=2, beta5=2, tau=1)
inits <- list(in1,in2,in3)
parameters <- c("beta0", "beta1", "beta2", "beta3", "beta4", "beta5", "tau")

sample.wb <- bugs(data, inits, parameters, 
model.file="C:/Program Files/R/R-2.7.2/model.txt",
n.chains = 3, n.iter = 1000, 
debug=TRUE,
bugs.directory = "c:/Program Files/WinBUGS14/",
working.directory = "C:/Program Files/R/R-2.7.2/")

print(sample.wb,digits=3)
plot(sample.wb)
