MCMC入門

記事の内容


理論から技術的な部分まで広く解説した記事を書いてみました. 間違いがあれば教えていただきたいです(Twitterで). サンプルコードは過去に書いたものを引っ張ってきました. 書き方がバラバラだったりします. 不具合・間違い等あれば, 時々更新するつもりです. 詳しく, かつ正しい内容は参考文献をご覧ください.

MCMCをどこで使うか

動機付け

MCMC(Markov Chain Monte Carlo; Markov連鎖モンテカルロ法)は, ある分布からのサンプル(分布の実現値たち)を計算する手法の総称です. 総称なので, 何種類もあります. 代表的なものは「具体的な手法」の節で取りあげます. 得られたサンプルはMarkov連鎖になっています. つまり, 直前に計算したサンプルにのみ依存します. ある分布というのは, データを解析する側が, サンプルを得ようとしている分布です. 以下, 「目的の分布」と呼ぶことにします.

以下は, MCMCを使いたい場面の例です.

  • Bayes推論を行っていて, 事後分布の平均を計算したい. しかし, 事後分布は式が複雑で, 手計算では対処できない.
  • 予測分布の対数の期待値を計算する必要が出てきた. しかし, 予測分布は式が複雑で, 手計算では対処できない.
  • 事後分布の標準偏差を計算することで, 推定の不確かさを評価したい. しかし, 事後分布は式が複雑で, 手計算では対処できない.

上の場面の共通点は, 手計算では対処できないくらい複雑な, 積分計算が必要という点です. こんな時に役立つのがMCMCです. MCMCは, 分布からの実現値を使って, 平均や期待値を近似計算します. うまくいけばそこそこの精度で計算できます. この記事では, MCMCの基本的な考え方, 具体的な計算方法等をざっくり解説します.

概要

MCMCの概要を説明します. 例えば, データを解析する人が, 次のような期待値を計算するとします. \(p\)は目的の分布の密度関数とします.

\begin{equation} \mathbb{E}\left[f(X) \right] = \int f(x)p(x)dx \end{equation}

もし, 目的の分布からのサンプル\(\{x^{(s)}\}_{s=1}^{S}\)が手に入れば, 次のような近似計算が可能です.

\begin{equation} \mathbb{E}\left[f(X) \right] = \int f(x)p(x)dx \simeq \frac{1}{S} \sum_{s=1}^{S} f(x^{(s)}) \end{equation}

近似の精度をよくするためには, 分布内を満遍なく探索してサンプルが欲しいところです. しかし, 計算時間や記憶量の観点から, なるべく効率的に計算できた方が良いです. そこで, MCMCに期待する挙動は以下の通りです.

  • 積分の値に大きく寄与する部分, つまり確率の大きい部分まで移動する.
  • 確率の大きい部分まで移動したら, その周辺をうろちょろする.

移動途中は初期値の影響を受けており, 目的の分布からのサンプルではありません. 段々と目的の分布に近づきます. 移動が完了して, ある程度サンプルが集まれば, 目的の分布に落ち着いていると考えられます. 分布が落ち着く時, この分布を定常分布と言います. また, 確率の大きい部分をうろちょろするという期待に応えるために, Markov連鎖を用います. Markov連鎖は直前の状態にのみ依存して移動します.

以上より, MCMCの基本方針は次のように要約できます. 「目的の分布を定常分布とするMarkov連鎖を構成する」ことです.

MCMCの基礎理論

Markov連鎖

まず, MCMCのアルゴリズムがどのように設計されていつかを概観します. 前節で触れたように, MCMCの基本方針は目的の分布を定常分布とするようなMarkov連鎖を構成することです.まずは, このMarkov連鎖について簡単に定義します. 曖昧で不正確な定義なので, 詳しくは参考文献をご覧ください.

【不正確な定義】

確率変数列\(\{X_n\}_{n=0}^{\infty}\)が遷移核\(K\)をもつMarkov連鎖であるとは, 任意の自然数\(n\), 任意の\(A\)に対して, 次式を満たすことである.

\begin{align} P\left( X_{n+1}\in A\mid X_n,\cdots,X_0 \right) &= P\left( X_n\in A\mid X_n=x\right) \\ &= \int_{A} K(x,y)dy \end{align}

上式を, \(K(x,A)\)で表す.

定義から分かるように, Markov連鎖に関する情報は遷移核に集約されています. したがって, MCMCの基本方針をより詳しく言えば, 目的の分布を定常分布とするようなMarkov連鎖の遷移核を構成することです. 遷移核は連鎖の移動ルールを確率的に定めたものです. そして, MCMCの各手法の違いは遷移核の違いです. この遷移核をうまく構成することで, 1つのMCMCアルゴリズムを得ることができます.

定常分布と極限分布

前節では, MCMCの設計指針について概観しました. ここでは, MCMCのユーザ側の視点に立って考えてみましょう. データの解析に際して, 設計されたMCMCアルゴリズムは信頼に足るものでしょうか. ユーザー側が信頼性を確認するにはどうすれば良いでのしょうか. まず, 手がかりとなる定常分布を定義します.

【不正確な定義】

分布\(\nu\)は, 任意の\(A\)に対して次式を満たすとき, 不変分布であるという. \begin{equation} \nu (A) = \int K(x,A)\nu (dx) \end{equation}

定常分布と似た概念に, 極限分布というものもあります. 定義は以下の通りです.

【不正確な定義】

分布\(\nu\)は, 任意の\(A\), 任意の\(x\)に対して次式を満たすとき, 極限分布であるという. \begin{equation} \nu (A) = \lim_{n\to \infty} K^n(x,A) \end{equation}

実は, ある条件の元で定常分布と極限分布は一致します. 一致すると, MCMCのユーザーにとって何が嬉しいのでしょうか. 極限分布は遷移核を何度も作用させた結果得られる分布です. 遷移核は移動のルールを定めたものですから, 極限分布とは, 同じルールを何度も適用した結果です. もし遷移核の構成方法(移動ルール)をアルゴリズムとして記述できれば, 反復することで定常分布が得られることになります. そもそも定常分布が目的の分布と一致するように遷移核を構成していました. したがって, 反復の結果として目的の分布が得られます. つまり, MCMCにより「目的の分布を定常分布とするMarkov連鎖」の極限分布として, 目的の分布を得ることができます. 手計算で扱うことが難しい分布であっても, 移動のルールをうまく定めれば計算機で計算できる. これがMCMCの威力です. ユーザー側としては, 「ある条件」とやらを確認しておけば安心です. その条件を次節で紹介します.

