Pythonで実装しながら緑本を学ぶ (第9章 GLMのベイズモデル化と事後分布の推定)

Python 統計学 読書メモ

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

データ解析のための統計モデリング入門――一般化線形モデル・階層ベイズモデル・MCMC (確率と情報の科学)

今回は第9章です。PyMC3を使って、GLMをベイズモデルで表現します。実装は以下で公開しています。

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

9 GLMのベイズモデル化と事後分布の推定

9.1 例題：種子数のポアソン回帰(個体差なし)

架空植物20個体において体サイズと種子数を計測した。個体iにおける体サイズを x_i 、種子数を y_i とします。得られたデータから体サイズ x_i と種子数 y_i の関係について調べる、というのが今回の例題です。

計測結果は以下にプロットしています。
なお著者サイトで提供されているサンプルデータはRDataですので、Pythonで扱えるようにするため、PypeRでロード後にDataFrame化しています。PypeRについては以前の投稿も参考にしてください。

f:id:ohke:20180303113433p:plain

種子数は上限が無い離散値ですので、ポアソン分布でばらつきを表現できそうです。
statsmodelsで最尤推定すると、ポアソン分布の平均 \lambda_i は下式で得られます(なお、このデータは作成するときには \lambda_i=\exp(1.5+0.1x_i) から生成されたとのことです)。

$$ \lambda_i = \exp(1.5661 + 0.0833x_i) $$

import statsmodels.formula.api as smf
result = smf.poisson('y ~ x', data=data).fit()
result.summary()

f:id:ohke:20180303120036p:plain

9.2 GLMのベイズモデル化

ベイズモデル化したとしても、中核はポアソン回帰のGLM。

  • 平均 \lambda_i のポアソン分布 p(y_i \mid \lambda_i) に従う
  • 線形予測子と対数リンクを使い、 \lambda_i =\exp(\beta_1 + \beta_2 x_i) で指定する

このモデルでの尤度関数Lは以下となります。 x_i を定数としている(Lのパラメータになっていない)ことに注意してください。

$$ L(\beta_1, \beta_2)=\prod_i p(y_i \mid \lambda_i)=\prod_i p(y_i \mid \beta_1, \beta_2, x_i) $$

ある \beta_1, \beta_2 において {\bf Y} ( {\bf Y}= \{ y_i \} )が得られる確率は、尤度関数と一致します(この関係は8.4で導入・説明しています)。

$$ p({\bf Y} \mid \beta_1, \beta_2)=L(\beta_1, \beta_2) $$

ベイズモデルの事後分布は、尤度×事前分布に比例するので、以下の関係が成り立ちます(こちらも8.4を参考にしてください)。

$$ p(\beta_1, \beta_2 \mid {\bf Y}) \propto p({\bf Y} \mid \beta_1, \beta_2)p(\beta_1)p(\beta_2) $$

9.3 無情報事前分布

まずは事前分布 p(\beta_1) と p(\beta_2) (まとめて p(\beta_{*}) )を設定する。

データ {\bf Y} が得られていない状態で決める事前分布なので、無情報事前分布と呼ばれます。つまり、どんな値( [-\infty, +\infty ] )でも良いのです。
こうした分布はの生成方法は2つあります。

  • 広い範囲 (例えば[-100000, +100000])の一様分布
  • 平均0で標準偏差が大きい(平べったい)正規分布

分散100の正規分布は以下のようにプロットできます。今回は標準偏差100の正規分布を、事前分として使います。

f:id:ohke:20180303131934p:plain

9.4 ベイズ統計モデルの事後分布の推定

事前分布が定まったので、事後分布 p(\beta_1, \beta_2 \mid {\bf Y}) をMCMCサンプリングで推定します。

書籍では、WinBUGS+R2WinBUGSで、RからMCMCサンプリンによる推定が行われています。
ベイズモデルでパラメータ推定できるPythonパッケージとして、PyStanやPyMCが有名です。今回はPyMC3を使います。いつも通り、pip install pymc3でインストールしておきます。

http://docs.pymc.io/

PyMC3によるGLMのパラメータ推定は、3ステップで行います。

  1. モデルを定義する
  2. サンプリングして推定する
  3. 得られた結果を確認する
import pymc3 as pm

