MCMC の基本

MCMC の目標とするところは与えられた確率分布 $p(z)$ からその分布に従ったデータの列 $z ^ {(i)}$ をサンプリングすることです．MCMCのアルゴリズムを上手に設計すれば高次元の分布からも効率的にサンプリング可能なことが知られています．

例えばマラソンでは観測結果 $X$ とそれを生成したパラメータ $\theta$ があった場合，その事後分布 $p(\theta | X) \propto p(\theta) p(X | \theta)$ から $\theta$ を MCMC でサンプリングすることで $\theta$ の分布を推定し，その分布下で最適な行動を出力するといったことに使われます．

MCMC でサンプリングされる列は 1 次マルコフ連鎖に従って生成されます．つまり $m+1$ 番目のサンプル $z ^ {(m+1)}$ は1 つ前のサンプルである $z ^ {(m)}$ のみに依存して順番にサンプリングされます．数式で書くと次のようになります．
$$p(z ^ {(m+1)} | z ^ {(1)}, z ^ {(2)}, \dots , z ^ {(m)}) = p(z ^ {(m+1)} | z ^ {(m)})$$
このとき遷移確率を $T(z^{(m)} , z ^ {(m+1)}) = p(z ^ {(m+1)} | z ^ {(m)})$ と置き，確率分布 $p ^ * (z)$ が各ステップで一定であるならば $p ^ * (z)$ は不変と呼ばれます．言い換えると $p ^ * (z) = \sum _ {z'} T(z', z) p ^ * (z')$ です．

確率分布 $p ^ * (z)$ が遷移確率 $T$ の下で不変になる十分条件に以下の詳細釣り合い条件があります．¹
$$p ^ * (z) T(z, z') = p ^ * (z') T(z', z)$$
必要な確率分布からサンプリングしたい場合はその分布がこの詳細釣り合い条件を満たすように遷移確率を設計すればよいことがわかります．²

メトロポリス・ヘイスティング法

メトロポリス・ヘイスティング法は最も基本的な MCMC アルゴリズムです．³

サンプリングしたい確率分布の定数倍 $\tilde{p} (z)$ がわかっており，それを正規化した分布である $p(z) = \frac{\tilde{p} (z)} {\int \tilde{p} (z) dz}$ からサンプリングをしたいとします．このアルゴリズムのマルコフ連鎖の過程では現在の状態が $z ^ {(m)}$ から次の状態の候補 $z ^ * $ を提案分布 $q(z ^ * | z ^ {(m)})$ から選び，以下の確率 $A(z ^ * , z ^ {(m)})$ で $z ^ * $ を受理して $z ^ {(m+1)} = z ^ * $ とし，受理しないときは $z ^ * $ を棄却して $z ^ {(m+1)} = z ^ {(m)}$ とします．
$$A(z ^ * , z ^ {(m)}) = \min \left(1, \frac{\tilde{p} (z ^ * ) q(z ^ {(m)} | z ^ * )} {\tilde{p} (z ^ {(m)} ) q(z ^ * | z ^ {(m)} )} \right)$$

これが詳細釣り合い条件を満たすことは簡単に確かめられます．

\begin{align}
p (z) q(z'| z) A(z', z)  &= \min \left(p (z) q(z'| z), p(z') q(z| z') \right)\\
&= p (z') q(z| z') A(z, z')\\
\end{align}

特に提案分布が対称，つまり $q(z'| z) = q(z| z')$ である場合はメトロポリス法と呼ばれ，$A(z ^ * , z ^ {(m)})$ は次のようになります．
$$A(z ^ * , z ^ {(m)}) = \min \left(1, \frac{\tilde{p} (z ^ * ) } {\tilde{p} (z ^ {(m)} ) } \right)$$

例として，$\tilde{p} (z) = \frac{1}{1 + z ^ 2}$ からサンプリングすることを考えます．⁴ 提案分布 $q(z ^ * | z ^ {(m)})$ を次のように定義するとこれは明らかに対称性を満たすので，メトロポリス法が使えます．

q(z ^ * | z ^ {(m)}) = 
\begin{cases}
1 & (|z ^ * - z ^ {(m)}| < 0.5) \\
0 & else
\end{cases}

Pythonで実装すると以下のようになります．

import random

def p_tilde(z):
    return 1 / (1 + z ** 2)

def metropolis(z0):
    z = z0

    while True:
        z_star = z + random.uniform(-0.5, 0.5)
        A = min(1, p_tilde(z_star) / p_tilde(z))

        # 確率 A で受理
        if random.uniform(0, 1) < A:
            z = z_star
        
        yield z

この結果をプロットすると次のようになり，実際にサンプリングできていることがわかります．

MCMC から焼きなまし法へ

ここからは MCMC をどのように最適化アルゴリズムとして使うのかを考えていきます．
もし MCMC の高次元の分布の探索能力を活用できればよい最適化ができると期待できます．

この章では関数 $f(x)$ を最小化する問題を考えます．このとき $f(x)$ が小さい $x$ ほど多くの回数探索をすれば効率が良くなるはずです．そこで以下の分布からサンプリングを行うこととします．
$$ \tilde{p}(x) = e ^ {-f(x)}$$

提案分布は対称であると仮定しメトロポリス法を使うと，受理確率 $A(z ^ * , z ^ {(m)})$ は以下のようになります．⁵
$$A(x ^ * , x ^ {(m)}) = \min \left(1, \frac{\tilde{p} (x ^ * ) } {\tilde{p} (x ^ {(m)} ) } \right) = \min\left( 1, e ^ {f\left(x ^ {(m)}\right) - f(x ^ * )}\right)$$

しかし実際にこの方法を試してみると最適値の周辺以外に多くの時間を費やしてしまったり，複数の山の間の何もない区間で遷移の棄却率が高まりうまく全体を探索できない問題が発生します．⁶

単純なアルゴリズムの改良

そこで，新たにパラメータとして温度 $T$ を導入し，$\tilde{p}(x)$ を以下のように書き換えます．これはボルツマン分布と呼ばれるものです．
$$ \tilde{p}(x) = e ^ {-f(x) / T}$$

受理確率 $A(z ^ * , z ^ {(m)})$ は以下のようになります．
$$A(x ^ * , x ^ {(m)}) = \min \left(1, \frac{\tilde{p} (x ^ * ) } {\tilde{p} (x ^ {(m)} ) } \right) = \min\left( 1, e ^ {\left(f\left(x ^ {(m)}\right) - f(x ^ * )\right) / T}\right)$$

$T$ を調整すれば $\tilde{p}(x)$ を探索にとって都合の良い形に変形することができます．下の例では裾野が広すぎた分布は温度を下げることで最適値周りの分布を集中させ，また複数の山がある分布は温度を上げることで山の間の溝がやわらげられています．

また，探索の初期では高い温度を使うことで広い範囲を探索し，探索の終盤では低い温度を使うことで最適値の周辺を集中的に探索することにすればよい解が得られそうです．

実はこれは焼きなまし法そのものになります．以下は上にあげたコードを少し変更してこのアルゴリズムをそのまま実装したものですが，よく見る焼きなまし法のコードになっていると思います．

import random
import math


def metropolis(z0):
    z = z0

    for i in range(total_iteration)[::-1]:
        z_star = z + random.uniform(-0.5, 0.5)

        t = start_temperature * (i + 1) / total_iteration
        A = min(1, math.exp((f(z) - f(z_star)) / t))

        # 確率 A で受理
        if random.uniform(0, 1) < A:
            z = z_star
        
    return z

これらから，実のところ焼きなまし法というものはボルツマン分布からメトロポリス法でサンプリングするアルゴリズムだと解釈できるということがわかりました．やったー！

最後に焼きなまし法の状態の推移を可視化した動画を載せます．感覚的にも焼きなまし法が $\tilde{p}(x) = e ^ {-f(x) / T}$ からサンプリングをしているということが理解できると思います．

参考文献

PRML 第 11 章サンプリング法

必要条件ではないです． ↩
このときサンプリングした分布が求めたい分布に収束するかについての議論も必要ですが，ここではできていません．気になる人はエルゴード性で調べてください． ↩
ギブスサンプリングやハミルトニアンモンテカルロ法もメトロポリス・ヘイスティング法の特殊な場合です． ↩
$\tilde{p} (z) = \frac{1}{1 + z ^ 2}$ の場合は $z = \tan\left(\frac{\pi}{2} \rm{uniform}(-1, 1)\right)$ としてもサンプリングできます． ↩
実用の際は遷移確率が非対称なこともありますが，出力の分布が多少ゆがむだけなので焼きなまし法だと問題にはなりにくいです．ただし MCMC でサンプリングをする際は気を付ける必要があります． ↩
複数の山がある分布(多峰性分布)でうまくサンプリングできない問題は MCMC でも発生しうる問題です．レプリカ交換モンテカルロ法ではいったん高い温度を経由させることでこの問題を解決しています． ↩