収束性

本節で, 定常分布と極限分布が(存在するとして,)一致するための十分条件を提示します. 詳しくは参考文献を参照してください. 一致するための条件は既約性と非周期性です.

【不正確な定義】

不変分布\(\nu\)をもつMarkov連鎖は, 任意の初期状態に対して, \(\nu\)が正である集合に属する確率が正であるとき, \(\nu -\)既約であるという.

Markov連鎖は, ある一定間隔で訪れる状態が存在するとき, 周期的であるという. 周期的でないとき, 非周期的であるという.

既約性は, どこから出発しても, 確率が正の場所へたどり着くことができる, ということです. 満遍なくサンプリングしたいなら当然の要請です. 非周期性の要請も, 偏ったサンプルを避けるためです.

【定理】

Markov連鎖が不変分布\(\nu\)をもち, \(\nu-\)既約かつ非周期的であるとき, 不変分布は一意的であり, 極限分布は不変分布に一致する.

以上より, MCMCの個々のアルゴリズムは, 以下の4つを満たしていればとりあえずOKです.

  • アルゴリズムによって得られる数またはベクトル列は, Markov連鎖である.
  • そのMarkov連鎖は目的の分布を定常分布とする.
  • そのMarkov連鎖は既約である.
  • そのMarkov連鎖は非周期的である.

詳細釣り合い条件

「MCMCアルゴリズムによって生成されたMarkov連鎖が, 目的の分布を定常分布とする」ことを保証する, 十分条件が存在します. これを詳細釣り合い条件といいます. MCMCアルゴリズムが詳細釣り合いを満たしてくれれば安心です.

【定理】

遷移核\(K\)に対して, 関数\(p\)が任意の\(x,y\)に対して, \begin{equation} K(y,x)p(y) = p(x)K(x,y) \end{equation} を満たすとき, \(p\)はこのMarkov連鎖の定常分布の密度関数である.

具体的な手法

Metropolis-Hastings法

Metropolis-Hastings法(以下, MH法)はMCMCの中でも基本的なアルゴリズムです. 導出過程の簡単な説明が[]にあります. アルゴリズムは以下の通りです.

【Metropolis-Hastings法(1変数)】

Initialize \(x^{(0)}\).
for \(s=1,\cdots,S-1\)
Sample \(x^{\prime}\) from \(q(x^{\prime}\mid x^{(s-1)})\).
Set \(\alpha = \min\left\{ 1,\frac{p(x^{\prime})q(x^{(s-1)}\mid x^{\prime})}{p(x^{(s-1)})q(x^{\prime}\mid x^{(s-1)})} \right\} \).
Sample \(u\) from uniform distribution.
if \(u \leq \alpha\)
Set \(x^{(s)}=x^{\prime}\).
else
Set \(x^{(s)}=x^{(s-1)}\).
end if
end for

MH法では, 目的の分布をよく近似するような提案分布\(q(x|y)\)を用います. \(y\)の部分にサンプルを代入して, 前回のサンプル値に依存したサンプルを提案します. 受容確率\(\alpha\)によって, 遷移のバランスをとるよう調節できます. \(\alpha\)の具体的な形は, 詳細釣り合い条件から得られます([7]). 提案分布の選び方に関しては様々な研究があります. 興味があれば調べてみてください.

Gibbsサンプラー

Gibbsサンプラーも基本的なアルゴリズムです. サンプルを行いたい変数が複数あるとき, 他を全て固定して, 条件付き分布から1変数ずつサンプルします. 条件付き分布からのサンプリングが簡単な場合に有効です. 階層モデルなど, 条件付き分布を用いるモデルではよく用いられます. アルゴリズムは以下の通りです. 以下のアルゴリズムでは, 変数を任意に2つのブロックに分けています. サンプリングしやすい分け方を勝手に採用すればOKです.

【Gibbsサンプラー(多変数)】

Initialize \(\boldsymbol{x}^{(0)} = \begin{bmatrix}\boldsymbol{x}_1^{(0)}\\ \boldsymbol{x}_2^{(0)}\end{bmatrix}\).
for \(s=1,\cdots,S-1\)
Sample \(\boldsymbol{x}_1^{(s)}\) from \(p(\boldsymbol{x}_1 \mid \boldsymbol{x}_2^{(s-1)})\).
Sample \(\boldsymbol{x}_2^{(s)}\) from \(p(\boldsymbol{x}_2 \mid \boldsymbol{x}_1^{(s)})\).
Set \(\boldsymbol{x}^{(s)} = \begin{bmatrix}\boldsymbol{x}_1^{(s)}\\ \boldsymbol{x}_2^{(s)}\end{bmatrix}\).
end for

Gibbsサンプラーのアルゴリズムは素直です. 全ての変数を動かすのが難しいなら, 他を全部固定して1つだけ動かせば良いわけです. また, GibbsサンプラーはMH法の特別バージョンです. MH法として見るならば, 受容確率は1です.

Hamiltonian Monte Carlo法

HMCは少々複雑ですが, 効率よくサンプルできます. Gibbsサンプラーとの比較を例2で扱います. アルゴリズムは以下の通りです.

【Hamiltonian Monte Carlo法(1変数)】

Initialize \(x^{(0)}\) and set step size \(h\) and terminated time \(T\) of leap-frog.
for \(s=1,\cdots,S-1\)
Sample \(v\) from \(\mathrm{N}(0,1)\).
Compute \((x^{\prime},v^{\prime}) = \mathrm{leapfrog}(x^{(s-1)}, v)\).
Compute \(\Delta H = H(x^{\prime}, v^{\prime})-H(x^{(s-1)}, v)\).
Set \(\alpha = \min\left\{ 1,e^{-\Delta H}\right\} \).
Sample \(u\) from uniform distribution.
if \(u \leq \alpha\)
Set \(x^{(s)}=x^{\prime}\).
else
Set \(x^{(s)}=x^{(s-1)}\).
end if
end for

まず, 基本的なアイデアを述べます.

次のような密度関数で表される分布からのサンプルが欲しいとします.

\begin{equation} p(x) = \frac{1}{Z}\tilde{p}(x) \end{equation}

