##自己紹介
データサイエンスを専攻している大学生です。コンピューターサイエンスをメインに学んでるわけではないので多少知識はあれども、明らかに知識不足です。
画像解析に興味を持ち勉強していますが、数理統計学の知識を深くしたいという思いから統計検定準一級合格を目指して勉強中です。
##本題
今回はベイズ統計学とシミュレーション技法を勉強中に触れたMCMC(マルコフ連鎖モンテカルロ法)、特にその中でもメトロポリスアルゴリズムについての理解を深めようと思い、試しにPythonでプログラミングしてみたので記事にしました。
具体的な内容はタイトル通り去年の統計検定準一級の問7で出題されたベイズ統計学の問題をメトロポリスアルゴリズムを用いて考えてみようというものになっています
「MCMCやベイズ統計学について勉強してみたけどよくわかんない」みたいな方が参考にできるような記事を書きたいなと思いますので、逆にベイズ統計学やMCMCを熟知した専門の方からすると面白くもなく、間違いなどが目立つかもしれません。そのようなときは是非ご指摘よろしくお願いします。
#環境
・windows
・python 3.5.6
##ベイズ統計とは
ベイズ統計学というのは統計学の分野の一つです。ベイズの定理の考え方を確率分布に応用することで成り立っている分野と私は解釈しています。（だいぶ簡潔にとらえすぎな気がします）
ベイズ統計学では「事後分布」と「事前分布」と「尤度」という言葉が出てきます。
事後分布\propto 事前分布 \times尤度
これがベイズ統計学のすべてとなる数式です。この事後分布というのはあるパラメータの確率分布になっています。このパラメータの確率分布を導出するのがベイズ統計学です。
では事前分布とは何か？ 事前分布とは事前に考えることができるパラメータの確率分布です。例えば過去のデータから「恐らくこのパラメータは正規分布に従うな」みたいな感じで予測したものでもいいですし
「常識的に考えてこの分布に従うだろ」という決めつけでもいわけです。
そして、尤度というのはパラメータに応じて与えられたデータが発生する確率を確率密度関数から導出したものです。
一番わかりやすい例え話は、コインを題材にしたものです。
ここにコインがあります。コインの表が出る確率としてふさわしいものを選びたいのですがコインを投げる回数は10回までとしか決められていません。ふさわしいものをどう選びましょう。
この問題で頻度主義的な最尤法を用いて考えるとどうなるでしょう？ 今回10回のベルヌーイ試行で3回表が出たとき、コインの表が出る確率をPとして考えてみましょう
L(P)={}_1{}_0C_3P^3(1-P)^7
こいつをPで微分し、極値を求めるとP=0.3
となるためコインが表になる確率は30％となるわけです。
しかし、この結果を改めてみると「ホンマか?」となるわけです。理由は簡単です。データ数が足りません。
コインなんて大体表になる確率が50％でもたまたま10回中3回しか表にならないなんてありえそうなものです。なのに最尤法を使うとパラメータの推定値は0.3と直感よりも少なめに出てきてしまいます。
そこでベイズ統計学の登場です。
先ほど記述した事後分布の比例式を見ると事前分布と尤度の積に比例していると書いてある。この時、事後分布というのはパラメータPの分布です。「パラメータP(つまり表の出やすさ）が0.5から0.6である確率は～ですよ。」というようなことわかる分布である必要があります。同様に事前分布もPの分布です。事前分布は一度書いたがなんでもいいのです。そこで、この場合β分布を定義します。具体的にはBe(50,50)を採用する。この分布は実際に記述してみます。
import matplotlib.pyplot as plt
from scipy import stats
import numpy as np
x=np.linspace(0,1,100)
beta_pdf1=stats.beta.pdf(x,50,50)
plt.plot(x,beta_pdf1,label="a=50,b=50")
plt.show()
これを見てもらうとわかると思いますがBe(50,50)というのはP=0.5付近での確率が高くなりがちな分布です。
つまり我々がコインの表のでやすさPを考えるとき事前に「大体50％くらいだろ」と考える感覚を反映した分布ともいえます。さて、ここからは上に書いた式を利用して事後分布を導出してみます。
\begin{align}
&事前分布=Be(50,50)=\frac{1}{B(50,50)}P^{50-1}(1-P)^{50-1}\\
&尤度={}_1{}_0C_3P^3(1-P)^7\\
\end{align}
よって、事後分布は以下のようになる
\begin{align}
事後分布&\propto \frac{1}{B(50,50)}P^{50-1}(1-P)^{50-1}{}_1{}_0C_3P^3(1-P)^7\\
&\propto P^{52}(1-P)^{56}
\end{align}
この時、最後の式部分以外は定数部分になるのであとで考えればよい。
重要なのはパラメータPを確率変数とした事後分布のPによって変化する部分は上の式で表すことができるということである。あとは事後分布を具体的に導出する
\begin{align}
\int_{0}^{1}t^{a-1}(1-t)^{b-1}dt &=B(a,b)\\
&=\frac{Γ(a)Γ(b)}{Γ(a+b)}
\end{align}
であるからして、事後分布は確率分布である(つまり定積分布した時面積が１でなければならない）ということを利用して
\begin{align}
事後分布&=\frac{1}{B(53,57)}P^{52}(1-P)^{56}\\
&=Be(53,57)
\end{align}
となることが分かり、事後分布もベータ分布となることがわかる。
ただ、「事後分布もベータ分布なんだ！」といわれても「は？だから何？何がうれしいの？」って感じだと思います。(ていうか私も思いました。） なので、この情報を視覚化します。
import matplotlib.pyplot as plt
from scipy import stats
import numpy as np
x=np.linspace(0,1,100)
beta_pdf1=stats.beta.pdf(x,50,50)
beta_pdf2=stats.beta.pdf(x,52,56)
plt.plot(x,beta_pdf1,label="a=50,b=50")
plt.plot(x,beta_pdf2,label="a=53,b=57")
plt.legend()
plt.show()
ほとんどさっきの使いまわしですが、事前分布と事後分布を同時に描画しています。
オレンジが事後分布です。どうでしょうか？少しだけ左に寄っているのがわかります。
つまり事後分布は事前分布より少し小さいPの値の付近で大きな値を出力するわけです。簡単な話、たった10回の試行で3回表が出たから安易に30％と考えるのでなく、コインというのは50%ぐらいで表が出るという事前情報を考慮しながら、10回の実験の結果も無視せずに反映しているような分布を作り出すことができているわけです。これがベイズ統計学です。
#メトロポリスヘイスティング法とは
次に、メトロポリスヘイスティング法(M-H法)というシミュレーション技法について私の理解できている範囲で書いていきます。 M-H法というのはマルコフ連鎖モンテカルロシミュレーション(MCMC)の一種です。MCMCっていうのは簡単に言えばある一定の確率Pでパラメータθが新しいθに更新されるような(マルコフ連鎖的な更新）を疑似乱数を使用したモンテカルロシミュレーションと組み合わせたものと自分はとらえています。
その中の一種であるM-H法のアルゴリズムは以下の通りです
\begin{align}
&(1)初期値q_0を設定\\
&(2)ランダムにq_1の候補を選ぶ\\
&(3)q_0の時の尤度L(q_0)とq_1の時の尤度L(q_1)を計算\\
&(4)L(q_1)>L(q_0)ならばq_1を採用して(1)に戻る\\
&(5)P=\frac{L(q_0)}{L(q_1)}として、Pの確率でq_1を採用して(1)に戻る
\end{align}
このアルゴリズムを繰り返すことにより、事後分布の標本を疑似的に手に入れることができるわけです。
つまり、M-H法をはじめとするMCMCは標本を疑似的に入手し、関数として導出することが困難な事後分布の形を疑似サンプルから描くことを目的としたシミュレーショなのです。
一見アルゴリズムが最適化手法に似ているからといって最適化手法と勘違いしてはいけません。(私は初め違いが分かりませんでした)
#コード
それでは実際に検定試験の問題からコードを作成していきます。
試験問題は統計検定のホームページで公開中なので見てください。
import numpy as np
from scipy import stats
import random
class MCMCmetropolic:
def __init__(self,x,n_iter=3000,premean=0,presigma=1,postsigma=2):
self.n_iter=n_iter
self.x=x
self.premean=premean
self.presigma=presigma
self.postsigma=postsigma
def pre(self,mean):
premean=self.premean
presigma=self.presigma
return stats.norm.pdf(x=mean,loc=premean,scale=presigma)
def log_likelihood(self,mean):
likelihood=0
x=self.x
if type(x)==list:
for i in range(len(x)):
likelihood=likelihood+np.log(stats.norm.pdf(x[i],mean))
else:
likelihood=np.log(stats.norm.pdf(x,mean))
return likelihood
def post(self,mean):
return self.log_likelihood(mean)+np.log(self.pre(mean))
def main(self):
mean=np.random.randn()
post=self.post
n_iter=self.n_iter
mean_list=[]
post_list=[]
for i in range(n_iter):
newmean=np.random.normal(loc=mean,scale=0.1)
selection_p_1=post(newmean)-post(mean)
alpha1=min(1,np.exp(selection_p_1))
r1=random.uniform(0,1)
if r1<alpha1:
mean=newmean
post_=np.exp(post(mean))
mean_list.append(mean)
post_list.append(post_)
return mean_list,post_list
以上がメトロポリスアルゴリズム
これに問7の説明文にある「例えば～」という部分の条件に当てはめてサンプリングすると
import matplotlib.pyplot as plt
from metropolic_hasting import MCMCmetropolic
def main():
x=[2]
metro=MCMCmetropolic(x=x,premean=0,presigma=1,postsigma=2)
mean,post=metro.main()
plt.plot(mean,post,"o")
plt.show()
if __name__=="__main__":
main()
[2]の条件に当てはめてみると
import matplotlib.pyplot as plt
from metropolic_hasting import MCMCmetropolic
def main():
x=[14.05,19.25,23.0,16.0,13.9,14.7,20.35,15.05,15.3,15.5]
metro=MCMCmetropolic(x=x,premean=13,presigma=2.7**2,postsigma=3.07**2)
mean,post=metro.main()
plt.plot(mean,post,"o")
plt.show()
if __name__=="__main__":
main()
出力結果です。
ただ、これだけ見てもM-H法のアルゴリズムの存在意義がいまいちピンとこないかもしれません。そこで、シミュレーションの動きを実際に見てもらおうと思います。
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.animation as animation
from metropolic_hasting import MCMCmetropolic
def main():
fig = plt.figure()
ims = []
x=[2]
metro=MCMCmetropolic(x=x,premean=0,presigma=1,postsigma=2)
mean,post=metro.main()
for i in range(len(mean)):
im=plt.plot(mean[i],post[i],"o")
ims.append(im) # グラフを配列 ims に追加
ani = animation.ArtistAnimation(fig, ims, interval=100)
ani.save("output.gif", writer="imagemagick")
if __name__=="__main__":
main()
どうでしょう？わかりますか？
このアルゴリズムでは山の頂上付近の値が頻繁に登場します。逆に山の裾部分の値はほとんど登場しません。
この動点の登場回数はいわば標本数なわけです。正規分布を仮定した標本調査を想像してほしいのですが、正規分布の平均付近のサンプル数は当たり前ですが多くなる傾向になることは何となく想像がつくと思います。
M-H法はそのサンプリングを疑似的に再現しているシミュレーションなわけです。
今回題材にした問題はベイズ統計学の「共役事前分布」の性質が絡んでくる問題で、簡単な話、事前分布と尤度が正規分布なら事後分布も正規分布という性質です。なので、サンプリングした結果は正規分布になっています。 この場合、解析的に解くことが可能なので実はM-H法を適用する必要がありません。しかし、事前分布と尤度の形によっては解析的に事後分布を求めることができない場合があります。その場合、疑似標本から分布の形を推測する必要があります。その際に活躍するのがM-H法ということなんですね。
では実際にM-H法を活用して、解析が困難な事後分布をサンプリングしてみましょう。
今回参考にするのは私がM-H法を勉強したデータ解析のための統計モデリング入門を使いたいと思います。
この本でM-H法について触れている章ではポアソン回帰のパラメータをベイズモデル化することで考えようとする方法をとっています。その際、当たり前ですが尤度はポアソン分布を仮定しています。事前分布に採用しているのはとても平べったい正規分布です。
描画すると
import matplotlib.pyplot as plt
from scipy import stats
import numpy as np
x=np.linspace(-1,1,100)
beta_pdf1=stats.norm.pdf(x,loc=0,scale=100)
beta_pdf2=stats.norm.pdf(x,loc=0,scale=0.1)
plt.plot(x,beta_pdf1,label="sigma2=100")
plt.plot(x,beta_pdf2,label="sigma2=0.1")
plt.legend()
plt.show()
青いグラフが事前分布です。正規分布には一瞬見えませんが立派な正規分布です。分散を大きくすることで平べったくしているわけです。なぜこのような分布を仮定するかというとパラメータの分布に見当がつかないからです。 上で書いたようなコインの話であれば表の出やすさの分布は何となく想像がつきますが、回帰パラメータの分布はあまり想像がつきません。なので、どんな値をとる可能性もほとんど等しいような(違いがほとんど出ない)ような分布を仮定しておきたいわけです。そこで採用されたのがこの平べったい正規分布です。これであれば-∞から∞までどんな値をとる可能性もほとんど等しいとみることができます。
こういう事前分布を無情報事前分布といいます。この平べったい正規分布以外にも多くの種類がるので調べてみてください。
では、早速この無情報事前分布を活用して、先ほどのコードを改変しながら見ていきます。
import numpy as np
from scipy import stats
import pandas as pd
import random
import matplotlib.pyplot as plt
class MCMCmetropolics:
def __init__(self,x,y,n_iter=3000,burn=1000):
self.n_iter=n_iter
self.x=x
self.y=y
self.burn=burn
self.mean=np.mean(x)
def pre_B1(self,B1):
return stats.norm.pdf(x=B1,loc=0,scale=100)
def pre_B2(self,B2):
return stats.norm.pdf(x=B2,loc=0,scale=100)
def lamda(self,B1,B2):
lamda=(np.exp(B1+B2*(self.x-self.mean)))
return lamda
def log_likelihood(self,B1,B2):
lamda=self.lamda(B1,B2)
likelihood=np.sum(np.log(stats.poisson.pmf(self.y,lamda)))
return likelihood
def post(self,B1,B2):
return self.log_likelihood(B1,B2)+np.log(self.pre_B1(B1))+np.log(self.pre_B2(B2))
def main(self):
post=self.post
n_iter=self.n_iter
burn=self.burn
B1_list=[]
B2_list=[]
post_list=[]
B1=np.random.randn()
B2=np.random.randn()
for j in range(n_iter):
newB1=np.random.normal(loc=B1,scale=0.01)
selection_p_1=post(newB1,B2)-post(B1,B2)
alpha1=min(1,np.exp(selection_p_1))
r1=random.uniform(0,1)
if r1<alpha1:
B1=newB1
newB2=np.random.normal(loc=B2,scale=0.01)
selection_p_2=post(B1,newB2)-post(B1,B2)
alpha2=min(1,np.exp(selection_p_2))
r2=random.uniform(0,1)
if r2<alpha2:
B2=newB2
post_=np.exp(post(B1,B2))
B1_list.append(B1)
B2_list.append(B2)
post_list.append(post_)
del B1_list[0:burn]
del B2_list[0:burn]
del post_list[0:burn]
return B1_list,B2_list,post_list
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from metropolics import MCMCmetropolics
data=pd.read_csv("data3.csv")
x=np.asarray(data["x"])
y=np.asarray(data["y"])
metro=MCMCmetropolics(x=x,y=y)
b1,b2,post=metro.main()
plt.plot(b1,post,"o")
少し気味の悪いグラフができましたが、別に失敗ではありません。(多分)
これはあくまでパラメータが二つであることに起因するもので、いわば同時確率分布を無理やり二次元平面上に書き出しているから気味があるいだけです。
本来なら積分することで周辺密度関数にする必要がありますが離散化されているので方法がわかりません。
ただここからでも十分パラメータの統計量はわかるので十分としておきましょう。
今後ももう少しベイズに詳しくなったらちょくちょく追記しています。