ベイズ及びMCMCの応用

ベイズ統計額が寄与している応用例を以下に簡単に紹介します。

ベイズ統計学はロジットモデル含む、一般化線形モデルのパラメータを確率分布に従うとしてその密度を推定することができる。モデル選択はAICみたいなBICという情報量もってして行う。(補足;AICの赤池さんの授業を聴講していたschwartsにより定義された)

階層ベイズとは事前分布の事前分布を考えるメタヒューリスティクな考え方。特にマーケティングの消費者間の異質性を考えたモデルの時に使われます。

項目反応理論(IRT)→ロジットみたいなもの。お勉強のテスト評価や、心理学のPOMS他あとで例題を解く。BMA:ベイジアンモデルアベレージ、モデルのアンサンブルに使われる。モデルのデータが与えられたときの事後分布をもって重み付けする。特にKDDcup(kaggleの前身)2009で単純ベイズモデルに対する変数選択で用いた論文あり。

BART:良く分からないが、回帰木のアンサンブルみたいなもの、gradient boosting のベイジアン的アナロジーか

参考文献

  • Marc BoulléCompression-Based Avergein of Selective Naive Bayes Classifiers』
  • 『ベイズモデリングによるマーケティング分析: 照井 伸彦』
  • Jim,Albet 『Bayesian Computetion with R』springer
  • 中妻照雄 (2007)『入門ベイズ統計学』朝倉書店
  • 伊庭幸人(2003)『ベイズ統計と統計物理』
  • 豊田秀樹(2008)統計ライブラリー『マルコフ連鎖モンテカルロ法』朝倉書店
  • C.M.ビショップ『パターン認識と機械学習 上』シュプリンガージャパン

ベイズ統計簡単入門

通常の頻度統計(正しい言い方ではないと思う)では最尤推定により点推定でパラメータを求めるがベイズはパラメータの確率分布を求めるという点で異なる。

例)ロジットを考える

データ集合:(Y,X)、Yはターゲット変数ベクトル、Xは説明変数行列

データ集合i行目のターゲットが1なる確率は

P(yi=1|xi)=1/(1+exp(-Z))

Z=α+β1*xi1+…+βn*xin

通常の場合→最尤法によりα、β1、βnを推定

パラメータの検定はp値とか用いる。

推定対象となる母パラメータは常にFIX(神のみが知る)されていて、95%信頼区間は100回区間つくると95回はそのfixされた値を含んでいるのでは?というニュアンス

母パラメータの区間推定は標準正規分布を仮定している。→だからいつも1.96とかどうのとかいう数字が出る。

Bayes→α、β1、…βnの確率分布を求める。

信頼区間はそのまま、∫(αの確率密度関数)dα=0.95となる区間

変な形の分布になる時もある。

事前分布、事後分布

上記確率分布はどのように求めるか?パラメータはある確率分布に従っていると仮定する

どんなのか分からないから、適当に決める。→事前分布(適当とはいい加減という意味ではない、ふさわしい事前分布)

データが与えられた時に、その適当に決めたパラメータの分布を逐次更新していく→更新後、事後分布と呼ぶ。

 

更新する方法。

 

・P(α|D):データが与えられたときのパラメータの事後分布

・P(α):事前分布←適当に決めた。(例:αがベルヌーイ分布のパラメータならば、α=1/2,α=3/4の2値をとるとしてみたときに、p(α=1/2)=1/2、p(α=3/4)=1/2)

・P(D|α):αが決まった時のデータの出現確率→尤度と呼ぶ

P(α|D)=P(D|α)*P(α)÷P(D)

ベイズの定理より!!!

分母P(D)は

P(D)=∫P(D|α)dα

と書き直す。

 

またp(D)はαの関数ではなく定数であるため、無視をして以下のように書く流儀がある。

P(α|D)∝P(D|α)*P(α)

 

補足)事前分布を一様にとると、もちろん上式は

P(α|D)∝P(D|α)

となる、従って、事後確率を最大にするαは最尤推定量となる。

 