このとき, \(E(x)=-\log \tilde{p(x)}\)とおくと, \begin{equation} p(x) = \frac{1}{Z}\mathrm{exp}\left\{-E(x)\right\} \end{equation} と書けます. ここで, 新しい変数\(v\)を導入して, \begin{equation} H(x,v) = \frac{1}{2}\|v\|^2 + E(x) \end{equation} と定めます. よって, \begin{equation} p(x) \propto \mathrm{exp}\left\{-\left( H(x,v)-\frac{1}{2}\|v\|^2 \right) \right\} = \mathrm{exp}\left\{\frac{1}{2}\|v\|^2\right\}\mathrm{exp}\left\{-H(x,v)\right\} \end{equation} となります. よって, サンプル\(\{x^{(s)}, v^{(s)}\}\)が得られれば, \(v\)のサンプルを無視して, サンプル\(\{x^{(s)}\}\)が得られます. 一見ややこしいように見えます. しかし, \(H\)を導入することで効率がよくなります. 鍵となるのは次のHamilton系です.

【Hamilton系】

次のHamiltonian \begin{equation} H(x,v) = \frac{1}{2}\|v\|^2 + E(x) \end{equation} に対応するHamilton系は, 次式で表される.

\begin{align} \frac{dx}{d\tau} &= p\\ \frac{dv}{d\tau} &= -\nabla E(x) \end{align}

ただし, \(0\leq \tau \leq T\)とする. 以下, Tを終端時間と呼ぶ.

このように, \(H\)をHamilton系と関連付けます. 物理的な系と解釈することで, 微分方程式のシミュレーションが可能になります. 微分方程式を解くことで, より遠くの状態にまでアクセスすることができます. もちろん, 微分方程式によりアクセスできる状態は一部分です. そこで, Hamiltonianに確率的な揺らぎを加えて, 反復ごとに異なるHamilton系をシミュレートします. これが, アルゴリズム冒頭の, 標準正規分布からのサンプルです. 変数\(v\)は, Hamiltonianの値を変える役割を担います. 問題は, どうやってHamilton系をシミュレートするかです.

微分方程式を解くのに, 数値解法の一種, leap-frog法(Störmer-Verlet法とも呼ぶ)を用います. 疑似コードでは関数として与えています. なぜleap-frog法が選ばれたのか. 理由は, leap-frog法がHamilton系の性質を保ったまま計算してくれるからです. ここでいうHamilton系の性質とは, 以下の3つです.

  • エネルギー保存性
  • 体積保存性
  • 時間反転性

エネルギー保存性については, 「微分方程式の数値解法」の節で解説します. (その他の性質についてはまた実験しようと思います.)

leap-frog法の更新式は以下の通りです. これを適当な刻み幅\(h\), 適当な終端時間\(T\)まで反復して, 次の候補を選びます. これが関数leap-frogです. 記号に関しては, 「微分方程式の数値解法」の節を参照してください.

【leap-frog法(Störmer-Verlet法)の1ステップ】

\begin{align} v_{n+\frac{1}{2}} &= v_n - \frac{h}{2}\nabla E(x_n) \\ x_{n+1} &= x_n + h v_{n+\frac{1}{2}} \\ v_{n+1} &= v_{n+\frac{1}{2}} - \frac{h}{2} \nabla E(x_{n+1}) \end{align}

技術的な注意点

burn-in

MCMCでは極限分布に収束するまでに時間がかかります. 最初のうちは, 初期値の影響を受けており, 定常分布に収束していません. この部分を定常分布からのサンプルと見做すのは危険です. このように, 初期値の影響を受けて, 定常分布に収束していない期間をburn-in期間といいます. 実際にMCMCを使う際には, burn-in期間のサンプルは捨てます. ここで, どれくらいの期間をbrun-in期間と見做すか, という問題が生じます. これはケースバイケースです. 個人的には, 全サンプル数の10分の1を捨てることにしています. 全部で10,000サンプルあれば, 1,000サンプルを捨てます. 残りの9,000サンプルを定常分布からの実現値と考えます.

収束判定

MCMCでは極限分布を計算します. 計算機の上ではサンプル数を無限大にすることはできないので, 収束を判定する指標が欲しいところです. 代表的なものをいくつか挙げます. なお, 収束診断に関してはいくつかサーベイ論文が出ているので, 検索してみて下さい.

  • 時系列プロット
  • 標本自己相関関数のプロット
  • Gewekeの診断方法
  • Gelman-Rubinの方法

時系列プロットは, 横軸に時間(反復数), サンプル値を縦軸にとって系列をプロットする方法です. うまくサンプルできていれば, 縦にギザギザした特徴的なプロットが得られます. 同じ部分を何度もサンプルしたり, 規則性があるなど, 何か異常があればすぐに分かります. 一方, 一度にプロットできるのは1変数分だけという制約があります. 多変数の場合は別の方法が必要です. 以下は, 例1(後ほど登場)で用いた時系列プロットです.

【コード10の実行結果】

標本自己相関関数は次のように定義されます.

【定義】

時系列データ\(\{x_k\}_{k=1}^n\)に対して, 標本自己共分散関数を次のように定義する. \(\tau\)をラグという.

\begin{equation} \hat{C}_\tau = \frac{1}{n-\tau} \sum_{k=\tau+1}^{n} (x_k-\hat{\mu})(x_{k-\tau}-\hat{\mu})\quad (\tau=0,1,\cdots,n-1) \end{equation}

ここで, \(\hat{\mu}\)は標本平均である. このとき, 標本自己相関関数を, 次式で定める.

\begin{equation} \hat{\rho}_\tau = \frac{C_\tau}{C_0} \end{equation}

自己相関関数を見ると, 数個前のサンプルの影響をどの程度受けているかを確認できます. 一般に, 自己相関が早く減衰するほど, 定常分布への収束は早いです. これも, プロットして確認します. 縦軸に自己相関の値, 横軸にラグをとります. また, 一度にプロットできるのは1変数分だけという制約があります. 以下は, 例1(後ほど登場)で用いた自己相関関数です. なお, サンプル間の相関を現象させるために, 適当に間引くこともあります.

【コード10の実行結果】

Gelman-Rubinの方法は, 複数のMarkov連鎖を計算する必要があり, 計算コストがかかります. かといって一本当たりのMarkov連鎖を短くすると, Gibbsサンプラーのような収束の遅いアルゴリズムでは十分なサンプルが得られない危険性があります. ここでは, Markov連鎖1本で計算できる, Gewekeの方法を紹介します.

