Pythonで実装しながら緑本を学ぶ (第8章マルコフ連鎖モンテカルロ(MCMC)法とベイズ統計モデル)

データ解析のための統計モデリング入門(通称、緑本)を読み進めています。
述べられている理論を整理しつつ、Rでの実装をPythonに置き換えた際のポイントなども深掘りしていきます。

今回は第8章です。実装は以下で公開しています。

introduction_to_machine_learning_with_python/chapter8.ipynb at master · ohke/introduction_to_machine_learning_with_python · GitHub

8 マルコフ連鎖モンテカルロ(MCMC)法とベイズ統計モデル

第7章では、ランダム効果を組み込んだ最尤推定を行いましたが、発生源が増えるとその分だけ多重積分が必要となり、計算が困難になります。

この章では、観測されたデータからある確率分布に従うランダムサンプルを取得するマルコフ連鎖モンテカルロ法(MCMC)と、パラメータを確率分布で表現するベイズ統計学を導入することで、上の問題を解くための準備をします。

8.1 例題：種子の生存確率(個体差なし)

架空植物20個体の種子8個の生死を調べることを想定する(第6章と同じ)。

(生存)種子数 $y_i$ が二項分布に従うと仮定すると、ある個体iの種子数が $y_i$ の確率と、尤度L(q)および対数尤度 $\log L(q)$ は以下の式で計算されます。

$$ p(y_i \mid q)={}_8\mathrm{C}_{y_i}q^{y_i}(1-q)^{8-y_i} $$

$$ L(q)=p({\bf x} \mid q)=\prod_i p(y_i \mid q) $$

$$ \log L(q)=\sum_i {y_i \log q + (8-y_i) \log(1-q)}+\mbox{定数} $$

$\frac{d \log L(q)}{dq}=0$ となるqを求めると、 $\hat{q}=0.46$ と推定されます。なおこのデータは $q=0.45$ で生成されています。

f:id:ohke:20180223092916p:plain

8.2 ふらふら試行錯誤による最尤推定

上のように解析的に $\hat{q}$ を求められない場合に、最尤推定を行なう方法を検討します。

非効率ですが1つの方法として、qを増減させながら、対数尤度が最大となる $\hat{q}$ を逐次的に探索するアルゴリズムが考えられます。

qを離散化する
- 例えば、0.01から0.99まで0.01刻みの値とする
スタート地点として、適当なqを選択して、その地点の対数尤度を求める
- q=0.30において-46.38
qの両隣の点からランダムに選択して、とする
- $q_{new}=0.29$ または $q_{new}=0.31$ が選択される(それぞれ選択される確率は0.5)
における対数尤度を求める
- $q_{new}=0.29$ で-47.62、 $q_{new}=0.31$ で-45.24
の方の対数尤度が大きければ、qをで更新する
- $q_{new}=0.29$ なら更新せず、 $q_{new}=0.31$ ならqを0.31で更新する
試行回数(例えば100回)だけ、3〜5を繰り返す

6.の指定回数を十分大きくすれば、やがて対数尤度が最大となる $\hat{q}$ で収束することがイメージできるかと思います。

先程のデータに対して、上のアルゴリズムで最尤推定した結果が以下となります。
スタート地点が0.3でも0.6でも、 $\hat{q}=0.45$ で安定していることがわかります。

f:id:ohke:20180223152430p:plain

8.3 MCMC アルゴリズムのひとつ：メトロポリス法

8.2の最尤推定アルゴリズムの手順5.に、以下の追加ルールを設けることで、MCMCのメトロポリス法となります。

$q_{new}$ の対数尤度の方が小さい場合でも、確率rでqを $q_{new}$ で更新する

この確率rは、qと $q_{new}$ の尤度比です。

$$ r=\frac{L(q_{new})}{L(q)}=\exp(\log L(q_{new}) - \log L(q)) $$

この変更により、対数尤度が悪くなる場合であっても、確率rでqが更新されるようになります(対数尤度が良くなる場合は、常に $q_{new}$ で更新されます)。
Pythonでは以下のような実装となります。

# 対数尤度
def loglikelihood(data, q):
    ll = 0
    
    for i, r in data.iterrows():
        ll = ll + math.log(scipy.misc.comb(r['N'], r['y'])) + r['y']*math.log(q) + (r['N'] - r['y'])*math.log(1 - q)
        
    return ll

# MCMC(メトロポリス法)
def mcmc_metropolis(data, q_start, number_of_samples):
    q_current = q_start
    ll_current = loglikelihood(data, q_current)
    
    q = [q_current]
    ll = [ll_current]
    
    for r1, r2 in zip(np.random.random(number_of_samples), np.random.random(number_of_samles)):
        q_new = q_current + 0.01 if r1 > 0.5 else q_current - 0.01
        if q_new <= 0.01:
            q_new = 0.02
        elif q_new >= 0.99:
            q_new = 0.98
        ll_new = loglikelihood(data, q_new)
        
        # 対数尤度が悪くなる場合でも、尤度比の確率でqを更新
        if ll_current < ll_new or (math.exp(ll_new - ll_current)  > r2):
            q_current = q_new
            ll_current = ll_new
            
        q.append(q_current)
        ll.append(ll_current)
    
    return q, ll

