#### 第2章 ####
#### 2.1節 ####

# パッケージMCMCpackの読み込み
library(MCMCpack)
# 図2.1
# dinvgamma(n, alpha, beta)
plot(dinvgamma(1:100,0.1,5),type="l", xlab="x",ylab=expression(InvGamma(a,b)),lwd=2,cex.axis=1.3,cex.lab=1.3)

# 図2.2
x <- seq(0,1,length=100)
plot(x,dbeta(x,0.5,0.5),type="l", xlab="x",ylab=expression(Beta(a,b)),lwd=2,cex.axis=1.3,cex.lab=1.3)

#### 2.2節 ####
#### 2.2.1 ####
# swissデータの読み込み
data(swiss)
summary(swiss)
# lm(被説明変数~説明変数, data=データ名)　
# .はその他のデータを省略することを意味する
summary(lm(Fertility~.,data=swiss)) 

kk<- ncol(swiss)
y<-swiss[,1]
X<-cbind(1,as.matrix(swiss[,2:kk]))
n<-nrow(X)
k<-ncol(X)

# βの不偏推定量
# solve() は逆行列を計算する
betahat <- solve(t(X)%*%(X))%*%t(X)%*%y
betahat
# σ^2の不偏推定量
S2=t(y-X%*%betahat)%*%(y-X%*%betahat)
sig2hat <- S2/(n-k)
sig2hat
# βの不偏分散
diag(as.real(sig2hat)*solve(t(X)%*%X))

#### 2.2.2 ####
# k x k 正定値対象行列
M <- diag(k)

# 事前情報を適当に与えた場合の事後平均と事後分散の計算
a <- 2.1
b <- 2
c <- 0.1
M=(1/c)*diag(k)
T=solve(solve(M)+solve(t(X)%*%X))
# βの事後平均
betavar <- solve(M+t(X)%*%X)%*%t(X)%*%y
betavar
# σ^2の事後平均
S2=t(y-X%*%betahat)%*%(y-X%*%betahat)
S2var <- (2*b+S2+t(betahat)%*%T%*%betahat)/(n+2*a-2)
S2var
# βの事後分散
betasig <- diag(as.real(2*b+S2+t(betahat)%*%T%*%betahat)/(n+2*a-2)*solve(M+t(X)%*%X))
betasig

#### 2.2.3 ####
#上限値（定数項を加えた説明変数の数mのもとで，自由度はn-m-1）
m=k
HPDUP<- betavar+qt(0.975,n-m-1)*sqrt(diag(as.real(S2var)*solve(t(X)%*%X)))
HPDUP
# 下限値
HPDLW<- betavar+qt(0.025,n-m-1)*sqrt(diag(as.real(S2var)*solve(t(X)%*%X)))
HPDLW

# 図2.3
bv <- seq(-0.3,0.5,length=100)
plot(bv,dnorm(bv,betavar[2],sqrt(betasig[2])),type="l",lwd=2,xlab="posterior of beta",ylab="density")
rect(HPDLW[2],-0.5,HPDUP[2],dnorm(HPDUP[2],betavar[2],sqrt(betasig[2])),col="grey",border=NA)

#### 2.2.4
# Jeffreys' prior
# betahat
# S2/(n-k-3)

#### 2.2.5
# ZELLNERs' G-prior
# βの事後平均
(c/(c+1))*betahat
# σ^2の事後平均
(S2+t(betahat)%*%t(X)%*%X%*%betahat/(c+1))/(n-2)
# βの事後分散
diag(c/n*(c+1)*as.real(S2+t(betahat)%*%t(X)%*%X%*%betahat/(c+1))*solve(t(X)%*%X))

#### 2.3節 ####
nt <- 15
yt <- y[1:nt]
yp <- y[(nt+1):n]
Xt <- X[1:nt,]
Xp <- X[(nt+1):n,]
summary(lm(yt~Xt-1))

ndraw <- 1000
bt <- solve(t(Xt)%*%(Xt))%*%t(Xt)%*%yt
s2t<-t(yt-Xt%*%bt)%*%(yt-Xt%*%bt)
s2t <- s2t/(nt-k)
tau <- 1/sqrt(s2t)
bv <- t(matrix(rnorm(ndraw*k,bt,1/sqrt(tau*sum(Xt^2))),nrow=k))
ypred<-matrix(0,ncol=n-nt,nrow=ndraw)
ypredvar <- rep(0,(n-nt))
for(i in 1:ndraw){
	ypred[i,] <- rnorm(n-nt,bv[i,]%*%t(Xp),1/sqrt(tau))}
for(i in 1:(n-nt)){ypredvar[i] <- mean(ypred[,i])}
ypredvar

