###4章 ベイズ流の柔軟な仮説の検証


###4.2 選挙の事前調査
###ディリクレ分布の乱数発生関数定義
rDir <- function(R,alpha) {#引数の型の検査は行っていないので注意！
  m <- length(alpha)#カテゴリー数
  g_mat <- t(matrix(rgamma(m*R,alpha,1),m,R))#ガンマ分布から乱数を作成
  return(g_mat/rowSums(g_mat))#R×カテゴリー数の行列形式の結果を返す
}
#100000回乱数を発生させる
P_sim <- rDir(R=100000,alpha=c(231+1,219+1,214+1,336+1))
#無投票を除き各仮説に合致しているかを判定
tb <- apply(P_sim[,-4],1,which.max)
#集計
out <- table(tb)/length(tb)
#ラベルを付ける
names(out) <- c("A氏が当選","B氏が当選","C氏が当選")
#結果を表示
print(out)


###4.3 ウェイト集計の仮説検定
#回数の指定
R <- 100000
#事後分布のシミュレート
p1_sim <- rbeta(R,60+1,40+1)#p1
p2_sim <- rbeta(R,40+1,60+1)#p2
p_sim <- 0.2*p1_sim+0.8*p2_sim#p
#仮説が合致した率を求める
mean(p_sim>0.5)


###4.4 重回帰モデルの偏回帰係数の比較
#データの取り込み
#現在のディレクトリに「満足度データ.csv」があるとする
Data <- read.csv("満足度データ.csv",row.names=1)

#目的変数
y <- Data[["総合満足度"]]
#説明変数
X <- Data[,c("品質満足度","魅力満足度","接客満足度","アフターサービス満足度")]

###重回帰モデルの最尤推定
#最尤法の推定結果
out_ML <- lm(y~.,data=X)
#結果のサマリー
summary(out_ML)

###重回帰モデルのベイズ推定
n <- nrow(X)#サンプルサイズ
m <- ncol(X)#説明変数の数
#事前分布のハイパーパラメータの設定
#sigma2(逆カイ二乗分布)
nu0 <- 3
c0 <- 2
#偏回帰係数(多変量正規分布)※切片含む
beta0 <- rep(0,m+1)
V0 <- diag(100^2,m+1)

#事後分布ための計算
V0_inv <- solve(V0)
X1 <- cbind(1,as.matrix(X)); XX1 <- crossprod(X1)
H1_inv <- XX1+V0_inv; H1 <- solve(H1_inv)
h2 <- crossprod(X1,y)+crossprod(V0_inv,beta0)
beta_hat <- H1%*%h2#betaの事後平均
n1 <- n+nu0#自由度
n1c1 <- as.vector(nu0*c0+crossprod(y)+crossprod(beta0,crossprod(V0_inv,beta0))-
                    t(beta_hat)%*%H1_inv%*%beta_hat)

#事後分布の要約
sd_beta <- sqrt(diag(solve((H1_inv)/(n1c1/n1))*(n1)/(n1-2)))
out_Bayes <- cbind(beta_hat,sd_beta,sd_beta*(qt(0.025,n1))+beta_hat,sd_beta*(qt(0.975,n1))+beta_hat)

#乱数で仮説検定
R <- 100000
sigma2 <- 1/rgamma(R,n1/2,n1c1/2)#逆カイ二乗分布の乱数発生
beta2 <- matrix(0,R,m+1)#偏回帰係数の乱数を入れる行列
for(r in 1:R) {
  beta2[r,] <- t(chol(sigma2[r]*H1))%*%rnorm(m+1)+beta_hat
}
#仮説「すべての偏回帰係数が正である」
mean(apply(beta2[,-1]>0,1,prod))
#仮説「一番影響度が高い説明変数」
res <- table(apply(beta2[,-1],1,which.max))/R
names(res) <- colnames(X)
res