ベイズ統計データ分析
R&WinBUGS
古谷知之著
羽藤研究室
修士2年 松村草也

目次

1.

ベイズアプローチの基本

2.

ベイズ推論

3.

1.

事前分布

2.

線形回帰モデルに対する事前分布と事後分布

3.

予測分布

4.

線形回帰モデルにおける事後密度の生成

5.

モデル選択

マルコフ連鎖モンテカルロ法1.

マルコフ連鎖2.

ギブズ・サンプラー

3.

メトロポリスーヘイスティング法

4.

収束判定

5.

線形回帰モデルへのギブズ・サンプラーの適用ベイズアプローチの基本
Basis for Bayesian Approachベイズアプローチの基本
複数のモデルでデータを再現１つのモデルを求める••モデル（パラメータセット）によって説明出来るデータが異なってくる．

•

データが観測されればそれは固定されることになる．（θの確率分布には影響しない）

P(観測データ)は同時確率P(観測データ|モデル)P(モデル)の和であり，モデルを特定しない場合に全てのデータが観測される「周辺確率」である．

ベイズの定理↑モデルが与えられた状態でデータが観測される尤もらしさ．

↑モデルを特定しない場合に全てのデータが観測される確率

：事後分布
：尤度関数
：事前分布
：周辺確率
↑データが与えられた時にモデルが成立する確率

従来のモデル推定方法とベイズアプローチの違い• データを得る前のモデル（パラメタ）に関する事前の知識としての事前情報をあらかじめ用意し，実験や調査などから得られたデータが成立するかどうかという仮説の確からしさとを掛け合わせる．

• 頻度主義的立場からは，最も好ましいモデルを求める．

• ベイズ主義的立場からは，パラメータの尤もらしさの分布を求める．

ベイズ推論Baysian Estimation

事前分布の与え方とその特徴• ベイズ的手法を用いたデータ分析では，何らかの事前情報確率p(θ)を与える事により，事後確率密度分布を計算している．

• その際，事前情報については何ら情報を持たないままに，あてずっぽで事前情報確率を与えている．

• そこで，「与える事前情報確率を変えたらθの事後確率も変化するのか？」という疑問がわいてきます。。

事前分布の与え方とその特徴
事後情報■事前分布の例
自然共役事前分布
natural conjugate prior distribution

事後分布と同じ分布が得られる

非正則事前分布
improper prior distributionθに対する確率分布を標本分布に関して積分しても収束しない分布．

無情報事前分布noninformative prior distribution

事前情報を定数と置く事で事前情報からの依存を無くすことが出来る分布．

階層事前分布hierarchical prior distribution

未知パラメータθを標本iごとに推定するモデル．
ハイパーパラメータλを用いてθを推定する．

尤度関数 事前情報■共役事前分布の例
変数Xの従う分布(尤度関数)

事前分布

事後分布

正規分布

正規分布

正規分布

指数分布

ガンマ分布

ガンマ分布

ポアソン分布

ガンマ分布

ガンマ分布

線形回帰モデル - 定義
■線形回帰モデルの定義■計算例のデータについて
Rに組み込まれている，スイスの1888年の社会経済データ．
地域ごとの出生力/男性の農業従事者割合/教育水準/カソリッ
クの割合/乳児死亡率など

線形回帰モデル - 最小二乗法による解法■Rのlmによる計算例■尤度関数としてとらえてみると？

# swissデータの読み込み
data(swiss)summary(swiss)# lm(被説明変数~説明変数, data=データ名)# .はその他のデータを省略することを意味する
summary(lm(Fertility~.,data=swiss))■β，σに着目した計算→これはσ^2に関する逆ガンマ関数．