# 2.3.3 事後予測分布
mu <- rnorm(ndraw,mean(y),1/sqrt(n))
yrep <- matrix(0,nrow=ndraw,ncol=n)
for(i in 1:ndraw){yrep[i,] <- rnorm(n, mu[i],1)}
yrep

#### 2.4節 ####
ndraw <- 1000	 # 繰り返し計算回数
#　τの事前情報を計算
tauval <- rgamma(ndraw,(n-k)/2,(n-k)*sig2hat/2)
bval <- matrix(0,nrow=ndraw,ncol=k)
# 多変量正規分布の関数mvrmormを使うためライブラリMASSを呼び出す
library(MASS)
for(i in 1:ndraw){
V <- (1/tauval[i])*solve(t(X)%*%(X))
bval[i,] <-mvrnorm(1,betahat, V)}
# パラメータの事後平均・分散・標準偏差を計算
bval.mean <- matrix(0,nrow=k,ncol=1)
for(i in 1:k){bval.mean[i,] <- mean(bval[,i])}
bval.mean
bval.var <- matrix(0,nrow=k,ncol=1)
for(i in 1:k){bval.var[i,] <- var(bval[,i])}
bval.var
bval.sd<-sqrt(bval.var)

# 計算結果をプロットする
par(mfrow=c(3,1))
# 3番目の変数（Examination）のパラメータ事後確率密度分布をプロットする
plot(density(bval[,3]),main=" ",ylab="確率密度",cex.axis=1.3,cex.lab=1.3)
# 3番目のパラメータのndraw回の計算結果をプロットする
plot(1:ndraw,bval[,3],type="l",ylab="パラメータ",xlab="計算回数",cex.axis=1.3,cex.lab=1.3)
# 事後確率密度分布の自己相関を計算する
bval.auto <- matrix(0,nrow=(ndraw-1),ncol=k)
for(i in 1:(ndraw-1)){bval.auto[i,] <- bval[i,]/bval[i+1,]}
# 3番目のパラメータの自己相関をプロットする
plot(1:(ndraw-1),bval.auto[,3],type="l",ylab="自己相関",
xlab="計算回数",cex.axis=1.3,cex.lab=1.3)

# 事後情報を要約する
s2=t(y-X%*%bval.mean)%*%(y-X%*%bval.mean)
s2hat <- S2/(n-k)
s2hat
# 95%信頼区間を計算する
SD <- matrix(0,ncol=1,nrow=k)
XXI <- solve(t(X)%*%(X))
for(j in 1:k){SD[j,] <- sqrt(s2hat*XXI[j,j])}
# 上限値
HPDUP <- bval.mean+qt(0.975,n-k)*SD
HPDUP
# 下限値
HPDLW <- bval.mean+qt(0.025,n-k)*SD
HPDLW

#### 2.5節 ####
c <- c0 <- 0.1; q <- 1
# 各パラメータのベイズ・ファクター
log10BF10 <- matrix(0,ncol=1,nrow=k)
for(i in 1:k){
X0 <- X[,-c(i)]
P0=X0%*%solve(t(X0)%*%X0)%*%t(X0)
marg.dist0=(c0+1)^(-(k+1-q/2))*(t(y)%*%y-c0/(c0+1)*t(y)%*%P0%*%y)^(-n/2)
P=X%*%solve(t(X)%*%X)%*%t(X)
marg.dist=(c+1)^(-k/2)*(t(y)%*%y-c/(c+1)*t(y)%*%P%*%y)^(-n/2)
log10BF10[i,] <- log10(marg.dist/marg.dist0)}
log10BF10

# モデル全体のベイズ・ファクター
q <- 2
X0 <- X[,-c(5,6)]
marg.dist0=X0%*%solve(t(X0)%*%X0)%*%t(X0)
lulu0=(c0+1)^(-(k+1-q/2))*(t(y)%*%y-c0/(c0+1)*t(y)%*%P0%*%y)^(-n/2)
P=X%*%solve(t(X)%*%X)%*%t(X)
marg.dist=(c+1)^(-k/2)*(t(y)%*%y-c/(c+1)*t(y)%*%P%*%y)^(-n/2)
log10(marg.dist/lulu0)

# BIC
# 定数項と最初の4つの説明変数を抽出
k1 <- 5
X1 <- X[,1:k1]
# パラメータの不偏推定量と誤差の二乗和を計算
betahat1 <- solve(t(X1)%*%(X1))%*%t(X1)%*%y
S12=t(y-X1%*%betahat1)%*%(y-X1%*%betahat1)
# BICを計算
BIClogM10 <- (n/2)*(log10(S12)-log10(S2))+(k1-k)/2*log10(n)
BIClogM10

