pymc3を使ったベイズモデリング入門 Bayesian-Modelling-in-Python　和訳（0-2章）

間違いがある可能性があるので見つけたら指摘していただけると幸いです。

今回紹介するのは機械学習において避けられない考え方であるベイズという考え方です。
下記の記事を基に解説を記述しています。

参考：
https://github.com/markdregan/Bayesian-Modelling-in-Python


#この記事で得られる知識

１：ベイズの考え方
２：ipython notebookを使用してベイズの理論を体感すること
３：ベイズ理論を用いて最適化したモデルの評価方法

#そもそもなぜベイズが必要か

機械学習ではデータからパターンを学習して、未知なるデータに対応することが目的ですが、完璧なデータは

”完全で、一貫していて、正しくて、内容が説明できる”

そのようなケースは非常にレアなので、事前に知識を付与して、変なデータに左右されづらいようにしようという考え方に近いです。

人間も何か行動をするときに予測を立てて行動して、フィードバックを改善していくと思います。

ただし、事前の予測が間違っていれば、よろしくない結果が待っているので、そこは注意が必要です。

例えば、女性にうける肉料理の店を探したい（独断と偏見に基づいています）

![150912426027-thumb-815xauto-19288.jpg](https://qiita-image-store.s3.amazonaws.com/0/10496/63b1a257-7c52-e8d5-90ae-f47f6c0818a5.jpeg)


>*20代前半の女性：未経験なオシャレな店や普段経験できないようなことを提供しているお店に連れて行けば良いだろう

>*20代後半の女性：オシャレな店はある程度、行き尽くしているため、単純にオシャレな店を探しても厳しい。隠れ家的なお店や行きつけのお店があると良いだろう

このように女性にうけると行っても年齢によって事前分布が異なるため、この事前分布を間違えるととんでもないことになります。

#チュートリアル

今回は数式を追っていくのではなく、実際のデータを使用してベイズの考えを理解するチュートリアルが英語版であったので勉強の一環で日本語に翻訳するとともに実践した備忘録を残しておきます。

3章以降はアドバンスドな内容なので２章まで抑えれば、モデルの作成から評価まで行えます。

本編が気になる方は下記をご覧下さい

https://github.com/markdregan/Bayesian-Modelling-in-Python

##データの取得:Section 0. Introduction.ipynb

このチュートリアルでは自身のGoogle Hangoutのデータを取得しています。
データの取得には時間がかかるので他の作業をしながらデータ取得することをお勧めします。

またpython環境が整っていない方は下記の`requirement.txt`を使用すると必要なライブラリはダウンロードできます。

```
numpy==1.9.2
ipython==4.0.0
notebook==4.0.4
jinja2==2.8
pyzmq==14.7.0
tornado==4.1
matplotlib
simplejson
pandas
seaborn
datetime
scipy
patsy
statsmodels
git+https://github.com/pymc-devs/pymc3
```

私はOSXの環境で試したのですが、下記のようなエラーが出ました。

```
RuntimeError**: Python is not installed as a framework. The Mac OS X backend will not be able to function correctly if Python is not installed as a framework. See the Python documentation for more information on installing Python as a framework on Mac OS X. Please either reinstall Python as a framework, or try one of the other backends.
```

対応方法は下記です。

http://stackoverflow.com/questions/21784641/installation-issue-with-matplotlib-python

得られるjsonファイルのネストがすごいので下記にデータフィールドの簡単な情報を載せておきます。

|Field |	Description |	Example
|:-----------|------------:|:------------:|
|conversation_id | 	Conversation id representing the chat thread |	Ugw5Xrm3ZO5mzAfKB7V4AaABAQ |
|participants | 	List of participants in the chat thread |	[Mark, Peter, John]|
|event_id |	Id representing an event such as chat message or video hangout |	7-H0Z7-FkyB7-H0au2avdw|
|timestamp |	Timestamp |	2014-08-15 01:54:12|
|message |	Content of the message sent |	Went to the local wedding photographer today|
|sender 	|Sender of the message |	Mark Regan|

チュートリアルではjsonのデータをパースして、メッセージごとにPandasのデータフレームワークに直しています。

下記のコードでデータを著者とAlison Darcyに絞っているため、ここはコメントアウトしないとデータが出ないので注意です。

```
messages = messages[(messages['sender'] == 'Mark Regan') & (messages['participants_str'] != 'Alison Darcy, Mark Regan')]
```

##解きたい問題

１：レスポンス時間は誰と話しているときに影響している
２：レスポンス時間に影響する要因は？
３：もっとも反応が悪い日は？

私の場合はGoogle Hangoutをあまり使用していないため、下記のような結果になりました。

<img width="798" alt="Screen Shot 2015-11-27 at 10.59.17.png" src="https://qiita-image-store.s3.amazonaws.com/0/10496/42a41312-ebc2-d73d-fcba-9436164acd52.png">

これは今回、解きたい問題には使用できないため、チュートリアルのデータを資料することにしました。
*注意：Export data for usage throughout tutorialで自身のデータを出力するため、チュートリアルのデータを使用したい場合はこの処理は動作させない方が良いです。

##モデルのパラメータの予測：Section 1: Estimating model parameters

ここからベイジアンのチュートリアルの本編です。
図を使用して、結果が分かりやすくなっており、面白いところです。

最初にベイジアンの考え方を述べています。
例として：

>男の子が１日に家の前を通る車の数を数えてノートに記述している。彼のノートには通った車の数が記述されている。

ベイズの考えでは観測されたデータはランダムに発生しているが、それはなんらかの確率分布で発生していると考えます。

例のような離散的なケースではポアソン分布を使用して考えます。

例では平均が5,20,40のケースを載せています。

<img width="713" alt="Screen Shot 2015-11-27 at 11.27.26.png" src="https://qiita-image-store.s3.amazonaws.com/0/10496/e73ed53f-b378-a744-726c-96f7fa3ca326.png">

緑が平均が５の確率分布、オレンジが平均が20の確率分布、ピンクが平均が40の確率分布であります。

前回のレスポンスの時間をポアソン分布の枠組みに当てはめて、ベイズでパラメータを予測することで今回湧いた疑問を解決にチャレンジしています。

<img width="707" alt="Screen Shot 2015-11-27 at 11.30.38.png" src="https://qiita-image-store.s3.amazonaws.com/0/10496/d31a2fba-4c94-b6b6-c54e-91935a631ef5.png">

ポアソン分布の平均値を最尤推定（対数）により推定しています。
今回の推定を行う尤度と推定すべきポアソン分布の平均値の値が下記で確認できます。

<img width="465" alt="Screen Shot 2015-11-27 at 12.52.24.png" src="https://qiita-image-store.s3.amazonaws.com/0/10496/f153f4ad-0460-9a8c-9e35-4aa3087202e1.png">

平均値が20に近い値で尤度がもっとも高い値を取っていることが確認できます。
レスポンス時間のポアソン分布は下記のようになり、18秒で返すことが一番多いです。

<img width="695" alt="Screen Shot 2015-11-27 at 12.55.32.png" src="https://qiita-image-store.s3.amazonaws.com/0/10496/b8b39f5d-5941-40c1-8560-ee2adc38c9ea.png">

今回の例でベイズを当てはめると事前に分かっている情報は10から60の範囲内にデータが収まるということ。それに対するポアソン分布を定義し、それを最尤推定により求めることが主題です。

<img width="341" alt="Screen Shot 2015-11-27 at 12.57.43.png" src="https://qiita-image-store.s3.amazonaws.com/0/10496/fca85b89-eca0-ff35-097f-f04f3b850821.png">

###MCMC

今回の平均値をデータに変動させていき、尤度が最大化するような値になるまで繰り返し行うテクニックです。この手法の良いところはデータがない状態でも事前分布から推定するパラメータを決めてランダムに値を推定していき、尤度が最大化したところで止めることが可能な点が良い点です。

ただし推定するパラメータが多い場合に収束しづらい点や事前分布が適切でない場合は効果を発揮しづらいデメリットがあります。

<img width="658" alt="Screen Shot 2015-11-27 at 13.13.56.png" src="https://qiita-image-store.s3.amazonaws.com/0/10496/5257709d-c942-be20-6c06-8955e73ef954.png">

ipython notebookで動作させるとデータを生成しながらパラメータ推定していく過程が分かるので是非、トライして様子をみてください。

MCMCで実際に推定した結果です。
データは17と19の間で発生し、平均値も18を少し超えただけなので、単純なベイズ推定と同程度の精度が出ています。

<img width="874" alt="Screen Shot 2015-11-27 at 13.20.00.png" src="https://qiita-image-store.s3.amazonaws.com/0/10496/b1359137-faeb-8ab3-407f-2cfbc2395d37.png">

MCMCで最大化されている尤度の軌跡は下記で確認できます。
推定された平均値が必ずしも予想通りに収束しているわけではないのでこのトレースをチェックすることでどのような変遷をおっているか確認をすることができます。　

<img width="688" alt="Screen Shot 2015-11-27 at 13.39.47.png" src="https://qiita-image-store.s3.amazonaws.com/0/10496/5b2186d2-c51d-a1ff-61c6-2f97b57662ad.png">

サンプルで出力する値と今まで出力した値の相関係数も把握する必要があります。

##モデルのチェック：Section 2: Model checking

チェックする観点は２つ

１：モデルはデータを表しているか
２：モデルの比較

データと予測された分布をチェックしてみます。
分布の最頻値と頻度の高いレスポンス時間が一致していません。今回のケースではこのモデルは適していないことが分かります。

<img width="757" alt="Screen Shot 2015-11-27 at 14.23.26.png" src="https://qiita-image-store.s3.amazonaws.com/0/10496/69ae1850-953a-4cb0-1303-d2a22150d6aa.png">

そこでポアソン分布と比較的似ている負の二項分布を使用すると平均だけでなく分散も扱えるようになるため、置き換えてみます。

下記のように分布は似ています。

<img width="736" alt="Screen Shot 2015-11-27 at 14.18.08.png" src="https://qiita-image-store.s3.amazonaws.com/0/10496/6968513a-2168-ea14-e1c6-292ca0bb9a9d.png">

負の二項分布は下記のようにαとμのパラメータを推定します。

<img width="424" alt="Screen Shot 2015-11-27 at 14.25.39.png" src="https://qiita-image-store.s3.amazonaws.com/0/10496/d1269c16-8f68-4ee7-0971-7b1bfa2bc1af.png">

予測された値が下記となっており、αが1.4から2.4の値をとっており、分散の役割によって分布の表現力が上がっています。

<img width="884" alt="Screen Shot 2015-11-27 at 14.29.36.png" src="https://qiita-image-store.s3.amazonaws.com/0/10496/a0355533-b11b-d350-5fdd-4f3b1aee81d2.png">

下記が先ほど推定したαとμのパラメータを使用して作成した分布とレスポンス時間の図である。
分布はレスポンス時間の分布と似ており、よりその特徴を表しています。

<img width="744" alt="Screen Shot 2015-11-27 at 14.34.35.png" src="https://qiita-image-store.s3.amazonaws.com/0/10496/942e7d73-ab78-4784-e3a0-de56715499f9.png">

ポアソンと負の二項分布を組み合わせて使用する方法も提案されています。

<img width="552" alt="Screen Shot 2015-11-27 at 14.46.22.png" src="https://qiita-image-store.s3.amazonaws.com/0/10496/d972070c-f096-aded-1da2-eaa4c1ad2ea4.png">

ベイズファクターを計算して、下記の基準にそってどっちのモデルを使用すれば良いか決めると記述されています。

<img width="366" alt="Screen Shot 2015-11-27 at 14.48.19.png" src="https://qiita-image-store.s3.amazonaws.com/0/10496/451afb9d-5719-9472-e013-6f4d21b6b84a.png">

今回は2章までの基本的な所を行いましたが3章からは発展版になるので是非チャレンジしてください。

https://github.com/markdregan/Bayesian-Modelling-in-Python/blob/master/Section%203.%20Hierarchical%20modelling.ipynb

#参考

ぱくたそ
https://www.pakutaso.com/

Stanで統計モデリングを学ぶ(2): そもそもMCMCって何だったっけ？
http://tjo.hatenablog.com/entry/2014/02/08/173324

ベイズ推定を知っているフリをするための知識
http://www.anlyznews.com/2012/01/blog-post_31.html

Bayesian-Modelling-in-Python
https://github.com/markdregan/Bayesian-Modelling-in-Python

#次回の記事

[tksarah](http://qiita.com "tksarah")さんによるOpenStack LIBERTY で DRBD Cinder Volume Driver 動作確認です!!