MCMCサンプリングのMetropolis–Hastings法
どうも皆さんこんにちは、アトラエでデータサイエンティストをしている芝原です!
今日はベイズ統計のMCMCサンプリングで出てくるMetropolis–Hastings法について書きたいと思います。参考文献はこちらです。
Understanding the Metropolis-Hastings Algorithm, The American Statistician, Vol. 49, No. 4 (Nov., 1995), pp. 327-335.
皆さんが使うときはコンピュータが全てやってくれるので、特に意識することもないかもしれませんが、先人の努力やアイデアが垣間見えて面白かったりします。そこからどこが問題なのか考える際に非常に有益なヒントを与えてくれたりもします。というわけで今回は、これどうやって動いているんだろう??生で動かしてみよう!!あれ!?うまくいかない!?みたいなことをやっていきたいと思います。
今回の記事ではある程度ベイズ統計やマルコフ連鎖を授業や独学などで少し触ったことがあるという人を対象にしています。なので、ベイズの定理などは軽く触れますが、大雑把な説明です。証明もしません。すいません。知っていることの復習だったり頭の整理に使っていただければ嬉しいです。
こちらはAtrae Advent calendar 2024の一環として書かれたものになります!是非興味が湧いた時に覗いていってください!
ベイズ統計の復習
今、$${\theta}$$を推定したいパラメーター、$${y}$$をデータとします。尤度関数、事前分布、事後分布をそれぞれ$${L(y \mid \theta), p^{\text{prior}}, p^{\text{post}}}$$と表すことにします。
ベイズの定理より
$$
p^{\text{post}}(\theta) \propto L(y \mid \theta) \times p^{\text{prior}}(\theta)
$$
となります。ここで関数$${f, g}$$に対して、$${f{\propto}g}$$は定数倍の誤差を許す記号で、数学的にはある$${C >0}$$が存在して$${f=Cg}$$となる、という意味です。この$${C }$$の役割は$${L(y \mid \theta) \times p^{\text{prior}}(\theta)}$$を正規化して確率密度関数にすることに他なりません。
積分が無理問題
ベイズ統計に慣れてくると、「事後分布を手計算で出してみよう!」や「Posterior Meanを計算しよう!」などなど色々計算したくなってきます。その際、確率分布として正規化している定数$${C}$$も必要になってきます。共役事前分布と呼ばれる都合のいいもの以外で立ちはだかる壁が「この$${C}$$が計算できない😇」というものです。
この$${C}$$が計算できない問題をマルコフ連鎖の発想を胸に、アルゴリズム的に(近似的に)解決してくれるのがMCMC (Markov Chain Monte Carlo)です。MCMCのMetropolis-Hastings法というものを用いて、このNoteの最終章で頑張って幾つかの分布を近似するので、楽しみにしておいてください!
マルコフ連鎖の復習
マルコフ連鎖とは
例えば、あなたが引越しを毎年するとします。今の居住地から色々な要素を踏まえて、次の引越し先を決定します。この選択には確率的な揺らぎがあるので、それを表現したのが大雑把にマルコフ連鎖です。
上で、Town $${i}$$からTown $${j}$$にいく確率を$${P_{ij}}$$と表しています。この$${P=({P_{ij}})_{i,j}}$$を確率遷移行列と言います。
つまり、(離散)マルコフ連鎖とは時刻$${t+1}$$での位置が時刻$${t}$$の位置のみによって決定され、それ以前の挙動とは無関係なもののことを指します。数学的には$${\{X_n\}_{n \in \mathbb{N}}}$$という状態空間$${S}$$に値をとる確率変数の列があった時、$${\{X_n\}_{n \in \mathbb{N}}}$$がマルコフ連鎖であるとは、全ての$${n \in \mathbb{N}}$$と$${x_1, \cdots, x_n, x \in S}$$に対し
$$
\text{Pr}(X_{n+1}=x \mid X_n=x_n, \cdots, X_1=x_1)=\text{Pr}(X_{n+1}=x \mid X_n=x_n)
$$
が成り立つことと定義します。ここで$${\text{Pr}(A\mid B)}$$とは$${A}$$の$${B}$$による条件付き確率のことです。
マルコフ連鎖の定常分布とDetailed Balance Equation
(状態空間$${S}$$が有限で既約な)マルコフ連鎖の$${n \to \infty}$$の挙動として重要なのが「定常分布$${\pi}$$を持つ」というものです。また町の例で説明してみます。あなたがTown 2にいたとして、確率遷移行列により、Town 3に移動したとします。この時あなたの分身をTown 2に残しておくとして、これを繰り返すと以下のようになります。
これを$${n}$$回繰り返した後、それぞれのTownにいる自分の分身の数を$${n}$$で割ってあげたものがおおよその定常分布$${\pi}$$になります。なぜならば、これ、もう1回同じ操作をしてこの分身たちを1回ずつずらしても、ほとんど分布は変わらないですよね。これを数式で書くと、
$$
\pi =\pi P \ \ \text{つまり任意の$j$に対して} \ \pi_{j}= \sum_{i} \pi_i P_{ij}.
$$
となります。これは「今Town $${j}$$にいる人の比$${\pi_j}$$は次の時に色んなTown $${i}$$からTown $${j}$$に流入してくる人の比の和$${\sum_{i} \pi_i P_{ij}}$$と同じ」ということになります。定常って感じですよね。これはBalance Equationとも呼ばれます。数学的には(状態空間$${S}$$が有限の場合は)定常分布$${\pi}$$は確率遷移行列の固有値1の固有ベクトルに他なりません。
これよりも強い条件として、Detailed Balance Equation
$$
\pi_j P_{ji}=\pi_i P_{ij}
$$
というものがあります。つまりTown $${j}$$からTown $${i}$$への流入がTown $${i}$$からTown $${j}$$への流入と釣り合っているということです。上で$${i}$$に関して和を取るとBalance Equationを満たすことがすぐに分かります。このDetailed Balance Equationを使ってMetropolis-Hastings法のアルゴリズムを解説していきます!
発想:マルコフ連鎖の定常分布として事後分布を見つける
ここでベイズ統計での問題を思い出します。事後分布$${p^{\text{post}}}$$を求める際に積分ができないというところが問題でした。なので、この積分計算をうまく避けたアイデアが必要となります。
クールな発想は「定常分布$${\pi}$$として事後分布$${p^{\text{post}}}$$を持つような、マルコフ連鎖を発生させれば良い!」ということです。でも求めたい分布が計算できなくて困っているのに、どうやってそれに対応するマルコフ連鎖を発生させるのでしょうか?少し現状を整理すると
私たちは事後分布$${p^{\text{post}}}$$をexactには計算できない。
でも、定数倍を覗いて$${p^{\text{post}}}$$を知っている(つまり$${L(y\mid \theta) \times p^{\text{pior}}}$$)。
といった状況です。なので、$${L(y\mid \theta) \times p^{\text{pior}}}$$のみを用いて定常分布$${\pi}$$として事後分布$${p^{\text{post}}}$$を持つマルコフ連鎖を発生させれば良い!
といった感じになります。ではこれはどうするのでしょうか?
Metropolis–Hastings法
確率遷移行列の出し方
今やりたいことは求めたい分布$${\pi=p^{\text{post}}}$$を実現するように確率遷移行列$${P=(P_{ij})_{ij}}$$を定めたい、ということです。
Detailed Balance Equationを軸に確率遷移行列$${P=(P_{ij})_{ij}}$$を発生させるのが、Metropolis–Hastings法というものです。詳しく説明していきます。まずDetailed Balance Equationから
$$
\pi_i P_{ij}=\pi_j P_{ji} \ \ \text{変形すると} \ \ \frac{\pi_i}{\pi_j}=\frac{P_{ji}}{P_{ij}}
$$
ここで、確率遷移行列の各成分$${P_{ij}}$$が「提案分布$${Q(j \mid i)}$$」と「受け入れ分布$${A(i, j)}$$」の積で書かれていると仮定します。つまり
$$
P_{ij}=Q(j \mid i) \times A(i, j)
$$
これの意味は割りと簡単で、($${i}$$から$${j}$$に移る確率$${P_{ij}}$$)=(神から$${i}$$から$${j}$$に移るよう提案される確率$${Q(j \mid i)}$$)$${\times}$$(それを自分が受け入れる確率$${A(i, j)}$$)と思うことができます。なので、確率的推移$${P_{ij}}$$が$${Q(j \mid i)}$$と$${A(i, j)}$$の2段階で決まっていく感じですね。今回は単純のために、$${Q(j \mid i)=Q(i \mid j)}$$とします。これとDetailed Balance Equationを合わせて
$$
\frac{\pi_i}{\pi_j}=\frac{P_{ji}}{P_{ij}}=\frac{Q(i \mid j) \times A(j, i)}{Q(j \mid i) \times A(i, j)}=\frac{A(j, i)}{A(i, j)}
$$
となります。なので、$${A(i, j)}$$の1つとして、
$$
A(i, j)=\min\Big\{\frac{\pi_j}{\pi_i}, \ 1\Big\}
$$
と取れば良いことが分かります。ここで1番重要なのが、「定数倍を覗いて$${\pi(=p^{\text{post}})}$$を知っている(つまり$${L(y\mid \theta) \times p^{\text{pior}}}$$)」ということなのです。つまり、$${A(i, j)}$$は現時点で$${L(y\mid \theta) \times p^{\text{pior}}}$$を用いて計算可能なのです!
Metropolis-Hastings法では、提案分布$${Q(j \mid i)}$$は自分で決めるので、これで求めたい確率遷移行列$${(P_{ij})_{ij}}$$が手に入ったことがわかります。
Metropolis–Hastings法の概要
上で確率遷移行列$${(P_{ij})_{ij}}$$を求めることができました。あとはその遷移に沿うようにサンプルしていけばよろしいです。ここまで来れば、アルゴリズムがどう動いているかを知るのは簡単です。
$${x_0=i}$$と適当に取る。
$${Q(x \mid x_0)}$$からランダムに$${x=j}$$をサンプルする。
$${\min\Big\{\frac{\pi_j}{\pi_i}, \ 1\Big\}}$$を計算する。
区間[0, 1]からランダムな数$${u}$$を発生させて、$${u}$$が$${\min\Big\{\frac{\pi_j}{\pi_i}, \ 1\Big\}}$$以下なら受け入れ($${x_1=j}$$に移り)、そうでなければ却下($${x_1=i}$$に留まり、サンプルとしてカウントしない)。
これを繰り返す。
数学的に深いマルコフ連鎖の一連の結果から始まりベイズ統計でこんなシンプルなアルゴリズムが活躍するって素敵ですよね!
実際試してみよう!
Metropolis-Hastings法がうまくいく例
ここでは簡単な例でMetropolis-Hastings法がうまくいっていることを確認します。今回実現したい分布$${\pi}$$を平均0、分散1の正規分布$${\mathcal{N}(0, 1)}$$として、とりあえず100000回Metropolis-Hastings法をぶん回してみます。提案分布$${Q(x' \mid x)}$$として$${\mathcal{N}(x, 1)}$$と設定しています。
def target_distribution(x):
return np.exp(-x**2 / 2) / np.sqrt(2 * np.pi)
def metropolis_hastings(initial, n_samples, proposal_width):
samples = []
current = initial
for _ in range(n_samples):
u=np.random.rand()
proposal = current + np.random.normal(0, proposal_width)
acceptance_ratio = target_distribution(proposal) / target_distribution(current)
if u < acceptance_ratio:
current = proposal
samples.append(current)
return np.array(samples)
n_samples = 100000
initial = 0
proposal_width = 1.0
samples = metropolis_hastings(initial, n_samples, proposal_width)
取り出したサンプルのヒストグラムが下の図です。いや〜。渋いですね。
Metropolis-Hastings法が(時々)うまくいかない例
アルゴリズムの理解は何かがうまくいかなかった時に真価を発揮します。ということで、今回は$${\mathcal{N}(-5, 1)}$$と$${\mathcal{N}(5, 1)}$$を合わせたBimodalな確率分布を考えます。今回は提案分布$${Q(x' \mid x)}$$として$${\mathcal{N}(x, 0.5)}$$と設定しています。
def target_distribution(x):
return return 0.5 * (np.exp(-(x+5)**2 / 2) / np.sqrt(2 * np.pi) + np.exp(-(x-5)**2 / 2) / np.sqrt(2 * np.pi))
def metropolis_hastings(initial, n_samples, proposal_width):
samples = []
current = initial
for _ in range(n_samples):
proposal = current + np.random.normal(0, proposal_width)
acceptance_ratio = target_distribution(proposal) / target_distribution(current)
if np.random.rand() < acceptance_ratio:
current = proposal
samples.append(current)
return np.array(samples)
n_samples = 100000
initial = 0
proposal_width = 0.5
samples = metropolis_hastings(initial, n_samples, proposal_width)
Metropolis-Hastings法で生成されたヒストグラムがこちらになります。激しく失敗しています。
なぜ左側から一切サンプルできなかったのでしょうか?これは提案分布$${Q(x' \mid x)}$$の取り方にあります。アルゴリズムではまず、提案分布$${Q(x' \mid x)}$$から次のサンプルの候補を出します。今、$${Q(x' \mid x)}$$は$${\mathcal{N}(x, 0.5)}$$なので、右側の山が現在地の場合、左側の山のどこかが提案されることは(ほぼ)ありません。それに加えて、提案された位置を受け入れるかどうかに関してはアルゴリズムにより、分布の山の高い方に動きがちなので、右側の部分からしかサンプルされなかった、ということになります。なるほどー。
解決策はシンプルで、提案分布を設定する段階で、もう少し遠くの位置も提案できるようにしてやればいいのです。なので下では$${Q(x' \mid x)}$$として、$${\mathcal{N}(x, 3)}$$として、Metropolis-Hastings法を動かしてみます(上のコードでproposal_width=3に変える)。すると、以下のような感じでした。セクシー!
終わりに
いかがでしたか?今回はベイズ統計のMCMCの必要性からMetropolis-Hastings法を実際に動かしてみることをやってみました。直感や理解の助けになれば嬉しいです!
採用情報及び会社情報はこちらになります。興味が湧いた方は是非コンタクトしてみてください!