Pythonによるマルコフ連鎖モンテカルロシミュレーション

　モデル関数のパラメータの最適解の推定方法は、モデルとサンプルの残差を最小にするようにモデルパラメーターを決定する最小二乗法や、ある確率分布からサンプルを生み出すと考えられる尤もらしさ(尤度関数)を最大化するようにモデルパラメーターを決定する最尤法などあるが、今回はマルコフ連鎖モンテカルロ法(MCMC : Markov chain Monte Carlo methods)によるサンプリング手法を用いて、サンプルとモデルの差により得られる尤度関数を最適化することにより、局所最適解を探索する回帰分析の手法を紹介する。

MCMC法(Markov chain Monte Carlo methods)とは

　MCMCは多次元の確率密度関数を事前分布として持つ標準集団に、詳細つり合い条件が成り立つことを仮定してマルコフ連鎖を繰り返すことで、サンプリングを行う手法である。

　ベイズ統計の基本概念として、確率は古典統計学のような頻度やモデルに基づく固定値ではなく、新しい情報が収集されると変化する動的な確率であるとしている。式にすると以下のようになる。

P(x_{t-1} |x_t ) = \frac{P(x_t|x_{t-1})Q(x_{t-1})}{Z} \\
Z= \sum_x P(x_t|x_{t-1})Q(x_{t-1})

　これはベイズの定理といい、人間が主観的に設定した事前分布と実際のデータによって得られる尤度の加重平均によって、新たに事後分布が得られるという条件である。ここで$P(xt-1 | xt)$ は事象x-1の後に事象xが起こる条件付き確率分布で、事後分布と呼ばれる。$P(xt | xt-1) $は$xt-1$ の値が決まっているときに事象$xt$が観測される確率分布で、提案分布と呼ばれる。$Q(xt-1)$ は$xt-1$が起こりうる確率分布、言い換えると$xt-1$のデータが観測されるかくりつで、事前分布と呼ばれる。分母 $Z$ は確率の総和を 1 とするための規格化定数である。

                            事後分布 = 尤度 × 事前分布

　またMCMCでは確率変数$x_t$の条件付き分布が、一つ前の確率変数$xt−1$の値にのみ依存するとしている。言い換えると未来の状態は、ひとつ前(もしくは現在の状態)より過去の状態には依存しないとする。これをマルコフ連鎖と呼ぶ。

 P(x_t∣x_{t−1},x_{t−2},…)=P(x_t  ∣x_{t−1} )

　古典統計学における最尤法では、尤度が最も大きくなるようモデルを決定するのに対して、ベイズ統計では尤度と事前分布から計算される事後分布が最大になるようにモデルを決定する。
　メトロポリス・ヘイスティングス法(M-H法)では，重み関数$w(x | y)$を与えて以下の条件が成り立つとしている。

P(x|y) = \frac{P(y|x)P(x)}{Z} w(x | y)\\
Z= \sum_x P(y|x)P(x)

重み関数w(x | y)は、提案分布P(x | y)により求められた事後分布をw(x | y)の確率で採択することを意味しており、これにより詳細つり合い条件が満たされるとしている。提案分布がガウス分布のような対象の分布を用いる場合は，単にメトロポリス法と呼ぶ。今回は提案分布に多次元正規分布を用いたメトロポリス法を用いた。

サンプルコード

　今回は挙動を理解するためmcmcのpython moduleを使わずに書いてみた。以下は、y=2x+2の関数に適当な乱数を振りサンプルを作成し、初期値をかなり離れたa=10000,b=-10000に設定して、メトロポリス法により尤度関数が最大になるaとbを求めている。最終的な出力はaとbの尤度分布と尤度関数の対数値である。

import sys
import os
import math
import numpy as np
import random

class Func:
    def __call__(self,x,a,b):
        y=a*x+b
        return y

class Likelihood:
    def __init__(self,par1,par2):
        self.a=par1
        self.b=par2
    def __call__(self,x,y,sigma):
        Func_instance=Func()
        Log_Likelihood=-(y-Func_instance(x,self.a,self.b))**2/(2.0*sigma**2)
        return Log_Likelihood