【Gewekeの診断方法】

サンプル列\(\{x_k\}_{k=1}^{n}\)に対して, 次のように定める.

\begin{equation} n_{10}=\frac{n}{10},\quad n_{50}=\frac{n}{2},\quad \hat{\mu}_{10}=\frac{1}{n_{10}}\sum_{k=1}^{n_{10}}x_k,\quad \hat{\mu}_{50}=\frac{1}{n_{50}}\sum_{k=n_{50}+1}^{n}x_k, \end{equation}

次に, 次の統計量を計算する.

\begin{equation} Z = \frac{\hat{\mu}_{10}-\hat{\mu}_{50}}{\sqrt{\frac{\hat{S}_{10}(0)}{n_{10}}+\frac{\hat{S}_{50}(0)}{n_{50}}}} \end{equation}

ここで, \(\hat{S}_{10}(0)\), \(\hat{S}_{50}(0)\)はそれぞれ, \(\{x_k\}_{k=1}^{n_{10}}\)と\(\{x_k\}_{k=n_{50}+1}^{n}\)のスペクトル密度推定値である.

\(z_{\alpha}\)を標準正規分布の\(100\alpha\)%点とするとき, \[ |Z| \leq z_{\frac{\alpha}{2}} \] ならば, サンプル列は収束していると判断する.

定理からわかるように, スペクトル計算を伴います. スペクトル計算にはいくつか選択肢があります. 個人的にはYule-Walker法を使っています. 手順を解説します. まず, MCMCにより得られた系列を, ARモデルでモデル化します. ARモデルは, 次のようなモデルです.

【定義】

次式を満たす時系列\(\{x_t\}_{t=1}^{\infty}\)を\(m\)次のARモデルという.

\begin{equation} x_t = \sum_{k=1}^{m} a_{k}x_{t-k} + \epsilon_t,\quad \epsilon_t\sim \mathrm{N}(0,\sigma^2) \end{equation}

係数\(\{a_k\}_{k=1}^{m}\)をAR係数という. ノイズ\(\epsilon_t\)は互いに無相関で, 過去の時系列データと独立とする.

次数を予め決めておく必要があります. AICを用いたモデル選択もできます. 個人的には次数を2に設定しています. このモデルに従う限り, サンプル系列のスペクトルは以下のように表されます.

【定理】

上記ARモデルのスペクトル密度は次のように表される. \begin{equation} S(f) = \frac{\sigma^2}{|1-\sum_{k=1}^{m} a_j e^{-2\pi i k f}|} \end{equation}

ARモデルに従うと仮定はしたものの, 肝心のAR係数とノイズ分散は未知です. そこでAR係数とノイズ分散を求めます. 求め方の1つがYule-Walker法です.

【Yule-Walker法】

上記\(m\)次のARモデルのAR係数の推定値は, 次式を解くことで得られる.

\begin{equation} \begin{bmatrix} \hat{C}_0 & \hat{C}_1 & \cdots & \hat{C}_{m-1} \\ \hat{C}_1 & \hat{C}_0 & \cdots & \hat{C}_{m-2} \\ \vdots & \vdots & \ddots & \vdots \\ \hat{C}_{m-1} & \hat{C}_{m-2} & \cdots & \hat{C}_{0} \\ \end{bmatrix} \begin{bmatrix} \hat{a_1} \\ \hat{a_2} \\ \vdots\\ \hat{a_m} \\ \end{bmatrix} = \begin{bmatrix} \hat{C}_1 \\ \hat{C}_2 \\ \vdots\\ \hat{C}_m \\ \end{bmatrix} \end{equation}

また, ノイズ分散推定値は次式で求められる.

\begin{equation} \hat{\sigma}^2 = \hat{C}_0 - \sum_{k=1}^m \hat{a}_k \hat{C}_k \end{equation}

この連立一次方程式を解くことで, AR係数, およびノイズ分散を求めることができます. これらの値が分かれば, スペクトルを算出できます. 特に周波数0のときの値が知りたいので, 次のような代入計算を行います.

\begin{equation} \hat{S}(0) = \frac{\hat{\sigma}^2}{|1-\sum_{k=1}^{m} \hat{a_k}|} \end{equation}

これでGewekeの診断方法が使えます. まとめると, 「サンプル系列がARモデルに従うと仮定する」→「標本分散共分散行列を求める」→「Yule-Walker法により, AR係数およびノイズ分散の推定値を求める」→「スペクトル密度の推定値を求める」→「Gewekeの方法の通りに計算する」.

微分方程式の数値解法

HMC法で用いたleap-frog法に関して補足しておきます. leap-frog法はなぜ必要なのでしょうか. leap-frog法は微分方程式の数値解法の1種です. まずは解法を整理しておきましょう. 次の1階自励系微分方程式を考えます.

\begin{equation} \frac{d}{dt}x = f(x),\quad x:\mathbb{R}\to \mathbb{R}^d,\quad f:\mathbb{R}^d\to \mathbb{R}^d \end{equation}

このような微分方程式の数値解法として有名なものを2つほど導入します. 以下, \(h\)を離散化の幅, \(t_n=nh\)とし, \(x_n\)を\(x(t_n)\)の近似値とします. 1つ目は, 陽的Euler法です.

【スキーム1: 陽的Euler法】

\begin{equation} x_{n+1} = x_n + hf(x_n) \end{equation}

2つ目は4段4次Runge-Kutta法です.

【スキーム2: 4段4次Runge-Kutta法】

\begin{equation} x_{n+1} = x_n + \frac{h}{6}(k_1 + 2k_2 + 2k_3 + k_4) \end{equation}

ただし, \begin{align} k_1 &= f(x_n) \\ k_2 &= f\left(x_n + \frac{h}{2}k_1\right) \\ k_3 &= f\left(x_n + \frac{h}{2}k_2\right) \\ k_4 &= f(x_n + hk_3) \end{align} とする.

微分方程式の数値解法の適用範囲を広げるため, 階数を1つあげます. つまり, 2階の微分方程式も考えることにします.

\begin{equation} \frac{d^2}{dt^2} q = f(q) \end{equation}

ここで, 次のような変数を導入します.

\begin{equation} v = \frac{d}{dt}q \end{equation}

このとき, 上の2階常微分方程式は以下のように書けます.

