ベイズファクターとブリッジ・サンプリング #R

本記事はClassi AdventCalendar2020 11日目の記事です。
昨日は@kasaharuさんが書いてくれました。
Angular アプリケーションの E2E テストのフレームワークを Cypress に変更するまで

こんにちは。データ/AI部で部長をしている@tetsuroito です。今年もこの季節がやってきましたね！
今年は暦と季節が非常にマッチしていた1年だったように思います。なので、冬は寒いですが、みなさまお元気でしょうか。

実はClassiでは、R&D活動として大学との共同研究を精力的に行っていて、過去にも様々なプレスリリースなどを出しています。

教育プラットフォーム「Classi」と大阪経済大学が社会人基礎力に関する調査を実施
 機械学習の専⾨家東京⼤学の佐藤⼀誠⽒がClassiのアドバイザーに就任
 AI×教育の第一人者電気通信大学の植野真臣氏がアドバイザーに就任

今年、私が特に関わっていたのは、東京大学の岡田先生の研究室との取り組みです。
岡田先生は東京大学の教育学部に所属されていますが、専攻が心理統計ということもあり、Classiのデータと心理統計を掛け合わせたテーマで共同研究を行っています。

今年の共同研究発表リスト

今年、共同研究を行った中で、対外的に発表したものは下記の通りです。

問題の解き直しの学業成績への効果-階層線形モデルによるビッグデータの分析- 菱山、伊藤、岡田行動計量学第47巻2号2020、
Bayesian penalization for variable selection in explanatory cognitive diagnostic models Yoshito Tan,Tetsuro Ito.Kensuke Okada International Meeting of the Psychometric Society 2020
回答者の背景情報を考慮した説明的認知診断モデル丹、伊藤、岡田　日本行動計量学会第48回全国大会
補償型の説明的認知診断モデルの提案　丹、伊藤、岡田日本分類学会第39回大会

上記の4本です。コロナによって、学会やイベントなどが中止やオンライン化されたため、現地での発表というものはなく、誌面発表やオンラインでの発表という体裁を取る珍しい1年だったと思います。日本分類学会では、共著者の丹くんが学生優秀発表賞を受賞するなど、評価もいただけました。
また、行動計量学のジャーナル投稿では査読のプロセスもあり、私個人としては新しいチャレンジもできたという意味では有意義な一年でした。

そんなわけで、今年一年はけっこうベイズ統計モデリングに触れる機会も多かったので、本記事ではベイズ統計モデリングのモデル評価に関する話題を紹介します。

モデリングとベイズ統計における分析のステップ

まずはじめに、用語の説明を少ししておきましょう。

モデル
- 観測したデータを生み出す確率的な過程を簡潔に記述したもの
モデリング
- 観測したデータをもとにモデルを作る行為。
ベイズ統計モデリング
- 統計的データ分析に基づく推論における不確実性を定量化するために、確率を明示的に扱うモデリング

モデリングを行う目的は観測されたデータからモデリングの過程を通して、現象の解釈や将来の予測を行いたいというモチベーションからきます。つまり、データ発生のメカニズムについて思いを馳せることができるかもしれないのです。
私たちはベイズ統計モデリングの性質を使って、人が学習する効果を適切にモデリングできないか、その現象はどのような解釈なのか、どうしたらよくなるのかなどを知りたいと思って使っています。

次にベイズ統計学に基づくデータ分析のステップを示します。

データやその収集過程に関する知識に基づき、確率モデルを作る
観測されたデータを条件としてモデルのパラメータや予測値などの事後確率分布を計算する
データに対するモデルの当てはまりや、計算された事後確率分布の評価を行う

大別すると、この3ステップに分けられます。
近年ではJAGSやStan、tensorflow probabilityなどの確率的プログラミング言語が用いられ、これらの確率モデルやシミュレーションが実行しやすい環境が整ってきています。
ここで重要になるのは、そのシミュレーションを実データと照らし合わせて確率モデルやパラメータを評価することです。

モデル評価のための指標

さて、確率モデルを評価するのはどうしたら良いのでしょうか。
良いモデルとは、「真の分布から生成されるデータをうまく説明、予測できるモデル」であるといえるでしょう。
これを比較することができるのが、情報量です。確率モデルと生成されたデータの分布の情報量を計算して、その大小を計算することで、その確率モデルの良し悪しを測ることができます。

2010年に東工大の渡辺澄夫先生がWAIC、WBICという特異モデルの指標を発表して10年が経過しましたが、今ではこれがデファクトになっています。

Asymptotic Equivalence of Bayes Cross Validation and Widely Applicable Information Criterion in Singular Learning Theory

また、最近の学会の発表ではよくベイズファクターという指標についての議論もよく出てきます。
このベイズファクターとは一体なんでしょうか。

ベイズファクターは2つのモデルの周辺尤度比として、

BF_{10}= \frac{p(x^n|M_1)}{p(x^n|M_0)}

と定義できます。また、自由エネルギーを使って

BF_{10}= exp(F_{M_0}-F_{M_1})

とも記述できます。
評価基準は特に定まっていませんが、3以上であればモデル1がやや優勢、20以上でかなり優勢と判断されたりするようです。(社会科学のためのベイズ統計モデリング朝倉書店より)
ただし、このベイズファクターは解析的に解くことが非常に難しく、計算が容易ではありません。そこで、近年ではブリッジ・サンプリングという手法での近似計算が主流になっています。

ブリッジ・サンプリングとは

ブリッジ・サンプリングはモンテカルロ積分を利用し、その乱数から積分を計算し、近似解を求める手法のようです。
理論的な紹介は、共同研究をしている岡田先生が詳細に書かれた論文が公開されているので、そちらに譲ります。

ベイズファクターによる心理学的仮説・モデルの評価

Rにはこのブリッジ・サンプリングを手軽に計算することができる{bridgesampling}パッケージが公開されています。CRANに登録されているため、install.packagesで利用できます。

簡単なstanコードを用意して、計算してみましょう。

data{
    int N;
    int X[N];
    real a;
    real b;
}

parameters{
    real<lower=0> lamda;
}

model{
    target += poisson_lpmf(X | lambda);
    target += gamma_lpdf(lamda | a,b);
}

ブリッジ・サンプリングの計算

library(bridgesampling)
model.fe <- stan_model("FE_poisson.stan")
fit.bs <- sampling(model.fe,data=list(N=n,x=x,a=a,b=b),iter=10000,warmup=1000,chaine=4)
bs <- bridge_sampler(fit.bs,method="wrap3")
-logml(bs)

注意として、予めサンプリングの数を決めておかないといけません。引数に設定するためです。
しかし、関数で一発で書くことができるため、非常に簡単ですね。
これによって、様々な確率モデルを評価していくことができます。

色々な事象のメカニズムを推測し、手元に取れるデータと真のモデルからのサンプリングしたモデルの情報量を比較して、あなたも明日からベイジアンモデラーになれますね！
Enjoy！！

あと！書籍が出ます！kindle版はすでに買えますので、ぜひ！

前著データサイエンティスト養成読本ビジネス活用編から2年が経ちましたが、12/21に次回作「AI・データ分析プロジェクトのすべて」を発売します！前回から2年の経験や学びを少しでもお伝えしたくて書きました！皆さん、ぜひよろしくお願いします！https://t.co/l7r0kB7jBX
— 機械の体を手に入れるのよ鉄郎 (@tetsuroito) November 29, 2020

明日は @kitaharamikiya さんの番です。楽しみですね！