🔖

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

2022/05/23に公開約4,500字

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

MCMCの基本条件

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

  • 探索するときの確率変数xの過程がマルコフ連鎖であること
  • 規約性:あらゆるx, x'は有限回のステップで移動できること
  • 非周期性:1ステップで任意のxに移動できること
  • 詳細平衡条件:マルコフ連鎖の遷移確率TP(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))/ZS(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の上で運動量pxを押し出し、運動方程式に従ってシミュレーションをしてx, pを更新していく。

シミュレーションは離散化したハミルトン方程式(運動方程式)に従ってx, pをリープフロッグ法で計算する。シミュレーション時間\tauで変化するx, px(\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_iT_{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_mx_{m+1}を入れ替える。ただし
    \Delta{S}=(\frac{1}{T_{m+1}}-\frac{1}{T_m})(f(x_m)-f(x_{m+1}))

まとめ

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

Discussion

ログインするとコメントできます