# モデルを定義する
with pymc3.Model() as model:
    # 事前分布をN(0, 100)の正規分布で設定
    beta1 = pymc3.Normal('beta1', mu=0, sd=100)
    beta2 = pymc3.Normal('beta2', mu=0, sd=100)
    
    # 線形予測子θをβ1+β2xで設定
    theta = beta1 + beta2*data['x'].values
    
    # ログリンク関数(log(μ)=θ⇔μ=exp(θ))を設定し、ポアソン分布で推定する
    y = pymc3.Poisson('y', mu=np.exp(theta), observed=data['y'].values)

# サンプリングして推定する
with model:
    # 101個目から3個置きでサンプルを取得するチェインを3つ作る
    # NUTSではburnとthinが効いていない？
    trace = pymc3.sample(1600, burn=100, thin=100, njobs=3, random_seed=0)

# 得られた結果を確認する
pymc3.traceplot(trace) # サンプリング過程を表示する
pymc3.summary(trace) # 推定結果を表示する

モデルを定義する

まずwith句でモデルクラス(pymc3.Model)を作成し、事前分布、線形予測子、リンク関数、確率分布を指定します。

  • \beta_1, \beta_2 の事前分布には、それぞれ平均0・標準偏差100の正規分布(pymc3.Normal)を使う
  • 線形予測子は \theta=\beta_1 + \beta_2 x
  • ポアソン分布(pymc3.Poisson)で尤度を計算
    • ログリンク関数 \log(\mu)=\theta \Leftrightarrow \mu=\exp(\theta) で平均 \mu を設定
import pymc3

# モデルを定義する
with pymc3.Model() as model:
    # 事前分布をN(0, 100)の正規分布で設定
    beta1 = pymc3.Normal('beta1', mu=0, sd=100)
    beta2 = pymc3.Normal('beta2', mu=0, sd=100)
    
    # 線形予測子θをβ1+β2xで設定
    theta = beta1 + beta2*data['x'].values
    
    # ログリンク関数(log(μ)=θ⇔μ=exp(θ))を設定し、ポアソン分布で推定する
    y = pymc3.Poisson('y', mu=np.exp(theta), observed=data['y'].values)

サンプリングしてパラメータ推定する

次に、pymc3.sampleメソッドでサンプラーの定義とサンプリングを行い、パラメータを推定する。

  • 最初の引数(draws)で、サンプル数を指定
  • stepで、サンプリングアルゴリズムを指定
    • Metropolis、HamiltonianMC、NUTSなどが選択できます(デフォルトはNUTS)
  • tuneで、先頭から捨てるサンプル数を指定する(WinBUGSではn.burninで指定する値)
    • サンプルの最初の方は、ランダムに選ばれた初期値の影響を大きく受けるため、捨てた方が良い
  • njobsで、チェイン数(サンプル列数)を指定
    • 3を指定すると、3つの異なる初期値からそれぞれサンプリングが行われる
  • observedに、観測された従属変数(ここではy)を渡します
  • 2個飛ばしでサンプリング(つまり合計500個)するために、スライス表記[::3]でサンプリング過程を取得
# ハミルトニアンモンテカルロ法
with model:
    # 101個目から3個置きでサンプルを取得するチェインを3つ作る
    trace = pymc3.sample(1500, step=pymc3.HamiltonianMC(), tune=100, njobs=3, random_seed=0)[::3]

得られたサンプリング過程(trace)は、添字でアクセスできます。

print('Trace type:', type(trace)) # Trace type: <class 'pymc3.backends.base.MultiTrace'>
print('Trace length:', len(trace)) # Trace length: 500
print('trace[0]:', trace[0]) # trace[0]: {'beta1': 2.0772965015391716, 'beta2': -0.02971672503615687}

得られた結果を確認する

サンプリング後には、初期値やサンプリング数などの設定値が適切であったかどうかを確認する必要があります。

pymc3.traceplotメソッドで、サンプリング過程がグラフ化できます。

pymc3.traceplot(trace)

メトロポリス法(上図)とハミルトニアンモンテカルロ法(下図、以降はHMCと表記)で並べて比較してみます。

右図のパラメータの推移を見ると、HMCではサンプル列同士が近づきつつありますがまだ不安定な状態にある(つまりサンプリング数が不足している)ことがわかります。対して、HMCではいずれのサンプル列同士も十分近づき、類似した波形となっています。
いずれも、3つのサンプル列がプロットされていることや、サンプル数が500(1500から2つ飛ばしでサンプリングしているため)となっていることに注意してください。

  • メトロポリス法
    f:id:ohke:20180303231830p:plain

  • HMC
    f:id:ohke:20180303231728p:plain

また、pymc3.summaryメソッドで、各パラメータの推定値とそれぞれの統計値を確認できます。