# サンプル数100のMCMCでqと対数尤度の遷移を計算
q_100, ll_100 = mcmc_metropolis(data, 0.3, 100)

試行回数(MCMCではサンプリング数と呼ぶ)を100、1000、100000とした時の、それぞれのqの遷移を見てみます。 100回程度では最尤推定値までまだまだ遠く、また最尤推定値に到達後も変動することがわかります。

f:id:ohke:20180223152515p:plain
f:id:ohke:20180223152527p:plain
f:id:ohke:20180223152545p:plain

qの度数からqの確率分布が得られます。サンプリング数を十分大きくし、増やしても変動しなくなった時の確率分布を、マルコフ連鎖の定常分布( $p(q \mid {\bf Y}$ )とする)と呼びます。
サンプル1000回と100000回の確率分布を比較します。1000回の方が偏りが大きく、またスタート地点( $q=0.3$ )に引きずられてピークが左側にあることがわかります。

f:id:ohke:20180223230620p:plain

どのようなqからスタートしても最終的には定常分布 $p(q \mid {\bf Y})$ に従うが、サンプリング数は十分大きくする必要があります。ナイーブなメトロポリス法にはいくつか改善ポイントがあります。

あるステップとその次のステップでサンプルされる値の相関が低いアルゴリズムを選ぶ
(でたらめに選んだスタート地点の影響が強い)サンプリングの最初部分の値を捨てる
異なるスタート地点を選んだ複数のサンプリングを足し合わせる

qが0.1, 0.3, 0.6, 0.9からスタートした場合の動きを見てますと、いずれの場合もサンプル数が増えると0.45に近づくことがわかります。

f:id:ohke:20180223231235p:plain

定常分布 $p(q \mid {\bf Y})$ は尤度 $L(q)$ に比例する確率分布です( $\sum_q L(q)$ は定数と見るため)。つまり、十分に長いMCMCサンプルは、メトロポリス法の定常分布 $p(q \mid {\bf Y}$ からのランダムサンプルとなります。

$$ p(q \mid {\bf Y})=\frac{L(q)}{\sum_q L(q)} $$

8.4 MCMCサンプリングとベイズ統計モデル

これまでGLMで推定してきたパラメータは、ただ1つに定まる値でした(統計学の枠組みでは頻度主義と呼ばれます)。

一方で、(上で求めたqのように)パラメータはある確率分布に従うものと表現されるものは、ベイズ統計学と呼ばれます。

上の例題をベイズ統計学で捉え直します。ベイズ統計学では、ベイズの公式で推論する統計学です。

$$ p(q \mid {\bf Y})=\frac{p({\bf Y} \mid q)p(q)}{\sum_q p({\bf Y} \mid q)p(q)} $$

$p(q \mid {\bf Y})$ は、データ ${\bf Y}$ が得られた時にqが従う確率分布(事後分布)
は、qの値の時にデータが観測される確率
- 二項分布の積である尤度が対応するので、 $p({\bf Y} \mid q)=L(q)$
はデータが無い時のqの確率分布(事前分布)
- 離散一様分布だろうか？(後の章の問題とする)
$\sum_q p({\bf Y} \mid q)p(q)$ は、条件付き確率の和なので、 $p({\bf Y})=\sum_q p({\bf Y} \mid q)p(q)$ で、qが不明の時にデータ ${\bf Y}$ が得られる確率(定数)

$p({\bf Y} \mid q)=L(q)$ なので、ベイズの公式を以下のように変形できます。もし事前分布 $p(q)$ がqによらない定数なら、メトロポリス法で得られた $p(q \mid {\bf Y}=\frac{L(q)}{\sum_q L(q)}$ と一致することがわかります。

$$ p(q \mid {\bf Y})=\frac{L(q)p(q)}{\sum_q L(q)p(q)} $$

8.5 補足説明

8.5.1 メトロポリス法と定常分布の関係

メトロポリス法で得られたqのサンプルが、定常分布 $p(q \mid {\bf Y})$ からのランダムサンプルであるためには、2つの条件を満たす必要がある。

qが任意の初期値から定常分布 $p(q \mid {\bf Y})$ に収束すること
あるqが $p(q \mid {\bf Y})$ に従っていて、次に得られた $q_{new}$ も $p(q \mid {\bf Y})$ に従っていること

8.5.2 ベイズの定理

先に登場した下式は、ベイズの定理と呼ばれます。

$$ p(q \mid {\bf Y})=\frac{p({\bf Y} \mid q)p(q)}{\sum_q p({\bf Y} \mid q)p(q)} $$

元々は条件付き確率と同時確率の関係を整理した、 $p(A \mid B)p(B)=p(A, B)$ という定義です。
$p({\bf Y}, q)=p({\bf Y} \mid q)p(q)$ 、 $p({\bf Y})=\sum_q p({\bf Y} \mid q)p(q)$ の関係を使うことで、上の式を導出しています。

$$ p(q \mid {\bf Y})=\frac{p(q, {\bf Y})}{p({\bf Y})}=\frac{p({\bf Y} \mid q)p(q)}{p({\bf Y})}=\frac{p({\bf Y} \mid q)p(q)}{\sum_q p({\bf Y} \mid q)p(q)} $$