a0,b0=2.0, 2.0 
xdata=np.zeros(1000)
ydata=np.zeros(1000)
sigmadata=np.zeros(1000)
for i in range(1000):
    tmp1=random.uniform(0,100)
    Make_Data=Func()
    xdata[i]=tmp1
    ydata[i]=Make_Data(tmp1,a0,b0)+np.random.normal(loc=0,scale=10)
    sigmadata[i]=abs(np.random.normal(loc=0,scale=1))

a,b=10000,-10000
FA = 0.0
cov0=[[0.1,0.0],[0.0,0.1]]
p=1
S0=np.random.multivariate_normal([0.0,0.0], p*cov0,1) 
m0,X0=np.array([a0,b0]),np.array([a0,b0])
Xn=X0 + S0
pprev,mprev,covprev= p,m0,cov0
Pre_Sum_Log_Likelihood = -1e100
FinalLikelihood = 0
ok=0.0
i=1

Nboot=10000 
para1=np.zeros(Nboot)
para2=np.zeros(Nboot)
Larray=np.zeros(Nboot)
while i < Nboot:
    mn = mprev+(10.0/(i+100.0))*(Xn-mprev)
    covn = covprev+(10.0/(i+100.0))*((Xn-mprev)*(Xn-mprev).transpose()-covprev)
    pn = pprev+(10.0/(i+100.0))*(FA-0.243)
    S0 = np.random.multivariate_normal([0.0,0.0], pn*covn,1)
    a_tmp=a+S0[0][0]
    b_tmp=b+S0[0][1]
    tmp=Likelihood(a_tmp,b_tmp) 
    All_Log_Likelihood=0.0
    All_Log_Likelihood=[tmp(xdata[j],ydata[j],sigmadata[j]) for j in range(np.size(xdata))]
    Sum_Log_Likelihood=np.sum(All_Log_Likelihood)
    Log_Likelihood_Ratio=Sum_Log_Likelihood-Pre_Sum_Log_Likelihood
    if Log_Likelihood_Ratio >= 0.0:
        a=a_tmp
        b=b_tmp
        FA=1.0
        FinalLikelihood=Sum_Log_Likelihood
        Pre_Sum_Log_Likelihood=Sum_Log_Likelihood
    else:
        pp = np.random.rand()*0.01
        if Log_Likelihood_Ratio >= math.log(pp):
            a=a_tmp
            b=b_tmp
            FA=1.0
            FinalLikelihood=Sum_Log_Likelihood
            Pre_Sum_Log_Likelihood=Sum_Log_Likelihood
        else:
            FinalLikelihood=Pre_Sum_Log_Likelihood
            FA = 0
    pprev = pn
    mprev = mn
    covprev = covn
    Xn = np.array([a,b])  
    para1[i]=a 
    para2[i]=b
    Larray[i]=FinalLikelihood 
    i += 1 

outputfile="./output.dat"
label="a,b,LogLikelihood"
line=np.array([para1,para2,Larray])
np.savetxt(outputfile, line.T, fmt='%0.8e',delimiter=',', newline='\n', header=label, comments="")

結果

　傾きaは2でサンプル生成したのに対して、1.94付近が分布の頂点を示した。切片bは真値2に対して、5.72付近が分布の頂点を示した。ちょっとbが離れすぎてる気もするので、修正の余地はかなりありそう。。。

参考文献

上記の記事を書くため、以下の素晴らしくまとまったサイトを参考にしました。
はじめてのMCMC (メトロポリス・ヘイスティングス法)
誰でもわかるマルコフ連鎖モンテカルロ法（MCMC）入門
 EMアルゴリズム徹底解説

余談

MCMC法は、EMアルゴリズム(期待値最大化法)におけるEステップの積分計算を、モンテカルロ法を用いてかつマルコフ性を仮定してサンプリングする方法である。