\begin{equation} \frac{d}{dt}\left[\begin{array}{c}q\\v\end{array}\right]  = \left[\begin{array}{c}v\\f(q)\end{array}\right] \end{equation}

このように, 2階の微分方程式であっても, 1階に帰着します. したがって, 先ほどの数値解法, 陽的Euler法, 4段4次Runge-Kutta法も適用できます. さらに2つほど数値解法を導入します. 3つ目の数値解法が, symplectic Euler法です.

【スキーム3: symplectic Euler法】

\begin{align} q_{n+1} &= q_n + hv_n \\ v_{n+1} &= v_n + hf(q_{n+1}) \\ \end{align}

4つ目の数値解法が, leap-frog法(Störmer-Verlet法)です.

【スキーム4: leap-frog法(Störmer-Verlet法)】

\begin{align} v_{n+\frac{1}{2}} &= v_n + \frac{h}{2}f(q_n) \\ q_{n+1} &= q_n + h v_{n+\frac{1}{2}} \\ v_{n+1} &= v_{n+\frac{1}{2}} + \frac{h}{2} f(q_{n+1}) \end{align}

数値解法を4つほど導入したところで, これらを比較します. 例題として, Kepler問題を考えます.

【Kepler問題】

次の微分方程式系を考える.

\begin{equation} \ddot{q_1} = -\frac{q_1}{(q_1^2 + q_2^2)^{\frac{3}{2}}}, \quad \ddot{q_2} = -\frac{q_2}{(q_1^2 + q_2^2)^{\frac{3}{2}}} \\ \end{equation}

この系の解の様子をそれぞれの解法について, プロットして確かめる. さらに, この系はHamilton系である. そこで, Hamiltonianの値の変化も確かめる.

Kepler問題も1階の系に帰着します. 次のようにおきます.

\begin{equation} q = \left[ \begin{array}{c} q_1 \\ q_2 \end{array} \right], \quad v = \left[ \begin{array}{c} \dot{q_1} \\ \dot{q_2} \end{array} \right] = \frac{d}{dt}q \end{equation}

よって, 次のように書けます.

\begin{equation} \frac{d}{dt}q = v, \quad \frac{d}{dt}v = -\frac{q}{\|q\|^3} \end{equation}

以上より, \begin{equation} \frac{d}{dt} \left[ \begin{array}{c} q \\ v \end{array} \right] = \left[ \begin{array}{c} v \\ f(q) \end{array} \right], \quad f(q) = -\frac{q}{\|q\|^3} \end{equation} と書けます. 上式に様々な解法を適用します. また, Hamiltonianは, \begin{equation} H(q,p) = \frac{1}{2}\|v\|^2 - \frac{1}{\|q\|} \end{equation} で与えられます. 下図は, ステップ数1,000で, ステップ幅0.05とした場合の軌道のアニメーションです. ただし, 表示しているのは最初の300ステップ分です. 初期値は, 以下のように設定しました.

\begin{equation} q_0 = \begin{bmatrix} 0.4 \\ 0.0 \end{bmatrix},\quad v_0 = \begin{bmatrix} 0.0 \\ 2.0 \end{bmatrix} \end{equation}

水色は陽的Euler法, 橙色は4段4次Runge-Kutta法, 緑色はsymplectic Euler法, 紫色はleap-frog法(Störmer-Verlet法)です. Euler法の軌道はかなり外れています. その他の解法は一定の軌道を描いています.

【コード8の実行結果】

Hamiltonianの誤差変化も調べます. 長時間軌道計算すると, Runge-Kutta法の誤差は大きくなります. 一方, symplectic Euler法とleap-frog法(Störmer-Verlet法)は誤差が一定です. 特に, 長時間計算した際に, 顕著な差が現れます. これらの解法は, エネルギー保存というHamilton系の性質を保っています(エネルギー保存性). この性質のおかげで, HMCでは(微分方程式内の仮想的な時間で)長時間計算しても, 棄却率を低く保つことに成功しています(\(\Delta H\)を非常に小さくできる). したがって, HMC内のleap-frog法の終端時間\(T\)を大きめに設定しておけば, 棄却率を低く保ちながら遠くに移動することができます. これらの理論的結果に関しては, [8]が詳しいです.

【コード7の実行結果】

MCMCの実装

例1: 簡易的なMetropolis-Hastings法の実験

1つ目の例です. ここでは以下の混合正規分布からサンプリングを行います.

\begin{equation} 0.6\mathrm{N}(-0.2,1.5) + 0.4\mathrm{N}(2.0,1.5) \end{equation}

提案分布は以下の正規分布とします.

\begin{equation} \mathrm{N}(x,2),\quad x\text{は前反復のサンプル} \end{equation}

サンプル数は10,000で, そのうちburn-in期間として1,000サンプルは除去します.

結果は下図の通りです. ヒストグラム, 時系列プロット, 自己相関関数共に大丈夫そうです.

【コード10の実行結果】

例2: GibbsサンプラーとHMCの比較アニメーション

2つ目の例です. ここではGibbsサンプラーとHMCの動きを比較します. 分散共分散行列が \[ \Sigma = \begin{bmatrix} 1.0 & \rho \\ \rho & 1.0 \end{bmatrix}, \quad \rho=0.99 \] の2次元正規分布からのサンプリングを行います. 1,000個のサンプルをとり, 100サンプルはburn-in期間として捨てます. 原点を初期値として, 収束診断はGewekeの方法を採用します. HMC内のleap-frog法の終端時間は100, ステップ幅は0.01とします. したがって, 10,000ステップ分計算しています.

以下のアニメーションを見れば分かるように, Gibbsサンプラーでは十分にサンプルできません. Gibbsサンプラーは一方の変数を固定して, 他方の変数を動かすため, 動きが縦と横方向のみに限定されます. このように, 縦と横の幅が狭い分布では違いが顕著に現れます.

Gibbsサンプラーの法は, 次のような正規分布からサンプルを行います.

\begin{align} x_1\mid x_2 &\sim \mathrm{N}\left(\rho x_2, \sqrt{1-\rho^2}\right) \\ x_2\mid x_1 &\sim \mathrm{N}\left(\rho x_1, \sqrt{1-\rho^2}\right) \\ \end{align}

また, HMCの方は, 次のようにします.

\begin{equation} E(x)=\frac{x_1^2-2\rho x_1x_2+x_2^2}{2(1-\rho^2)} \end{equation}
【コード13の実行結果】
【コード15の実行結果】

