読書メモ: ゼロからできるMCMC

ゼロからできるMCMCの読書メモ。

MCMCの基本条件

MCMCの目的は解析的に計算できない複雑な確率分布の期待値（積分計算）を計算することである。基本的には分布の確率変数を沢山生成してそのときの分布の値を計算し、平均することで期待値を求める。確率変数$x$の確率分布$P(x)$について、MCMCは以下の条件を満たしている必要がある。

• 探索するときの確率変数$x$の過程がマルコフ連鎖であること
• 規約性：あらゆる$x, x'$は有限回のステップで移動できること
• 非周期性：1ステップで任意の$x$に移動できること
• 詳細平衡条件：マルコフ連鎖の遷移確率$T$が$P(x)T(x{\rightarrow}x') = P(x') T(x'{\rightarrow}x)$を満たしていること。遷移に偏りがないこと
  • 厳密には平衡条件: $E[P(x)T(x{\rightarrow}x')] = E[P(x')T(x'{\rightarrow}x)]$のみでもよい
• 非周期性：1ステップで任意の$x$に移動できること

自己相関

• マルコフ過程にしたがう確率変数xをサンプリングする際に前ステップの値との相関（自己相関）がおきうる
• 自己相関がおきると計算される期待値に偏りがでてしまう
• 自己相関を評価する方法：ジャックナイフ法
  • xの系列を等分割したときの個々グループの平均値の標準誤差
  • この標準誤差が変化する＝自己相関がある、と言える

メトロポリス法

• $P(x) = exp(-S(x))/Z$の期待値を計算する
  • $Z$は分配関数
  • $S$は実変数$x$の連続関数。作用または対数尤度と呼ばれる
• ${\Delta}x$をランダムに選んで $x'=x+{\Delta}x$とする
• 詳細平衡条件を満たすために${\Delta}x$と$-{\Delta}x$が同じ確率で生成される必要がある
• 一様乱数 $r{\sim}U(0,1)$を使って$r<\exp(S(x)-S(x'))$であれば$x'$を採用、そうでなければ棄却でそのまま
  「確率$\min(1, \exp(S(x)-S(x')))$で採用する」と同義
• xが多変数の場合は要素を1つずつ更新するか、まとめて更新するか
  • まとめて更新は${\Delta}x$が大きいと棄却される確率が高くなる
  • 変数ごとに${\Delta}x$の大きさを変えることも出来る

Hyblid Monte Carlo(HMC)法

メトロポリス法の派生版。Hamiltonian Monte Carloともいう。 $P(x)=\exp(-S(x))/Z$の$S(x)$が小さい$x$が期待値計算に寄与するが、事前にどの$x$が良いのかはわからない。そのため$x$の候補を逐次的に更新して良い候補に変換する。

• $x$の共役運動量$p$をガウス分布に従ってサンプル
• 初期時刻$\tau=0$でのハミルトニアン$H_{init}=S(x)+\frac{1}{2}\sum_i{p_i^2}$を求める
• リープフロッグ法を使って時刻$\tau$に沿って$x, p$を計算する
• 最終時刻のハミルトニアン$H_{fin}$を計算
• メトロポリステスト：確率$\min(1, \exp(H_{init}-H_{fin}))$で$x$を採用して更新

ハミルトニアン$H$は位置エネルギー＋運動量の形になっている。つまり地形$S$の上で運動量$p$で$x$を押し出し、運動方程式に従ってシミュレーションをして$x, p$を更新していく。

シミュレーションは離散化したハミルトン方程式(運動方程式）に従って$x, p$をリープフロッグ法で計算する。シミュレーション時間$\tau$で変化する$x, p$を$x(\tau), p(\tau)$と書くことにする：

1. $x(\frac{{\Delta}\tau}{2}) = x(0) + p(0)\frac{{\Delta}\tau}{2}$
2. $n = 1, 2, \cdots, N_{\tau}-1$について以下を繰り返す
   
   $$p(n{\Delta}\tau) = p((n-1){\Delta}\tau) - \frac{{\partial}S}{\partial x}((n-\frac{1}{2}){\Delta}\tau){\Delta}\tau$$
   
   
   $$x((n+\frac{1}{2}){\Delta}\tau) = x((n-\frac{1}{2}{\Delta}\tau) + p(n{\Delta}\tau){\Delta}\tau$$
3. 最終時刻は
   
   $$p(N_{\tau}\Delta\tau) = p((N_{\tau}-1)\Delta\tau)-\frac{\partial S}{\partial x}((N_{\tau}-\frac{1}{2})\Delta\tau)\Delta\tau$$
   
   
   $$x((N_{\tau} + \frac{1}{2})\Delta\tau) = x((N_{\tau}-\frac{1}{2}\Delta\tau) + p(N_{\tau}\Delta\tau)\frac{1}{2}\Delta\tau$$

$x, p$を0.5ステップずつ交互に更新していくイメージ。$N_{\tau}, \Delta \tau$はハイパラで、計算量と$x$のステップ幅（幅が大きいと自己相関が少ない）とのトレードオフがある。リープフロッグ法は時間の可逆性がある。つまり時間$\tau$を逆戻しに計算しても同じ結果が得られる。この可逆性がないとMCMCの基本要件である詳細平衡条件が成り立たなくなってしまう。

Gibbs Sampling

多変数の$x={x_1, x_2, \cdots, x_n}$について、$p(x_i|x_{j{\neq}i})$で個々の$x_i$をサンプリングして更新していく。確率分布$p(x_i|x_{j{\neq}i})$が計算できる場合に使える。メトロポリス法と異なり、$x$の候補が棄却されることがない。

Metropolis Hastings(MH)法

メトロポリス法では詳細平衡条件を満たすために$\Delta x$と$-\Delta x$が同じ確率で生成される必要があったが、そうではない場合に拡張する。$x$を更新する提案$x'$を生成する確率$f(x \rightarrow x')>0$が与えられているとき、メトロポリステストを

$$\min(1, \frac{\exp(-S(x'))f(x' \rightarrow x)}{\exp(-S(x))f(x \rightarrow x')})$$

の確率で行う。

$f$をうまく選べば棄却率を下げることができる。Gibbs SamplingやHMC法はこのMH法の特別な場合である。

ベイズ更新への応用

パラメータ$\theta$をもつ事前分布（確率モデル）について観測$\{x_1, \cdots, x_n\}$が与えられた場合の尤もらしいパラメータをMCMCで計算することができる。事後確率$P(\{x_1, \cdots, x_n\}|\theta)$とパラメータの事前確率$P(\theta)$の積を$\theta$の確率分布とみなして上記のようなMCMCを適応すればよい。ただし$\theta$の定義に従って候補を生成する必要があることに注意。たとえばメトロポリス法であれば$\theta$の範囲が間違っていれば棄却する手順を加える。

イジング模型における臨界減速とWolfアルゴリズム

$\{+1, -1\}$のいずれかを状態$s_i$とする格子$i$が並んだイジング模型のエネルギー

$$E=-\sum_{\{i, j\}}s_i s_j-h\sum_i s_i$$

をMCMCでシミュレーションすることを考える。${\{i, j\}}$は隣接する格子インデックスの組み合わせで、$h$は外部エネルギー。メトロポリス法やギブスサンプリングでシミュレーション自体はできるが、$s$全体の振る舞いが大きく変わる相転移のときに自己相関が大きくなってしまう。この現象を臨界減速という。

臨界減速を防げるMCMCの手法としてWolfアルゴリズムがある。このアルゴリズムはランダムに1つの格子を選び、隣接する同じ状態の格子をクラスターに追加、追加した格子に隣接する同じ状態の格子を追加・・・というようにクラスターを広げられなくなるまで大きくする。そしてそのクラスターのすべての格子の状態をまとめて反転する。このように複数の変数をクラスタリングしてまとめて更新する方法をクラスター法という。

Simmulated annealing(焼きなまし)法

メトロポリス法においては作用$S(x)$のとる値が小さい領域に$x$があるとメトロポリステストで棄却されることが多くなり、収束が遅くなるという問題がある。そこで$S(x)$の代わりに温度パラメーター$T$を導入し、$S(x)/T$を作用としてMCMCを行う。そして$T$を段階的に小さくしていくことで収束を効率化する。

レプリカ交換法

焼きなまし法では温度パラメータ$T$のスケジューリングが難しいという問題がある。複数の温度パラメーターを並行して使う:

1. 対象としたい変数$x$と同じ変数の組$x_1, x_2, \cdots, x_M$、これらとペアになる温度パラメータ$T_i$を$T_{m-1}>T_m$となるように用意し、作用$S(x_1, x_2, \cdots, x_m)=\sum_{m=1}^M\frac{f(x_m)}{T_m}$を定義する。
2. 上記の変数$x_1, x_2, \cdots, x_m$と作用$S$を使ってMCMCによる更新を行う
3. 確率$\min(1,\exp(-\Delta{S}))$で$x_m$と$x_{m+1}$を入れ替える。ただし
   
   $$\Delta{S}=(\frac{1}{T_{m+1}}-\frac{1}{T_m})(f(x_m)-f(x_{m+1}))$$

まとめ

MCMCは解析的に計算できない複雑な確率分布の積分を変数をマルコフ連鎖で生成することで計算する方法。変数をランダムに生成してしまうと分布の形次第で高い棄却率や自己相関で効率が悪い。そのため変数の系列（この本では配位と呼ぶ）をMCMCの基本条件を壊さずいかに生成するかで手法が分かれる。