一般的にサンプル数が少ないときに有効といわれているが実際のところ良く分かりません。

 

例1:map→mlのNが大きくなる場合の関係)

ある集団の反応率パラメータは何か?ターゲットは1or0でありベルヌーイ分布に従う。

求めたいのは

B(n,α)のα

αの事前分布は経験上こんな感じと仮定する。

α,  事前分布         post[1,] 0.05 0.06250 2.186938e-02[2,] 0.15 0.12500 3.830522e-01[3,] 0.25 0.25000 4.836483e-01[4,] 0.35 0.25000 1.061900e-01[5,] 0.45 0.12500 5.135920e-03[6,] 0.55 0.06250 1.035584e-04[7,] 0.65 0.03125 6.629316e-07[8,] 0.75 0.03125 1.404428e-09[9,] 0.85 0.03125 8.471401e-14[10,] 0.95 0.03125 3.790957e-23

 

今,D=(0,0,1,0,0,0,0,0,1,0,0,0,0,0,1,0,1,0,0,0)が観測された。

最尤法で求める→pハット=4/20=0.2

bayesで確率分布求める。

Nが段々でかくなると、多分確率密度最大点は、最尤推定量に近づく。

薄い点線が事前分布、濃いのが事後分布

補足)Rのplotは結構キレイ!

MCMCとの関連

MCMCはアルゴリズムであり、特定のモデル名ではない。

事後分布から分布の平均を求めたり予測する場合、積分が必要な場合がある。

①:パラメータの事後分布は

P(α|D)=P(α|D)*P(α)/規格化定数

規格化定数=P(D)=パラメータのとりうる値での積分

②:事後分布の平均

∫αP(α|D)d(α)

を結構積分する必要がある。

この積分を簡単にする工夫→共役事前分布:=事後分布と事前分布が同じ分布となる。

例えば、二項分布の場合はパラメータαの分布としてベータ分布を仮定すると、

事後分布もベータ分布になって、計算が楽。回帰分析の係数パラメータの事前分布を正規分布にすると事後分布も正規分布。

しかし、これは計算の都合で事前分布を仮定している!

→もっと柔軟に求めたい。

MC(Monte Calro)→積分不可みたいな、関数の積分を石ころなげて、図形内におちた回数の割合で近似する。円周率もとめるのは有名。

MCMC(Markov Chain Monte Calro)→マルコフ連鎖にしたがうMC.定常分布が求める事後分布となるようにする。(マルコフ連鎖が規約かつ再帰的である必要がある)

ここまで理解すれば少なくともMCMC言われて、知ったかぶらずにすむ。

R2WinBUGS(現在はStanが主流)

RのMCMCパッケージを用いて例題を解いてみよう!

必要なもの

  • winBUGS:MCMC計算ツール
  • R2winbugs:winbugsをRから操作するパッケージ
  • BRugs:同上(今回は使わない)
  • Coda :グラフ作成パッケージ

補足1)アルゴリズムはgibbs sampler と書いてある。

補足2)WinBugsはライセンスをHPからコピーしロックを解除しないと、複雑なモデルガ動かないというエラーが起こる。

IRT項目反応理論への応用

参考URL

http://idecide.mskcc.org/stats/multi/irt_winbugs.php

「http://www.jartest.jp/pdf/3-4kawabata.pdf」MCMCによるテスト項目特性曲線の計算の方法:早大大学院の人のweb上でみれるペーパー。

5つのテスト項目へのresponse(正解、不正解)と被験者の特性を推定する。

2^5で32通りなので行列は32×5とする。

 

 

結果推定量

特性曲線

サンプルされたパラメターの平均値を使用して、特性曲線Pij(jごとのtheta[I]連続関数)書いてみると下記のようになる。

 

正答率75%の項目5はあるていど学力の低い層では他の項目とあまり正答率が変わらなかった。

項目1は賢い層での正答率は高いが、バカな層ではとても低い(二項定理とハサミウチを使う問題みたいなものか?)