ソースコード

インポート

【Juliaコード1; インポート】
using Plots
using Random
using Statistics
using Distributions
using LinearAlgebra
using Measures
pyplot()

微分方程式の数値解法

【Juliaコード2; 微分方程式の数値解法: 共通のコード】
#終端時間
T = 1000;

#ステップ幅
h = 0.05;

#ステップ数
n = Int(T/h);

#初期値
e = 0.6;
q₀ = [1-e; 0];
v₀ = [0;sqrt((1+e)/(1-e))];

#右辺の関数f
function f(q::Array{Float64,1})
    return -q ./norm(q)^3
end

#Hamiltonianの計算
function H(q::Array{Float64,1}, v::Array{Float64,1})
    return norm(v)^2/2 - 1/norm(q)
end

#hamiltonianの真値
H_true = H(q₀, v₀);
【Juliaコード3; 陽的Euler法】
#保存用配列
qvec_E = zeros(2, n);
vvec_E = zeros(2, n);

#初期値の保存
q_old_E = q₀;
v_old_E = v₀;
qvec_E[:,1] = q_old_E;
vvec_E[:,1] = v_old_E;

#陽的Euler法の1ステップ
function myEuler(q_old::Array{Float64,1}, v_old::Array{Float64,1}, h::Float64)
    q_new = q_old + h*v_old;
    v_new = v_old + h*f(q_old);
    return q_new, v_new
end

#変数の初期化
q_new_E, v_new_E = q_old_E, v_old_E;

#反復計算
for k in 1:n-1
    #q_oldとv_oldの参照
    q_old_E = qvec_E[:,k];
    v_old_E = vvec_E[:,k];

    #qとvの更新
    q_new_E, v_new_E = myEuler(q_old_E, v_old_E, h);

    #更新値の保存
    qvec_E[:,k+1] = q_new_E;
    vvec_E[:,k+1] = v_new_E;
end

#Hamiltonianの計算
H_Euler = [H(qvec_E[:,k], vvec_E[:,k]) for k in 1:n];
【Juliaコード4; 4段4次Runge-Kutta法】
#保存用配列
qvec_RK = zeros(2, n);
vvec_RK = zeros(2, n);

#初期値の保存
q_old_RK = q₀;
v_old_RK = v₀;
qvec_RK[:,1] = q_old_RK;
vvec_RK[:,1] = v_old_RK;

#4段4次RK法の1ステップ
function myRK(q_old::Array{Float64,1}, v_old::Array{Float64,1}, h::Float64)
    #kの更新
    kq₁ = v_old;
    kv₁ = f(q_old);

    kq₂ = v_old + h*kv₁/2;
    kv₂ = f(q_old + h*kq₁/2);

    kq₃ = v_old + h*kv₂/2;
    kv₃ = f(q_old + h*kq₂/2);

    kq₄ = v_old + h*kv₃;
    kv₄ = f(q_old + h*kq₃);

    #qとvの更新
    q_new = q_old + h*(kq₁+2*kq₂+2*kq₃+kq₄)/6;
    v_new = v_old + h*(kv₁+2*kv₂+2*kv₃+kv₄)/6;
    return q_new, v_new
end

#変数の初期化
q_new_RK, v_new_RK = q_old_RK, v_old_RK;

#反復計算
for k in 1:n-1
    #q_oldとv_oldの参照
    q_old_RK = qvec_RK[:,k];
    v_old_RK = vvec_RK[:,k];

    #qとvの更新
    q_new_RK, v_new_RK = myRK(q_old_RK, v_old_RK, h);

    #更新値の保存
    qvec_RK[:,k+1] = q_new_RK;
    vvec_RK[:,k+1] = v_new_RK;
end

#Hamiltonianの計算
H_RK = [H(qvec_RK[:,k], vvec_RK[:,k]) for k in 1:n];
【Juliaコード5; Symplectic Euler法】
#保存用配列
qvec_SE = zeros(2, n);
vvec_SE = zeros(2, n);

#初期値の保存
q_old_SE = q₀;
v_old_SE = v₀;
qvec_SE[:,1] = q_old_SE;
vvec_SE[:,1] = v_old_SE;

#Symplectic Euler法の1ステップ
function mySymplecticEuler(q_old::Array{Float64,1}, v_old::Array{Float64,1}, h::Float64)
    q_new = q_old + h*v_old;
    v_new = v_old + h*f(q_new);
    return q_new, v_new
end

#変数の初期化
q_new_SE, v_new_SE = q_old_SE, v_old_SE;

#反復計算
for k in 1:n-1
    #q_oldとv_oldの参照
    q_old_SE = qvec_SE[:,k];
    v_old_SE = vvec_SE[:,k];

    #qとvの更新
    q_new_SE, v_new_SE = mySymplecticEuler(q_old_SE, v_old_SE, h);

    #更新値の保存
    qvec_SE[:,k+1] = q_new_SE;
    vvec_SE[:,k+1] = v_new_SE;
end

#Hamiltonianの計算
H_symplecticEuler = [H(qvec_SE[:,k], vvec_SE[:,k]) for k in 1:n];
【Juliaコード6; Störmer-Verlet法 (Leapfrog法)】
#保存用配列
qvec_SV = zeros(2, n);
vvec_SV = zeros(2, n);

#e
e = 0.6;

#初期値の保存
q_old_SV = q₀;
v_old_SV = v₀;
qvec_SV[:,1] = q_old_SV;
vvec_SV[:,1] = v_old_SV;

#Störmer-Verlet法の1ステップ
function myStörmerVerlet(q_old::Array{Float64,1}, v_old::Array{Float64,1}, h::Float64)
    v_mid = v_old + h*f(q_old)/2;
    q_new = q_old + h*v_mid;
    v_new = v_mid + h*f(q_new)/2;
    return q_new, v_new
end

#変数の初期化
q_new_SV, v_new_SV = q_old_SV, v_old_SV;

#反復計算
for k in 1:n-1
    #q_oldとv_oldの参照
    q_old_SV = qvec_SV[:,k];
    v_old_SV = vvec_SV[:,k];

    #qとvの更新
    q_new_SV, v_new_SV = myStörmerVerlet(q_old_SV, v_old_SV, h);

    #更新値の保存
    qvec_SV[:,k+1] = q_new_SV;
    vvec_SV[:,k+1] = v_new_SV;
end