k<- ncol(swiss)y<-swiss[,1]X<-cbind(1,as.matrix(swiss[,2:k]))n<-nrow(X)k<-ncol(X)# βの不偏推定量# solve() は逆行列を計算する
betahat <- solve(t(X)%*%(X))%*%t(X)%*%ybetahat# σ^2の不偏推定量
S2=t(y-X%*%betahat)%*%(y-X%*%betahat)
sig2hat <- S2/(n-k)
sig2hat# βの不偏分散
diag(as.real(sig2hat)*solve(t(X)%*%X))

Mはk*kの正定値対象行列
JGは逆ガンマ関数

線形回帰モデル - 最小二乗法による解法
Call:lm(formula = Fertility ~ ., data = swiss)
Residuals:
Min
1Q Median
3QMax-15.2743 -5.2617 0.5032 4.1198 15.3213
Coefficients:
Estimate Std. Error t value Pr(>|t|)(Intercept)
66.91518 10.70604 6.250 1.91e-07 ***Agriculture-0.17211 0.07030 -2.448 0.01873 *Examination-0.25801 0.25388 -1.016 0.31546
Education-0.87094 0.18303 -4.758 2.43e-05 ***Catholic0.10412 0.03526 2.953 0.00519 **Infant.Mortality 1.07705 0.38172 2.822 0.00734 **--Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 7.165 on 41 degrees of freedom
Multiple R-squared: 0.7067,
Adjusted R-squared: 0.671
F-statistic: 19.76 on 5 and 41 DF, p-value: 5.594e-10•βとσの事前分布について何も情報を持たない状況で、適当にパラメータを設定して得られたもの。

線形回帰モデルのベイズ推定 - 尤度関数として扱う
β(ガウス型関数),σ(逆ガンマ関数) の事後平均β(ガウス型関数)の事後分散を計算し，比較を行う．

#### 2.2.2 ####M <- diag(k)a <- 2.1
b <- 2c <- 0.1
M=(1/c)*diag(k)T=solve(solve(M)+solve(t(X)%*%X))# βの事後平均
betavar <- solve(M+t(X)%*%X)%*%t(X)%*%ybetavar# σ^2の事後平均
S2=t(y-X%*%betahat)%*%(y-X%*%betahat)
S2var <- (2*b+S2+t(betahat)%*%T%*%betahat)/(n+2*a2)
S2var# βの事後分散
betasig <diag(as.real(2*b+S2+t(betahat)%*%T%*%betahat)/(n+2*a-2)*solve(M+t(X)%*%X))betasigについては情報がないので適当に設定．

→ 事前情報への依存度が高い．

信頼区間と最高事後密度
ところで，得られた事後分布はどの程度信頼出来る結果なのか？
■信頼区間

この区間内にθが入る確率が1-αであるため，100(1-a)％信頼区間と呼ぶ．
■最高事後密度(HPD:Highest Posterior Density)t検定量を定義．T_{n-k-1}は標本数n，自由度n-k-1のt分布である．ωiiはXTXの(i,i)要素を意味する．

経験ベイズによる推定

（１）
（２）

•
•

誤差項の分散σ^2を精度τ＝1/σ^2に置き換える。

•

モデルパラメータを求めるためには、βとτの事前分布が必要である。単純な例
として（２）のような非正則一様分布を与える。

τのように、事前分布の形状を決定するようなパラメタのことを超パラメタと呼ぶ。

経験ベイズによる推定• 線形回帰モデルのパラメータを，予測分布をシミュレーションする方法を用いてベイズ推定する．（この方法はMCMCの基礎となるものである．）• 同じ線形回帰モデルを使用．
#### 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.varbval.sd<-sqrt(bval.var)

τ→ε→βを繰り返し求める．
βとτの事前分布はどうするか？
→簡単な下記の非正則一様分布を与える．

経験ベイズによる推定

swissデータの線形回帰モデルのベイズ推定結果
平均

標準偏差

97.5%

2.5%

定数

67.2

1.29E+02

88.82945.586

農業-0.172

5.32E-03-0.030-0.314

出⽣⼒-0.246

6.68E-02

