メトロポリス・ヘイスティング法をPythonで簡単に実装してみる
メトロポリス・ヘイスティングス法とは?
- ある特定の確率分布から、その確率密度の大小にしたがってサンプル値を抽出するための基本的なアルゴリズムの一種
- サンプリング自体は、ベイズ統計で、事前分布と尤度の掛け合わせからなる事後分布からサンプル抽出を行いたい時によく用いられる。
- しかし、一般にメトロポリス・ヘイスティングス法によるサンプリングは収束に時間がかかるという課題があるため、改良版として、ハミルトニアンモンテカルロ(HMC)という改良版がある。しかし、基礎は共通。
サンプリング手法は沢山あるが、メトロポリス・ヘイスティング法は諸々の手法のベースになっているので、抑えておく必要がある。
ステップ数が無限ならば、正確に対象の確率分布をトレースしたサンプリングを実施できる事は理論的に証明されているらしいです。
アルゴリズム のアイデア
直感的に書くと以下のような仕組みでサンプリングを行う事で、対象とする確率分布に従うサンプル群を得ることができます。
対象とする確率変数を とし、確率分布を
とおきます。また、サンプリング結果を
と表します。(例えば、2番目のサンプリング結果は
と表す事にします。)
STEP1: 初期値を決め、確率密度の値
を求めておきます
STEP2: に、ランダムな値を足し、遷移先の値候補
及び、確率密度の値
を求めておきます。
STEP3: ならば、
とし
をサンプルとして採用(受容と言います。) 一方で、
ならば、確率密度の比率
の確率で受容しますが、それ以外の場合は
として、既存の値をサンプル値として採用します。
STEP4: STEP2に戻り、決められたステップ数の数だけSTEP2->3を繰り返す。
現在位置よりも、移動先の確率密度が低い場合であっても、一定確率で受容する事で、確率分布に沿ったサンプリングを実行できるのが肝です。より詳しい事を知りたい場合は、以下の書籍や動画が理解の助けになるかと思います。
Pythonによるベイズ統計モデリング: PyMCでのデータ分析実践ガイド
- 作者:マーティン,オズワルド
- 発売日: 2018/06/22
- メディア: 単行本
データ解析のための統計モデリング入門――一般化線形モデル・階層ベイズモデル・MCMC (確率と情報の科学)
- 作者:久保 拓弥
- 発売日: 2012/05/19
- メディア: 単行本
Pythonで実装してみる
サンプリングを実行する関数metropolis()の定義
実装したコードは以下になります。なお、本コードは先に紹介している書籍の実装を参考にしています。というよりかは、ほぼそのままなので、是非本書をご覧になってください。
def metropolis(func, steps = 10000):
# サンプリング結果保存用のリスト
samples = np.zeros(steps)
# STEP1: 初期値の設定 及び 確率密度の計算
# 与えられた確率分布の期待値で初期値に設定
old_x = func.mean()
old_prob = func.pdf(old_x)
# stepsの数だけループ
for i in range(steps):
# STEP2: 正規分布に従う乱数を足して、新しい遷移先候補を決定
new_x = old_x + np.random.normal(0, 0.5)
# 遷移先候補の確率密度を計算
new_prob = func.pdf(new_x)
# STEP3: 確率密度の比を計算
acceptance = new_prob / old_prob
# acceptanceを基に受容 or 棄却を決定
if acceptance >= np.random.random():
samples[i] = new_x
old_x = new_x
old_prob = new_prob
else:
samples[i] = old_x
return samples
文章で書くとややこしそうでしたが、とても短いコードでサンプリングを実行する関数を記述できました。acceptanceの判定に、np.random.random()で[0~1]の乱数を用いている点がポイントかと思います。
STEP2では、移動先の提案に正規分布に従う乱数を用いていますが、この場合は、正規分布を提案分布として採用している事になります。提案分布の選択は、幾分かサンプリング結果にも影響を与える事になります。
サンプリングの実行
サンプリング対象となる確率分布を設定 -> 関数に渡してあげます。
func = stats.beta(1.2,3) samples = metropolis(func, steps = 10000)
今回の例では、
のベータ分布からのサンプリングを実行しています。stats.~で定義する確率分布を正規分布や、ポアソン分布やその他の分布に変える事で、同様にサンプリングを行う事ができるでしょう。
次は?
より効率良くサンプリングを行うための、HMC法やNUTSなどのアルゴリズムについても書きたいと思います。
