マルコフ連鎖
モンテカルロ法
乱数を用いた数値計算手法について
Created by Kazuhito Takeuchi
このスライドの目的
- （学部）モンテカルロ法に触れてもらう
- （修士）MCMCについてなんとなく理解してもらう
目次
- モンテカルロ法とは？
- モンテカルロ法の限界
- マルコフ連鎖モンテカルロ法
- 不変分布の一意性
- まとめ
モンテカルロ法とは？
Monte Carlo (MC)
- 乱数を使った計算手法の総称
- 数値計算（シミュレーション）では擬似乱数を使う
- 線形合同法（一番単純なやつ・低コスト）
- メルセンヌ・ツイスタ法（すごいやつ・若干高コスト）
- xorshift（そこそこすごくて低コスト）
- など…
- 要はサイコロを転がして色々計算しようという手法
円周率の計算
$\pi = 0$
$ x, y \in [0, 1] $
while $ N $ :
if $x^2 + y^2 \leq 1 $ :
$\pi += 1 $
試行回数100回：$ \pi = 2.96 $
円周率の計算
$\pi = 0$
$ x, y \in [0, 1] $
while $ N $ :
if $x^2 + y^2 \leq 1 $ :
$\pi += 1 $
試行回数1000回：$ \pi = 3.168 $
円周率の計算
$\pi = 0$
$ x, y \in [0, 1] $
while $ N $ :
if $x^2 + y^2 \leq 1 $ :
$\pi += 1 $
試行回数10000回：$ \pi = 3.1356 $
円周率の計算
試行回数を増やすと誤差が減る
単純なモンテカルロ法(MC)の限界
半径1の$N$次元球の体積$$ V_N = \frac{\pi^{N/2}}{(N/2)!} $$
球を囲む$N$次元立方体の体積は$2^N$
球にヒットする確率$$ \frac{\pi^{N/2}}{2^N(N/2)!} $$
高次元・多変量になるほど計算効率が著しく低下する
例えば…2次元円と3次元球を比較すると分かりやすいかも
合金（スピン系）でのMC
確率高
確率低（ほぼゼロ）
高次元球と同じような問題が合金でも起こる
合金：$ \vec{\sigma} = \{ \sigma_i, \sigma_j, \sigma_k,... \} \quad ( \sigma_i = {\rm A \ or \ B} ) $
確率：$ P(\vec{\sigma}, \beta) \propto \exp (-\beta E(\vec{\sigma})) $
エネルギーの期待値：$ E(\beta) = \sum_{\vec{\sigma}} E(\vec{\sigma}) P(\vec{\sigma}, \beta) $
どこが問題なの？
- エネルギーの期待値を正しく評価したい
- 重みが効くところをサンプルしたい
- 例) 秩序的な配列を持つとエネルギーが下がる合金の場合
- 適当に金属を並べた時，規則正しく配列する確率はどれくらいか？
- ２元系だと $\propto \frac{1}{2^N}$でとても低い
- 円周率の次元/原子の数の多さ = 多自由度系における困難
- つまり，適当にサンプルすると無駄が多い...
- いい感じにサンプルするアルゴリズムを探そう！
- 期待値を計算する手法としての「マルコフ連鎖モンテカルロ法」
マルコフ連鎖モンテカルロ法(MCMC)
マルコフ連鎖：現在の状態は１つ前の状態のみによる
- 【例】今日は中華だったから,明日は和食にする確率は$1/3$
- 今日は中華だったから,明日も中華にする確率は$2/3$
- 今日は和食だったから,明日も和食にする確率は$1/4$
- 今日は和食だったから,明日は中華にする確率は$3/4$
- おとといの食事は考えない
- $A_t$と$B_t$はそれぞれ和食・中華にする確率
マルコフ連鎖モンテカルロ法(MCMC)
この漸化式を解くと不変分布$A_{\infty}, B_{\infty}$が得られる
$$ \begin{pmatrix} A_{\infty} \\ B_{\infty} \\ \end{pmatrix} = \begin{pmatrix} \frac{3}{11} \\ \frac{8}{11} \\ \end{pmatrix} \quad \begin{pmatrix} A_{\infty} \\ B_{\infty} \\ \end{pmatrix} = \begin{pmatrix} \frac{1}{3} & \frac{1}{4} \\ \frac{2}{3} & \frac{3}{4} \\ \end{pmatrix} \begin{pmatrix} A_{\infty} \\ B_{\infty} \\ \end{pmatrix} $$
不変分布が分かれば期待値が計算できる
不変分布は固有ベクトルに対応
マルコフ連鎖モンテカルロ法(MCMC)
今度は漸化式を解かずに，実際にシミュレーションすると
サンプル数の割合が不変分布に従っていることがわかる
| ステップ数 | 1 | 2 | ... | t | t+1 | ... |
|---|---|---|---|---|---|---|
| 状態 | A | B | ... | B | A | ... |
MCMCまとめ
- 一つ前の状態にのみ依存するモデル
- (期待値に効く近くの状態も効くだろうという仮定)
- 不変分布（確率分布）を明らかにして期待値を計算
- 完全な状態行列があれば漸化式（固有値問題）を解く
- なければサンプルした時系列から期待値を計算
- （これが非常に強力 サイコロを転がすだけで良い）
合金へMCMCを応用しよう
確率（再掲）：$ P(\vec{\sigma}, \beta) \propto \exp (-\beta E(\vec{\sigma})) $
これをカノニカル分布と呼ぶ
目的：MCMCを使って期待値（内部エネルギー）を計算する
- さっきと同じように漸化式を解けば良い？
- 状態の数が$2^N$あるので大きい系では現実的に無理…
- でも，サンプルの時系列が不変分布に対応している
- 不変分布がカノニカル分布になるように遷移確率を設計できれば，サンプルの時系列から期待値を計算できる
合金へMCMCを応用しよう
遷移確率を以下のように設定すると，不変分布はカノニカル分布になる（メトロポリス法）
$$ \pi (\vec{\sigma} \to \vec{\sigma}') = {\rm min} \{ 1, \exp( - \beta( E(\vec{\sigma}') - E(\vec{\sigma}) ) ) \} $$後述する以下の詳細釣り合い条件を満たすように設計されている
$$ P(\vec{\sigma}) \pi (\vec{\sigma} \to \vec{\sigma}') = P(\vec{\sigma}') \pi (\vec{\sigma}' \to \vec{\sigma}) $$今の原子配置から適当な異種原子同士を交換して１ステップとする
| ステップ数 | 1 | 2 | ... | t | t+1 | ... |
|---|---|---|---|---|---|---|
| 合金 | $\vec{\sigma}_1$ | $\vec{\sigma}_2$ | ... | $\vec{\sigma}_t$ | $\vec{\sigma}_{t+1}$ | ... |
| エネルギー | $E(\vec{\sigma}_1)$ | $E(\vec{\sigma}_2)$ | ... | $E(\vec{\sigma}_t)$ | $E(\vec{\sigma}_{t+1})$ | ... |
期待値 ： $\langle E \rangle = \sum_{t=1}^N E(\vec{\sigma}_t) / N$
CuAu合金
@ T=600[K]
期待値 ： $ E(\beta) = \sum_{\vec{\sigma}} E(\vec{\sigma}) P(\vec{\sigma}, \beta) $
CuAu合金
高温状態と低温状態
各温度で期待値をプロット
一次の規則不規則相転移をMCMCで表現できた
MCMCは常に使えるの？
- いいえ
- とはいえかなり応用範囲が広い
- 不変分布が一意に存在するためには以下の条件が必要
- エルゴード性（全ての状態へ遷移できるか）
- 釣り合いの条件（確率の収支があっているか）
- スピンフリップ/原子の交換ではエルゴード性は満たされる
- より厳しい詳細釣り合い条件がよく使われる（設計しやすい）
不変分布の一意性
証明のポイントは初期分布$P_0(x)$から$m$回遷移したあとの$P_m(x)$について，不変分布以外の部分がゼロに収束するように表現する点
$$ P_1(x) = cP(x) + (1-c)R_1(x) $$ $$ P_2(x) = \{ 1- (1-c)^2 \} P(x) + (1-c)^2 R_2(x) $$......
ここで
- $P(x)$は不変分布
- $ 0 \lt c \lt 1 $
- $ R_1(x) $ は$x$の確率分布
- $ m \to \infty$ で$P_m(x)$は不変分布に収束
それっぽい証明の手順
すごく大雑把なので細かい話は調べて下さい
- 分布をベクトル，分布の遷移を行列として扱える
- （状態と遷移行列はグラフ理論の点と隣接行列に対応）
- $P$を不変分布（ベクトル），$T$を遷移行列として$P = TP$
- 固有値を$\{ \lambda_\alpha \}$，固有ベクトルを$\{ V_\alpha \}$とする
- 任意の初期分布を正規直交系の固有ベクトルで展開
- ここで不変分布$P$の固有値は1であることに注意して
- $P_0 = P + \sum_{\alpha} c_\alpha V_\alpha$
- $P_m = T^m P_0 = P + \sum_{\alpha} c_\alpha \lambda_\alpha^m V_\alpha$
- ペロン・フロベニウスの定理より$|\lambda_\alpha|<1$
MCMCの利点と欠点
- 応用範囲が広い/実装が簡単
- 期待値が計算できる
- 分布の局所性を仮定しているので，多峰性の分布のサンプルに向かない
- 解決策：flat-histogram method (e.g., multi-canonical, Wang-Landau...)
- モーメントしか計算できない（期待値，分散など）
- グラフ上のノードにボトルネックが存在すると収束が遅くなる
- i.e., $\{ \lambda_\alpha \}$の一番高い固有値の大きさによって収束速度が変わる
- 解決策：詳細釣り合いを満たさないMCMC，Event-chain MC...
まとめ
- モンテカルロ法は乱数を使った数値計算手法
- マルコフ連鎖モンテカルロ法(MCMC)は多自由度の計算が得意
- MCMC法は不変分布からのサンプリング手法
- エルゴード性
- （詳細）釣り合いの条件
- 遷移確率を設計して，不変分布をカノニカル分布する
MCMCを使う時は，どの不変分布からサンプルしているかを常に意識しよう！
参考文献
- 伊庭幸人ほか(2005)『計算統計 2 マルコフ連鎖モンテカルロ法とその周辺 (統計科学のフロンティア 12)』岩波書店．
- teramonagi「マルコフ連鎖モンテカルロ法入門-１」，