MCMCについて学ぶ

今日からMCMCについて学んでいこうと思う。今授業でベイズを学んでいてベイズphilosophyに感動している。 経済学とMCMCの関連性も高いという事でしっかりと学んでいきたい。プログラミング・線形代数・統計に関しての知識は初心者レベルなので一つ一つ丁寧に理解していきたいと思います。

ベイズに惚れた理由

現在ビッグデータが支配する社会となり、大量にあるデータから、適切に必要なデータを取得しinsightを得るかが重要になっております。この時代背景が、私が統計学を始めたきっかけです。その一方で今世の中にあふれているデータは果たして分析に向いているでしょうか？big data = dark data = irrevant data となるケースが良く見られます。また、データ化はされていないが長年の経験から蓄積され人間の頭の中に存在するノウハウ（知恵）の活用も見逃せません。
こういった点を考慮し、事前分布の設定（ノウハウ）で不十分なデータを補うことができる可能性があるのがベイズ統計（ベイズモデリング）の素晴らしさです。またベイズはgiven dataを最大限にrespect しています。従来の統計学とは異なりデータが加わるごとに事後分布を更新します。変化の激しい社会だからこそ、頻度論が提唱する真の値から、データを予測するモデルではなく、未知の値を確率変数と捉えて観測されたデータ（結果）から母集団（原因）を推定するベイズ統計は注目されていくべきです。

NOTE

なぜ最尤推定法が好きではないのか。

→過学習の可能性がある。事前分布（ノウハウ）を考慮していない。

点推定はそもそもおかしい（適切な推論は、予測結果だけでなく、そのばらつき（分散）も考慮する必要がある）

MCMCがなぜ登場した？

授業ではposterior distribution をprior とlikelihood から求める計算をしてきました。また適切な事前分布の決め方として、無情報事前分布（Jeffery prior）や自然共役事前分布を学んできました。モデルの設定次第では事後分布はnamed distribution(betaやnormal)になることもありますが、ほとんどはunnamed distribution であり、実務面で利用することは難しいです。

そこでコンピュータが発達した現代では乱数シミュレーションや最適化を用いた近似値推論が主要な解法となってます。その中でも代表てきなMCMC法を用いれば、事後分布に従う乱数を生成することができ、複雑な積分計算を回避できます。

数式や導出過程についてまとめると沼にはまるのでgoogle検索で確認する方が好ましいです。実際にpythonを用いてどうなっているのか見ていきましょう。

# サンプルデータの作成
import numpy as np
import matplotlib.pyplot as plt
np.random.seed(112)

# サンプルデータ生成のパラメーター
N = 100
alpha_real = 2
beta_real = 1.2
eps_real = np.random.normal(0, 0.5, size=N) #ノイズに関して平均0 分散0.5^2 の正規分布

# モデル: y = α + βx + Noise
x = np.random.normal(10, 1, N)　#平均10分散1の正規分布からランダムにsampling 
y_real = alpha_real + beta_real *x
y = y_real + eps_real

# データの可視化
plt.plot(x, y, '.')
plt.xlabel('x')
plt.ylabel('y')
plt.plot(x, y_real, 'orange')
plt.show()

線形モデルは シンプルなy = α + β*xで定義しており、α=2、β=1.2が今回のパラメータになってます。

MCMC手順
①モデルの定義　②MCMCの実行
⓷パラメータの事後分布の確認・評価
⓸事後分布から推論（credible interval の可視化）

pymcがpython で使う事が出来る

import pymc3 as pm    

# モデルの定義
with pm.Model() as model:
   # 1. 事前分布
   alpha = pm.Normal('alpha', mu=0, sd=10)# αには平均0, 標準偏差10の正規分布を仮定
   beta = pm.Normal('beta', mu=0, sd=1)# βには平均0, 標準偏差1の正規分布を仮定
   epsilon = pm.HalfCauchy('epcilon', 5)#  noiseには最頻値5の半コーシー分布を仮定
   # 2. 尤度の計算
   y_pred = pm.Normal('y_pred', mu=alpha+beta*x, sd=epsilon, observed=y)#yは観測データ→これを基に更新していく
   # 3. MCMCの実行
   trace = pm.sample(11000, chains=1)　＃iterationの数, チェーン数の指定

# パラメータの事後分布の確認・評価（グラフ）
trace_n = trace[1000:]# 初めの1000件を捨てる <- burn in の指定
pm.traceplot(trace_n)

# 各パラメータの事後分布の確認・評価（統計量）
pm.summary(trace_n)

iteration: どれだけの数サンプルを取るのか

chain : MCMCを何回実装するのか

burn-in: 最初の数千のサンプルはinitial設定に影響されており好ましくないのでcutする

確認すべき3つのグラフ

①ts.plot(トレースプロット)：MCMCの乱数を順番にプロットしたもの
• burn-inの期間やMCMCの質の良さを確認する

②histgram or posteior pdf:  基本的にはbell shapeになっているはず

⓷acf(自己相関図):横軸にラグ、縦軸に自己相関をとる
• MCMCのサンプルはiidなので、横軸の増加とともに自己相関が低くなるはずである。

今回の結果

alphaの推定があまり上手く行っていません。　

トレース図から定常状態への収束は確認できました。また自己相関も問題がありませんね。

itrationを増やしてみます。

同じ結果で上手く推定できていません。観測データ100が十分ではないのでしょうか？100->1000に変えてみます。

近づきましたが未だに少し遠いですね。今回は原因が分かりませんでした（もう少し勉強します）

まとめ：

alphaに関してはサンプル100では上手く推定できなかったが、betaに関しては近似推定ができた。MCMCを用いればサンプルを無限に取り出すことで、必要なデータがなくて、DXを推進できないと悩む企業もデータ活用を推進できるのではないでしょうか。また今回は正規分布を仮定しが、それ以外にも多くの事前分布の設定がある。多くのlegacy企業が大切にしている長年の培った知識・知恵をデータと反映させることが出来るのではないでしょうか？

MCMCの基礎を勉強してみました。シンプルなモデルから入ったのでまだまだ学ぶことがたくさんです。じっくり学んでいきたいです。

参考：https://qiita.com/takubb/items/a2e02fe61aad80f8f2f3