#Hamiltonianの計算
H_StörmerVerlet = [H(qvec_SV[:,k], vvec_SV[:,k]) for k in 1:n];
【Juliaコード7; Hamiltonianの比較】
p_Herror = plot([
    abs.(H_Euler.-H_true),
    abs.(H_RK.-H_true),
    abs.(H_symplecticEuler.-H_true),
    abs.(H_StörmerVerlet.-H_true)],
    title="Error of H(q,v)",
    xlabel="time",
    ylabel="abs error",
    xscale=:log10,
    yscale=:log10,
    labels=["E" "RK" "SE" "SV"],
    ylim=[1e-6,1e0]
);
#savefig(p_Herror, "figs-oed/Herror.png");
【Juliaコード8; 軌道のアニメーション】
plot(xlabel="q₁", ylabel="q₂", xlim=[-2,1.5], ylim=[-1,1.5], legends=false);
anim = @animate for t in 1:300
    plot!(qvec_E[1,1:t], qvec_E[2,1:t], color=:skyblue);
    plot!(qvec_RK[1,1:t], qvec_RK[2,1:t], color=:orange);
    plot!(qvec_SE[1,1:t], qvec_SE[2,1:t], color=:green);
    plot!(qvec_SV[1,1:t], qvec_SV[2,1:t], color=:purple);
end
#gif(anim, "figs-oed/orbit.gif");

例1のコード

【Juliaコード9; 例1のMH法】
#1つ目の正規分布
μ₁ = -2.0
σ₁ = sqrt(1.5)
normal1 = Normal(μ₁, σ₁)

#2つ目の正規分布
μ₂ = 2.0
σ₂ = sqrt(1.5)
normal2 = Normal(μ₂, σ₂)

#混合正規分布
c = 0.6
p(x) = c*pdf(normal1,x) + (1-c)*pdf(normal2,x)

#自己相関関数
function auto_corr_func(samps::Array{Float64}, τ::Int64)
    #平均と分散
    μ = mean(samps)
    C₀ = var(samps)

    #共分散関数
    Ck = mean((samps[1:end-τ].-μ) .* (samps[τ+1:end].-μ))
    return Ck/C₀
end

#提案分布
q(x,x_samp) = pdf(Normal(x_samp, sqrt(2)), x)

#MH法
function myMH(init::Float64, n_samps::Int64, n_burnin::Int64)
    #サンプル保存用配列と初期値の保存
    x_samps = zeros(n_samps)
    x = init
    x_samps[1] = x

    #カウンター
    counter = 0

    for i in 2:n_samps
        x_cand = x+sqrt(2)*randn()
        α = p(x_cand)*q(x,x_cand)/p(x)/q(x_cand,x)
        if rand() < α
            #採択
            x = x_cand
            counter += 1
        end
        x_samps[i] = x
    end

    #burn-inを除去
    x_samps = x_samps[n_burnin:end]
    return x_samps, counter
end
【Juliaコード10; 例1のMH法の実行】
Random.seed!(42)
init = 0.0 #初期値
n_samps = 10000 #サンプル数
n_burnin = Int(n_samps/10) #burn-in
x_samps, counter = myMH(init, n_samps, n_burnin)

#自己相関の計算
auto_corr(τ) = auto_corr_func(x_samps, τ)

#ヒストグラム, 時系列プロット, 自己相関関数
true_pdf_plot = plot(p, -5, 5, xlabel="x", ylabel="prob density", label="true pdf")
histogram!(x_samps, bins=100, norm=true, alpha=0.5, label="mcmc sample", margin=5mm)
time_series_plot = plot(x_samps, lw=0.5, xlabel="iter", ylabel="x", label="mcmc sample", margin=5mm)
auto_corr_plot = bar(0:20, auto_corr, xlabel="lag", ylabel="auto corr", label=false, margin=5mm)
result_plot = plot(true_pdf_plot, time_series_plot, auto_corr_plot, layout=(3,1), size=(800, 1000))

#保存
#savefig(result_plot, "figs-MH/result.png")
#savefig(time_series_plot, "figs-MH/ts-plot.png")
#savefig(auto_corr_plot, "figs-MH/auto-corr-plot.png")

例2のコード

【Juliaコード11; 例2の共通のコード】
#相関係数
const ρ = 0.99;

#正規分布
x₁_arr = -2.5:0.05:2.5
x₂_arr = -2.5:0.05:2.5
Σ = [
    1.0 ρ;
    ρ 1.0
]
normal_dist = MvNormal(zeros(2), Σ)
normal2d(x₁,x₂) = pdf(normal_dist, [x₁;x₂])

#Yule-Waler方程式を解いて, AR係数を求める
function myYuleWalker(samps::Array{Float64,1})
#サンプルの長さ
L = length(samps);

#サンプル平均 
sample_mean = mean(samps);

#自己共分散関数C₀, C₁, C₂を求める
C₀ = var(samps);
C₁ = mean((samps[2:end] .- sample_mean).*(samps[1:L-1] .- sample_mean));
C₂ = mean((samps[3:end] .- sample_mean).*(samps[1:L-2] .- sample_mean));

#Yule-Walker方程式を解く
coeff_mat = [C₀ C₁; C₁ C₀];
rhs_vec = [C₁;C₂];
avec = coeff_mat\rhs_vec;

#ノイズ分散を求める
σsq = C₀ - dot(avec, [C₁,C₂]);

return avec, σsq, sample_mean
end

#Gewekeの診断方法(収束診断)
function myGeweke(samps::Array{Float64,1})
#サンプルの長さ
L = length(samps);

#sampsの前半10%の平均とスペクトルを求める
N10 = Int(floor(L/10));
samps10 = samps[1:N10];
a10, σ10_sq, mean10 = myYuleWalker(samps10);
spec10 = σ10_sq/(1-sum(a10))^2;

#sampsの後半50%の平均とスペクトルを求める
N50 = Int(floor(L/2));
samps50 = samps[N50+1:end];
a50, σ50_sq, mean50 = myYuleWalker(samps50);
spec50 = σ50_sq/(1-sum(a50))^2;

#診断
R = (mean10-mean50)/sqrt(spec10/N10 + spec50/N50);
is_convergence = abs(R) < 1.96