0.266-0.759教育-0.883

3.26E-02-0.513-1.252

カソリック

0.105

1.19E-03

0.177

0.034

幼児死亡率

1.059

1.52E-01

1.83

0.288

・「農業」はパラメタが負で有る可能性が高い．
・「出生力」はどちらかわからない．
→推定精度が低い．最小二乗法におけるパラメタ推定
結果ではt値が5%水準に達しておらず．

マルコフ連鎖モンテカルロ法
Markov Chain Monte Carlo method

はじめに• ベイズ推論では，与えられた事前情報に対して事後分布をサンプリングするシミュレーション計算により，予想事後分布を得る事が出来る．

• シミュレーションの計算回数が少ないと事前情報に強く依存する．

• 「何回程度計算すれば事前情報の影響を除去出来るか？」• 「膨大なシミュレーション計算を効率的に実行するには？」マルコフ連鎖について■推移確率(transition probability)■推移核(transition kernel)■時系列回帰モデル（マルコフ連鎖の代表例）

時系列回帰モデル - 数値計算結果•••

５通りの乱数を生成し、標本平均をプロットしている．

•

定常状態に至るまでの繰り返し計算を棄却し，定常状態にある推定値を用いて分析結果を議論した方が良いと考えられる．

•繰り返し計算を始めてからサンプリングを棄却するまでの期間を「稼働検査期間」(burn-in-period)という．

このような単純なマルコフ連鎖に限らず，MCMCでは事後分布の標本経路は初期値に依存する．
しかし，初期値がどのような値でも，計算回数を重ねるごとに標本平均の値が定常状態に近づく事が図からみてとれる．

ギブズ・サンプラー
マルコフ連鎖の推移核を構築するより一般的な方法．
■ギブズ・サンプラーのアルゴリズム
1.ドロー回数をsとする．s=(0,...,ndraw)
2.s=0に対して，初期値を決める．
3.xsからxs+1を生成する．

4.s<ndrawのとき，3に戻る．s=ndrawの時計算終了# 初期値(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.3par(mfrow=c(1,2))
plot(y1,y2,type="o",xlim=c(-10,10),ylim=c(-10,10))plot(y3,y4,type="o",xlim=c(-10,10),ylim=c(-10,10))ギブズ・サンプラー

初期値：(10,10)

初期値：(10,-10)↑初期値によらず定常状態に近づくことがわかる．

収束判定
MCMCでは，標本経路が初期値に依存しない定常状態であるかを見ることによって，収束判定を行う．

モデル間での差異・モデル内での差異
burn-in
ndraw

名称

メソッド名

概要

Gelman-Rubin統計量

gelman.diag()

モデル間の標本時系列の平均，モデル内の平均を利用した指標R

Gewekeの判定方法

geweke.diag()

平均値の差に関するZ検定を応用

Raftery-Lewisの診断方法

raftery.diag()

マルコフ連鎖の四分位偏差を用いる

離散選択モデルのサンプルプログラム■サンプルデータについて
目的変数：満足度2
説明変数：性別、年齢、年収のMNLモデルをMCMCにより、burnin1000回、ドロー10000回にて推定を行った結果を下記に示す。

Iterations = 1001:11000
Thinning interval = 1
Number of chains = 1
Sample size per chain = 10000
1. Empirical mean and standard deviation for each variable,
plus standard error of the mean:
Mean SD Naive SE Time-series SE
gend 0.45403 0.2286 0.002286
0.006974
age -0.13106 0.0886 0.000886
0.003121budget -0.06319 0.0921 0.000921
0.002944
2. Quantiles for each variable:
2.5% 25%50%
75% 97.5%
gend 0.002459 0.3033 0.45453 0.6001149 0.91065
age -0.306707 -0.1906 -0.13053 -0.0736864 0.04396
budget -0.243822 -0.1260 -0.06335 0.0005509 0.11695

推定結果