pymc3.summary(trace)

こちらもメトロポリス法(上表)とHMC(下表)で並べてみます。
各項目の詳細は次節で見ていきますが、このうちRhatで表記されている \hat{R} 指数はサンプル列間のばらつきを表す値で、パラメータ毎に求められます。この値が1に近いほどサンプル列間のばらつきよりも列内のばらつきが大きくなるので、収束していると言えます。経験的には \hat{R} \le 1.1 が1つの目安ですが、 \hat{R} も十分なサンプル数でない場合は安定した結果が得られないので注意が必要です。
メトロポリス法では1.007、HMCでは0.999なので、後者の方が僅かですがより収束していると言えそうです。

  • メトロポリス法
    f:id:ohke:20180303232045p:plain

  • HMC
    f:id:ohke:20180303232202p:plain

9.5 MCMCサンプルから事後分布を推定

得られた推定結果から事後分布を確認します。

HMCで得られたサンプリング過程を再掲します。
左図は、 \beta_1 と \beta_2 の周辺事後分布で、カーネル密度推定で近似された確率密度関数で表現されています。周辺事後分布は、あるパラメータ1つに関する事後分布で、ここでは p(\beta_1 \mid {\bf Y}) と p(\beta_2 \mid {\bf Y}) となります。

f:id:ohke:20180304125331p:plain

周辺事後分布に対して、 p(\beta_1, \beta_2 \mid {\bf Y}) は同時事後分布と呼ばれます。
パラメータの組み合わせを散布図にプロットします。 \beta_1 と \beta_2 の相関がかなり強いサンプリングが行われていることがわかります。本書の傾向と大きく異なっており、サンプリングアルゴリズム(本書ではWinBUGSのギブスサンプリング実装)の違いによるものと思われます。

f:id:ohke:20180304141715p:plain

なお、サンプリングした値はtrace(MultiTraceクラス)からget_valuesメソッド(返り値の型はnumpy.ndarray)で取得できます。

  • 1番目の引数(varname)にパラメータ名を指定
  • chainsオプションに、サンプル列のインデックスを渡す
# サンプル列数分だけ繰り返す
for i in trace.chains:
    # 各サンプル列のパラメータの平均値を計算
    beta1_averages += trace.get_values('beta1', chains=i) / trace.nchains
    beta2_averages += trace.get_values('beta2', chains=i) / trace.nchains

次に、HMCの統計量から事後分布の推定値を確認します。

  • パラメータ \beta_1 は平均1.5599で、95%信用区間は0.8630〜2.2622
    • \beta_1 は95%の事後確率で0.8630〜2.2622に収まる、と解釈できる
  • パラメータ \beta_2 は平均0.0826で、95%信用区間は0.0017〜0.2145
  • n_effは有効なサンプルサイズで、サンプル間の相関が高いと、この値が小さくなります

f:id:ohke:20180303232202p:plain

9.6 複数パラメーターのMCMCサンプリング

8章では、1パラメータについてメトロポリス法でMCMCサンプリングする方法を説明しました。

メトロポリス法には、更新前と更新後の値の相関が強く、なかなか収束しないという問題がありました。
この問題を解決するサンプリング方法の1つにギブスサンプリングがあります。ギブスサンプリングは、「新しい値の確率分布を作ってその確率分布からランダムに選択する」という方法で値を更新することで、更新前後の値の相関を弱くします。

さらに今回の例題のように、複数のパラメータ( \beta_1, \beta_2 )のMCMCサンプリングを考える必要があります。
こうした場合、全てのパラメータを同時に更新するよりも、 \beta_1 と \beta_2 を交互に少しずつ更新していく方が簡単です。

2パラメータのギブスサンプリングをまとめると以下のアルゴリズムとなります。

  • (1) {\beta_1, \beta_2} の適当な初期値を設定する
    • 例えば {\beta_1, \beta_2}={1.5, 0} と置く
  • (2) p(\beta_1 \mid {\bf Y}, \beta_2) に従う乱数を発生させ、得られた値を新しい \beta_1 とする
    • 他の変量を全て定数とする一変量確率分布で、全条件付き分布(FCD)と呼ばれる
    • FCDに従う乱数1つを発生させる( \beta_1^{new}=2.052 が得られたとする)
      • FCDからのサンプリングする方法は理解できませんでした

