こんにちは、情報科学技術研究所のy_kawasakiです。春になり天気予報が微妙にずれる（時間的に）のが気になる今日このごろです。
ベイズ統計モデリングとは
ベイズ統計学による、統計モデリングのことです！それぞれの意味は個々でググってください！
準備
ベイズ統計モデリングソフト
ベイズ統計モデリングをやるには、MCMC（マルコフ連鎖モンテカルロ法）で解く必要があるので様々なパッケージが用意されています。たとえば、BUGS言語を使ったソフトや、Stanというツールなど。ここでは、最近流行りのStanを使うことにします。Pythonから使う必要が（自分には）あるため、PyStanを利用します。
PyStanのインストール
> pip install PyStan
難しいことはありません。
実例
弊社では本社があるビルが手狭になったことから一部の部署が切り出されて、近くのビルに放り出されました（以下、本社とサテライト）。そこで、時々、本社と行き来する必要が発生して、片道どれくらいかかるのかを情報科学技術研究所的には調査する必要性が発生しました。（使命感）
そこで、被験者3名に計測を依頼したところ、A7件、B2件、C4件、計13件の下記の結果を得ることができました。なお、本社、サテライト間には、信号とエレベータが障害として存在しています。
| 試行回 | 被験者 | タイム[s] |
|---|---|---|
| 1 | A | 204.2 |
| 2 | A | 228.5 |
| 3 | A | 256.2 |
| 4 | B | 199.5 |
| 5 | B | 222.2 |
| 6 | A | 214.1 |
| 7 | A | 294.5 |
| 8 | C | 304.6 |
| 9 | C | 197.1 |
| 10 | A | 318.9 |
| 11 | A | 255.8 |
| 12 | C | 351.2 |
| 13 | C | 314.5 |
基礎統計的データ
今回13回の試行から、平均258.6、標準偏差50.5、最大351.2、最小197.1を得ました。今回は、この結果を真の統計分布とします。
結果
事前分布として、無情報分布を仮定し、事後分布として正規分布に従うと仮定して、ベイズ統計モデリングを行ってみました。
| 試行回 | 平均 | 標準偏差 | ベイズ推定平均 | ベイズ推定標準偏差 | ベイズ推定95%区間 |
|---|---|---|---|---|---|
| 3 | 229.6 | 21.2 | 245.3 | 152.3 | 24.0 - 698.8 |
| 4 | 222.1 | 22.5 | 222.1 | 53.0 | 161.0 - 277.9 |
| 5 | 222.1 | 20.2 | 222.2 | 37.1 | 181.4 - 268.8 |
| 6 | 220.8 | 18.7 | 221.3 | 29.5 | 194.5 - 247.0 |
| 7 | 231.3 | 31.0 | 232.7 | 45.2 | 192.6 - 275.8 |
| 8 | 240.5 | 37.8 | 240.6 | 49.7 | 199.9 - 275.3 |
| 9 | 235.7 | 38.2 | 235.6 | 48.2 | 201.0 - 266.3 |
| 10 | 244.0 | 43.9 | 244.2 | 54.4 | 209.1 - 277.1 |
| 11 | 245.1 | 42.1 | 245.1 | 50.6 | 215.1 - 275.3 |
| 12 | 253.9 | 49.8 | 254.8 | 59.0 | 221.1 - 288.7 |
| 13 | 258.6 | 50.5 | 258.9 | 59.0 | 224.7 - 291.2 |
この表で示されるように、今回のような単純な場合は、ほぼ古典的な統計手法と同様の結果を得ることができました。それだけでは、何の嬉しさもないんですが、注目すべきは、95%区間がベイズ推定によって与えらることです。
それによれば、4分弱から5分弱で到着できることがわかります。そういう意味で言うと標本数が4の時点で161秒から278秒であることを示しているのは驚異的であるように感じます。
このように、ベイズ推定モデリングは標本数が少ないときにもそれなりの結果を示してくれるという、利点を持っています。
結果2
人間には感覚的にわかることがあります。例えば、本社サテライト間は５分前後だな、と。というわけでそれをモデルに組み込んであげましょう。そういうことができるんです、そうベイズ推定モデリングならね！
（ただし、無情報分布を事前分布として設定することが推奨されているようないないような。。。）
| 試行回 | ベイズ推定平均 | ベイズ推定標準偏差 | ベイズ推定95%区間 |
|---|---|---|---|
| 2 | 261.1 | 375.0 | 185.9 - 378.5 |
| 3 | 251.9 | 86.6 | 187.7 - 348.4 |
| 4 | 234.0 | 54.5 | 189.1 - 310.8 |
| 5 | 228.1 | 37.3 | 196.1 - 274.6 |
| ...以下略 |
のように、それなりな事前分布を与えてあげることで、早いうちにそれなりの結果を与えてくれます。ただし、恣意的に操作していると考えられなくもないことから、無情報分布を与えたほうがいいという考え方があるのだと思われます。
今後の課題
信号待ちやエレベータ待ちを考慮したモデリングも面白そうです。
参考図書
StanとRでベイズ統計モデリング (Wonderful R)
- 作者: 松浦健太郎,石田基広
- 出版社/メーカー: 共立出版
- 発売日: 2016/10/25
- メディア: 単行本
- この商品を含むブログ (8件) を見る
- 作者: 岩波データサイエンス刊行委員会
- 出版社/メーカー: 岩波書店
- 発売日: 2015/10/08
- メディア: 単行本（ソフトカバー）
- この商品を含むブログ (9件) を見る
付録
import pystan stan_code = """ data { int <lower=1> N; real T[N]; } parameters { real avg; real <lower=0> sigma; } model { //avg ~ normal(300, 60); T ~ normal(avg, sigma); } """ model = pystan.StanModel(model_code=stan_code) d = [204.2, 228.5, 256.2, 199.5, 222.2, 214.1, 294.5, 304.6, 197.1, 318.9, 255.8, 351.2,314.5] stan_data = { "N": len(d), "T": d, } fit = model.sampling(data=stan_data, iter=1000, chains=4,seed=1234) print(fit)