#もし収束していたらtrueを返す
return is_convergence, R
end
【Juliaコード12; 例2のGibbsサンプラー】
#Gibbs sampler
function myGibbs(init::Array{Float64,N}, maxsamps::Int64, n_burnin::Int64, diag=false, counter=false) where N
    #パラメータの次元
    d = length(init);

    #保存用配列
    params_array = zeros(d, maxsamps);

    #初期値の保存
    params_array[:,1] = init;

    #変数の初期化
    x₁_samp, x₂_samp = init[1], init[2];

    #サンプル
    for k in 1:maxsamps-1
        #μのサンプル
        x₁_samp = sqrt(1-ρ^2)*(randn()+ρ*x₂_samp);

        #λのサンプル
        x₂_samp = sqrt(1-ρ^2)*(randn()+ρ*x₁_samp);

        #保存
        params_array[:,k+1] = [x₁_samp, x₂_samp]';
    end

    #burn-in期間の除去
    params_array = params_array[:,n_burnin:end];

    #Gewekeの診断
    if diag == true
        for i in 1:d
            samps = params_array[i,:]
            is_convergence, R = myGeweke(samps);
            if is_convergence == false
                println("警告: $(i)番目のMCMCサンプルは収束していません. R=$(R).");
            end
        end
    end

    #カウンターの有無
    if counter == true
        return params_array, n_accept
    else
        return params_array
    end
end
【Juliaコード13; 例2のGibbsサンプラーの実行】
Random.seed!(42)
maxsamps = 1000 #サンプル数
n_burnin = Int(maxsamps/10) #バーンイン期間
gibbs_samps = myGibbs(zeros(2), maxsamps, n_burnin, true, false);

anim = @animate for k in 1:(maxsamps-n_burnin)
    p = plot(x₁_arr, x₂_arr, normal2d, st=:contour)
    scatter!(
        gibbs_samps[1,1:k],
        gibbs_samps[2,1:k],
        xlim=[-3,3],
        ylim=[-3,3],
        markercolor=:green,
        markerstrokealpha=0,
        alpha=0.4,
        label=false,
        xlabel="x₁",
        ylabel="x₂",
        title="Gibbs sampler"
    )
end
#gif(anim, "figs-mcmc/gibbs-anim.gif");
【Juliaコード14; 例2のHMC法】
#Störmer-Verlet法の1ステップ
function myStörmerVerlet(q_old::Array{Float64, 1}, v_old::Array{Float64, 1}, ϵ::Float64, neg∇E)
    v_mid = v_old .+ ϵ*neg∇E(q_old)/2;
    q_new = q_old .+ ϵ*v_mid;
    v_new = v_mid .+ ϵ*neg∇E(q_new)/2;
    return q_new, v_new
end

#HMC
function myHMC(init::Array{Float64,N}, ϵ::Float64, maxsteps::Int64, maxsamps::Int64, E, neg∇E, n_burnin::Int64, diag=false, counter=false) where N
    #パラメータの次元
    d = length(init);

    #サンプル保存用配列
    params_array = zeros(d,maxsamps);

    #初期値の保存
    params_array[:,1] = init;

    #採択された回数
    n_accept = 0;

    for k in 1:maxsamps-1
        #正規分布からのサンプル
        r₀ = randn(d);

        #パラメータ格納用変数
        params = params_array[:,k];

        #leap-frog法
        r = r₀;
        for _ in 1:maxsteps
            params, r = myStörmerVerlet(params, r, ϵ, neg∇E);
        end

        #採択または棄却
        α = minimum([1,exp(-E(params)-norm(r)^2/2)/exp(-E(params_array[:,k])-norm(r₀)^2/2)]);
        if rand() < α
            #採択
            params_array[:,k+1] = params;
            n_accept += 1;
        else
            #棄却
            params_array[:,k+1] = params_array[:,k];
        end
    end

    #burn-inを除去
    params_array = params_array[:, n_burnin:end];

    #Gewekeの診断
    if diag == true
        for i in 1:d
            samps = params_array[i,:]
            is_convergence, R = myGeweke(samps);
            if is_convergence == false
                println("警告: $(i)番目のMCMCサンプルは収束していません. R=$(R).");
            end
        end
    end

    #カウンターの有無
    if counter == true
        return params_array, n_accept
    else
        return params_array
    end
end
【Juliaコード15; 例2のHMC法の実行】
Random.seed!(42)
#ポテンシャル
E(x) = (x[1]^2-2*ρ*x[1]*x[2]+x[2]^2)/2/(1-ρ^2)
neg∇E(x) = -[
    (x[1]-ρ*x[2])/(1-ρ^2);
    (x[2]-ρ*x[1])/(1-ρ^2)
]

maxsamps = 1000 #サンプル数
n_burnin = Int(maxsamps/10) #バーンイン期間
ϵ = 0.01 #ステップ数
maxsteps = 100 #最大ステップ数
hmc_samps = myHMC(zeros(2), ϵ, maxsteps, maxsamps, E, neg∇E, n_burnin, true);

anim = @animate for k in 1:(maxsamps-n_burnin)
    p = plot(x₁_arr, x₂_arr, normal2d, st=:contour)
    scatter!(
        hmc_samps[1,1:k],
        hmc_samps[2,1:k],
        xlim=[-3,3],
        ylim=[-3,3],
        markercolor=:green,
        markerstrokealpha=0,
        alpha=0.4,
        label=false,
        xlabel="x₁",
        ylabel="x₂",
        title="HMC sampling"
    )
end
#gif(anim, "figs-mcmc/hmc-anim.gif");
参考文献
      [1] Robert, Casella, Monte Carlo Statistical Methods, Springer, New York, 2004
      [2] Neal, Bayesian Learning for Neural Networks, Springer, New York, 1996
      [3] 渡辺澄夫, ベイズ統計の理論と方法, コロナ社, 2012
      [4] 北川源四郎, Rによる時系列モデリング入門, 岩波書店, 2020
      [5] Tierney, Markov Chains for Exploring Posterior Distributions, The Annals of Statistics, 1994, 22-4, pp.1701-1728
      [6] Andrieu, Freitas, Doucet, Jordan, An Introduction to MCMC for Machine Learning, Machine Learning, 2003, 50, pp.5-43
      [7] Chib, Greenberg, Understanding the Metropolis-Hastings Algorithm, 1995, 49-4, pp.327-335
      [8] Bou-Rabee, Sanz-Serna, Geometric integrators and the Hamiltonian Monte Carlo method, Acta Numerica, 2018, 27, pp.113-206