$$ p(\beta_1 \mid {\bf Y}, \beta_2=0.0) \propto \prod_i \frac{\lambda_i^{y_i} \exp(-\lambda_i)}{y_i !}p(\beta_1) $$

  • (3) p(\beta_2 \mid {\bf Y}, \beta_1) に従う乱数を発生させ、得られた値を新しい \beta_2 とする (2の逆)

$$ p(\beta_1 \mid {\bf Y}, \beta_1=2.052) \propto \prod_i \frac{\lambda_i^{y_i} \exp(-\lambda_i)}{y_i !}p(\beta_2) $$

  • (4) この新しい {\beta_1, \beta_2} を記録する
  • (5) 十分なサンプル数が得られるまで(2)〜(4)を繰り返す

Python: foliumでJupyter Notebookに地図を描画する

Python

Jupyter Notebook上で、緯度経度の情報を地図へ簡単にプロットできる方法を探していたところ、foliumの使い勝手が良かったので紹介します。

folium

PythonからLeaflet.jsで地図をプロットするパッケージです。

https://github.com/python-visualization/folium

使用前にpip install foliumまたはconda install foliumでインストールします。

地図を表示する

まずは、千葉県千葉市の地図を表示します。

Mapで生成できます。

  • locationオプションに中心の緯度・経度を渡しています
  • Mapインスタンスを出力すると、leaflet.jsが読み込まれて表示されます
  • デフォルトではOpenStreetMapの地図が使われます(tilesオプションでMapboxやStamenにも変更できます)
  • zoom_startで表示する時のズームを変えられます(デフォルト10で、大きくするとより拡大されます)
import folium

chiba_map = folium.Map(location=[35.607451, 140.106340])

chiba_map

f:id:ohke:20180228092918p:plain

マーカをプロットする

次に千葉県の4つの市にマークします(緯度経度は市役所所在地、人口は千葉県のサイトから2018/1/1時点のデータを参照)。

市 緯度 経度 人口
千葉市 35.607451 140.106340 975,535
館山市 34.996596 139.869906 46,349
銚子市 35.734795 140.826926 61,674
浦安市 35.653146 139.902058 168,169

Markerで地図上にマーカをプロットできます。

  • popupオプションで、マーカをクリックした時表示されるテキストを設定できます
  • add_toメソッドでマップ上にプロットします
import pandas as pd

chiba_cities = pd.DataFrame({
    'city': ['千葉市', '館山市', '銚子市', '浦安市'],
    'latitude': [35.607451, 34.996596, 35.734795, 35.653146],
    'longtude': [140.106340, 139.869906, 140.826926, 139.902058],
    'population': [975535, 46349, 61674, 168169]
})

chiba_map = folium.Map(location=[35.607451, 140.106340], zoom_start=9)

for i, r in chiba_cities.iterrows():
    folium.Marker(location=[r['latitude'], r['longtude']], popup=r['city']).add_to(chiba_map)
    
chiba_map

f:id:ohke:20180228093011p:plain

マーカはいくつか種類があります。

CircleMarkerを使って、人口比を可視化してみます。

f:id:ohke:20180228093047p:plain

緯度経度を含む情報をぱぱっと可視化したいときに便利でしたので紹介しました。

GeoJSONを描画することもできます。詳しくはこちらの方の記事が参考になるかと思います。

国土数値情報とfoliumで市区町村区切りのコロプレス図を描く

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

Python 統計学 読書メモ

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

データ解析のための統計モデリング入門――一般化線形モデル・階層ベイズモデル・MCMC (確率と情報の科学)

今回は第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} を逐次的に探索するアルゴリズムが考えられます。

  1. qを離散化する
    • 例えば、0.01から0.99まで0.01刻みの値とする
  2. スタート地点として、適当なqを選択して、その地点の対数尤度を求める
    • q=0.30において-46.38
  3. qの両隣の点からランダムに選択して、 q_{new} とする
    • q_{new}=0.29 または q_{new}=0.31 が選択される(それぞれ選択される確率は0.5)
  4. q_{new} における対数尤度を求める
    • q_{new}=0.29 で-47.62、 q_{new}=0.31 で-45.24
  5. q_{new} の方の対数尤度が大きければ、qを q_{new} で更新する
    • q_{new}=0.29 なら更新せず、 q_{new}=0.31 ならqを0.31で更新する
  6. 試行回数(例えば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が従う確率分布(事後分布)
  • p({\bf Y} \mid q) は、qの値の時にデータ {\bf Y} が観測される確率
    • 二項分布の積である尤度が対応するので、 p({\bf Y} \mid q)=L(q)
  • p(q) はデータ {\bf Y} が無い時の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)} $$