Science

モンテカルロ法の基本をわかりやすく説明してみた(マルコフ連鎖モンテカルロ法まで)

モンテカルロ法の基本をわかりやすく説明してみた(マルコフ連鎖モンテカルロ法まで)

 

本記事では、モンテカルロ法の基本を初心者でも理解できるように説明しました。

具体的には、モンテカルロ法の基本的な考え方からマルコフ連鎖モンテカルロ法までをPythonによるシミュレーションを適宜使用しながら説明しました。

 

 

モンテカルロ法の基本知識

 

モンテカルロ法とは、乱数を用いた数値計算手法の総称のことを表します。

そうは言っても『乱数を用いた数値計算手法?』と思う方が多いと思うので、モンテカルロ法の基本的な考え方を具体例を通して説明していきます。

  1. モンテカルロ法による期待値評価

 

モンテカルロ法による期待値評価

 

ある\(P(\mathbf{X})\)に従う確率変数の期待値は、以下のように定義されます。

期待値の定義

  • 連続変数の期待値

$$\mathbb{E}[f(\mathbf{X})] = \int d\mathbf{X} f(\mathbf{X}) P(\mathbf{X})$$

  • 離散変数の期待値

$$\mathbb{E}[f(\mathbf{X})] = \sum_{\mathbf{X}} f(\mathbf{X}) P(\mathbf{X})$$

ここで、\(\sum_{\mathbf{X}}\cdots\)は全ての\(\mathbf{X}\)に関する和を取ることを意味します。

 

例えば、\(\mathbf{X} = (X_{1}, X_{2}, \ldots, X_{n})\)のように\(n\)次元ベクトルで表され、\(X_{i}\)が\(2\)種類の離散状態を取るケースでは、全ての状態に関する和は\(2^{n}\)通りになってしまい、計算するのが困難になります。

どのくらいのスピードで状態の数が増加するかというと以下のようなスピードで増加していきます。

計算量爆発ベクトルの次元と総状態数の関係

 

モンテカルロ法のアイデアは、もとの分布\(P(\mathbf{X})\)に従って生成されるサンプル列\(\{\mathbf{X}\}_{t=1}^{T}\)を生成し、そのサンプルの平均で期待値を置き換えることで期待値を近似的に評価しようというものです。

具体的には以下のように近似を行います。

モンテカルロ法による期待値評価

$$\mathbb{E}[f(\mathbf{X})]  \approx \sum_{t=1}^{T} f(\mathbf{X}^{t})$$

確率論の定理(中心極限定理)から、モンテカルロ法の期待値評価は、サンプル数\(T \to \infty\)の極限で期待値に収束します。

また、モンテカルロ法の期待値評価による誤差は、変数の次元に依存せずに減少していきます。

期待値評価問題を解析的に計算する問題から『確率分布からのサンプリング問題』に変換し、平均を取る問題に変換されました。

 

モンテカルロ法のイメージ

 

モンテカルロ法は、変数の確率分布\(P(\mathbf{X})\)からサンプル列\(\{\mathbf{X}^{t}\}_{t=1}^{T}\)を実際に生成し、そのサンプル列の平均を取ることで期待値を評価する手法でした。

いまいち、イメージが掴めない方のためにサイコロを例にして説明します。

サイコロは、\(\{1, 2, 3, 4, 5, 6\}\)の目の数を取ることを仮定し、サイコロの出る目(\(X\))の期待値を評価することを目指します。

サイコロの例の場合、解析的に期待値を計算でき、計算結果は以下のようになります。

$$1\times \frac{1}{6} + 2\times \frac{1}{6}+3\times \frac{1}{6}+4\times \frac{1}{6}+4\times \frac{1}{6}+5\times \frac{1}{6}+6\times \frac{1}{6}=\frac{7}{2}$$

 

まず、モンテカルロ法で近似的に出る目の期待値を評価するためには、出る目の確率分布(1〜6が出る一様分布)から実際にサンプルを生成します。

具体的に1000個のサンプルを得た場合、以下の式で期待値を近似的に評価するのがモンテカルロ法の基本的なアイデアです。

$$\mathbb{E}[X] \approx \frac{1}{100} \sum_{t=1}^{100} X^{t}$$

ここで、各\(X^{i}\)は実際に確率分布からサンプルされた値で、サイコロの場合は、1〜6の値になります。

具体的にPythonによるサイコロのシミュレーション結果を以下に示します。

サイコロのシミュレーション横軸は試行回数(サイコロを振った回数)、縦軸は得られたサンプル列による出た目の平均を意味する。

 

試行回数が多くなるにつれて、分散は小さくなり厳密な期待値である3.5に近づく様子が確認できますね。

サイコロの場合は、一様分布なので簡単に確率分布からのサンプルを得ることができました。

しかし、より複雑な確率分布\(P(\mathbf{X})\)からサンプルを得るためにはどうすれば良いでしょうか?

この方法を次に説明します。

 

確率分布からのサンプリング法

 

ここまでの話から、確率分布\(P(\mathbf{X})\)からのサンプルを十分な数、得ることできれば期待値を近似的に評価できることがわかりました。

具体的に確率分布\(P(\mathbf{X})\)からサンプルを得るためには主に二つの方法が使用されます。

  1. 直接的に独立サンプルを得る方法
  2. マルコフ連鎖モンテカルロ法

 

以下では詳しく説明します。

 

直接的に独立サンプルを得る方法

 

直接的に独立サンプルを得る方法は、一般にサンプルを得たい確率分布\(P(\mathbf{X})\)が低次元の確率分布の場合のみ有効で、高次元の確率分布からサンプルを得たい場合は、サンプル点を一個獲得する計算コストが多くなるというデメリットがあります。

具体的には、以下のような方法が主に使用されます。

  1. 逆関数法
  2. 棄却法
  3. 重点サンプリング

これらの手法の説明に関しては、気になる方は調べてみてください。

基本的に上記の方法で得られるサンプルは、サンプル列がお互い独立なので、中心極限定理と呼ばれる性質に従い誤差が減少していきます。

しかし、高次元の場合は、一つのサンプルを得るための計算時間が大幅に大きくなり、計算コスト的に問題が生じてしまいます。

そのため、独立サンプルを諦めてマルコフ連鎖を利用し、効率的にサンプル列を得ることを目的とするのが、次に説明するマルコフ連鎖モンテカルロ法(MCMC)です。

 

マルコフ連鎖モンテカルロ法

 

ここからは、マルコフ連鎖モンテカルロ法について説明していきます。

これまで説明していた『直接的に独立サンプルを得る方法』は変数の次元が大きくなると問題が生じてしまうことを説明しました。

マルコフ連鎖モンテカルロ法では、この問題点を解決するためにマルコフ連鎖の考え方を利用します。

具体的には、初学者の方も理解できるように以下の流れで説明していきます。

  1. マルコフ連鎖とは
  2. マルコフ連鎖が定常状態となる条件
  3. マルコフ連鎖の具体例
  4. マルコフ連鎖モンテカルロ法(MCMC)
  5. マルコフ連鎖モンテカルロ法(MCMC)の注意点

 

マルコフ連鎖とは

 

マルコフ連鎖とは、時刻\(t+1\)の状態が、\(t\)の状態のみに依存するモデルのことです。

より具体的には、以下が成り立つような確率変数の列\(\{\mathbf{X}^{t}\}\)のことをマルコフ連鎖と言います。

$$P(\mathbf{X}^{t+1} \mid \mathbf{X}^{t}, \ldots, \mathbf{X}^{1}) = P(\mathbf{X}^{t+1} \mid \mathbf{X}^{t})$$

一般に\(\mathbf{X}^{t}\)のとりうる値の集合\(\mathcal{S}\)を状態空間といい、\(P(\mathbf{X}^{t+1} \mid \mathbf{X}^{t})\)を遷移確率といい\(W(\mathbf{X}^{t+1} \mid \mathbf{X}^{t}) \equiv P(\mathbf{X}^{t+1} \mid \mathbf{X}^{t}) \)と表します(遷移確率が時間に依存しない場合)

この遷移確率を利用することで、確率分布の時間発展を以下のように記述することができます。

$$P_{t+1}(\mathbf{X}) = \sum_{\mathbf{X}^{\prime}} W(\mathbf{X} \mid \mathbf{X}^{\prime}) P_{t}(\mathbf{X})$$

ここで、遷移確率は以下の条件を一般に満たします。

遷移確率の条件

  1. \(W(\mathbf{X} \mid \mathbf{X}) \ge 0\)
  2. \(\sum_{\mathbf{X}} W(\mathbf{X} \mid \mathbf{X}^{\prime}) \) 

 

マルコフ連鎖が定常状態となる条件

 

以下の条件を満たすような状態を定常状態といいます。

定常状態

$$P_{\infty}(\mathbf{X}) = \sum_{\mathbf{X}^{\prime}} W(\mathbf{X} \mid \mathbf{X}^{\prime}) P_{\infty}(\mathbf{X})$$

 

このような定常状態に到達するためには遷移確率は以下の三つの条件を満たす必要があります。

定常状態に到達する遷移確率の条件

  • つりあい条件 : 以下のようにある状態\(\mathbf{X}\)に対して確率の流入と流出の合計が等しい
$$\sum_{\mathbf{X} \neq \mathbf{X}^{\prime}} W(\mathbf{X} \mid \mathbf{X}^{\prime}) P(\mathbf{X}^{\prime}) = \sum_{\mathbf{X}\neq \mathbf{X}^{\prime}} W(\mathbf{X}^{\prime} \mid \mathbf{X}) P(\mathbf{X}) $$
  • エルゴード条件 : 任意の二つの状態\(\mathbf{X}, \mathbf{X}^{\prime}\)が有限のステップで遷移可能
  • 周期解に落ち込まない

 

この三つの条件を満たせば、初期条件に依存しない唯一の定常分布が存在することが示せます。

 

マルコフ連鎖の具体例

 

状態空間\(\mathcal{S}\)を頂点で表し、遷移確率をエッジで表したものを状態遷移図といいます。

具体的な状態空間として(晴れ、くもり、雨)のケースを考え、遷移確率を以下の状態遷移図で与えられるような場合を考えてみましょう。

マルコフ連鎖の具体例、エッジ上の数字は遷移確率を表すマルコフ連鎖の具体例、エッジ上の数字は遷移確率を表します。

 

このとき、最初に確率1で晴れとなる状態からスタートしたマルコフ連鎖は、以下のようになります。

縦軸は状態のとりうる確率を表し、横軸は状態遷移の回数を表します。初期条件として晴れになる確率を1として行ったシミュレーション。縦軸は状態のとりうる確率を表し、横軸は状態遷移の回数を表します。

 

次に初期条件として、雨になりうる確率が1となる状態から始めた例を示します。

 

マルコフ連鎖の具体例 初期条件として雨の確率1初期条件として雨になる確率を1として行ったシミュレーション。縦軸は状態のとりうる確率を表し、横軸は状態遷移の回数を表します。

 

このように最終的に異なる初期値から始めても、同じ定常状態に収束することが確認できます。

下記にシミュレーションに使用した関数(Python)を共有します。

def mc(x0=[[1.0], [0.0], [0.0]], times=100):
    '''
    INPUT
    x0 : 初期条件
    times : 状態の更新回数
    OUTPUT
    state : マルコフ連鎖のサンプル列
    '''
    x = np.array(x0)
    state = np.zeros((times, 3))
    state[0] = x.ravel()
    for i in range(times-1):
        x = trans_prob@x
        state[i+1] = x.ravel()
    return state

 

適当な初期値から始めて皆さんも遊んでみてください。

 

マルコフ連鎖モンテカルロ法(MCMC)

 

サンプリングしたい確率分布\(P(\mathbf{X})\)が定常状態となるようなマルコフ連鎖を構築すれば、マルコフ連鎖を利用してサンプリングを行うことが可能です。

このようにマルコフ連鎖を利用してサンプル列を得る方法を『マルコフ連鎖モンテカルロ法(MCMC)』と言います。

具体的に、サンプルしたい確率分布が定常状態となるような遷移確率を構築するには、定常状態に到達する条件の一つであるつりあい条件をヒントに構築すれば良いです。

$$\sum_{\mathbf{X} \neq \mathbf{X}^{\prime}} W(\mathbf{X} \mid \mathbf{X}^{\prime}) P(\mathbf{X}^{\prime}) = \sum_{\mathbf{X}\neq \mathbf{X}^{\prime}} W(\mathbf{X}^{\prime} \mid \mathbf{X}) P(\mathbf{X}) $$

 

しかし、遷移確率を構築するためのヒントを得ることはできましたが、具体的にどのように遷移確率を構築するべきかはまだよくわかりません。

そこで、つりあい条件の十分条件である以下の詳細つりあい条件を満たすように遷移確率を構築する方法が一般に使用されます。

$$ W(\mathbf{X} \mid \mathbf{X}^{\prime}) P(\mathbf{X}^{\prime}) = W(\mathbf{X}^{\prime} \mid \mathbf{X}) P(\mathbf{X}) $$

具体的に詳細つりあい条件を満たすマルコフ連鎖モンテカルロ法は以下のような種類があります。

  1. メトロポリス法
  2. ギブスサンプリング(熱浴法)

この二つの手法を以下で簡単に説明します。

 

メトロポリス法

 

メトロポリス法では、遷移確率としていて以下のような確率分布を設定します。

$$W(\mathbf{X}^{\prime} \mid \mathbf{X}) = P_{p}(\mathbf{X}^{\prime} \mid \mathbf{X}) A(\mathbf{X}, \mathbf{X}^{\prime}) = P_{p}(\mathbf{X}^{\prime} \min \left(1, \frac{P(\mathbf{X}^{\prime}}{P(\mathbf{X})}  \right)$$

 

ここで、\(P_{p}(\mathbf{X}^{\prime} \mid \mathbf{X})\)を提案分布、\(A(\mathbf{X}, \mathbf{X}^{\prime})\)と呼びます。

この遷移確率を利用し、適当な初期条件\(\mathbf{X}^{0}\)を与えマルコフ連鎖を構築します。

この遷移確率は、詳細つりあい条件を満たすことを示せます。

 

ギブスサンプリング(熱浴法)

 

ギブスサンプリング(熱浴法)の場合、適当な初期状態\(\mathbf{X}^{0}\)を与え、現状態\(\mathbf{X}^{t} = (X_{1}, \ldots, X_{n})\)のうち一つの状態を選び、各状態\(X_{i}\)を以下の遷移確率で状態遷移します。

$$x_{i}^{t+1} \sim P(X_{i} \mid X^{t}_{1}, \ldots, X^{t}_{t-1}, X^{t+1}_{1}, \ldots, X^{t}_{n})$$

 

これを繰り返すことで、サンプルしたい確率分布を定常状態とするマルコフ連鎖を構築することができます。

 

マルコフ連鎖モンテカルロ法(MCMC)の注意点

 

マルコフ連鎖モンテカルロ法の注意点は以下の二つです。

  1. マルコフ連鎖の最初の数サンプルは定常状態からのサンプルではない
  2. 各サンプルは過去のサンプルに依存するためサンプル間に統計的な相関がある

①の対処法としては、マルコフ連鎖の最初の数ステップを期待値近似のためのサンプル列に含まなければ良いです。

一般に定常分布に到達するまでの時間をバーンインタイムと言います。

②の対処法としては、各サンプルが独立になるように一定間隔を空けて、期待値近似のためのサンプル列を取得すれば良いです。

または、異なる初期値から独立に複数のマルコフ連鎖を走らせて各マルコフ連鎖からサンプルを得るという方法も使用されます。

さらに、相関を弱めるための拡張された手法を利用することもあります。

参考資料

 

本記事と私が勉強するために使用した資料を説明していきます。

参考文献

 

参考文献を紹介します。

計算統計 II-マルコフ連鎖モンテカルロ法とその周辺 (統計科学のフロンティア)

 

 

理論から実装上の注意点がまとまっています。

 

ゼロからできるMCMC マルコフ連鎖モンテカルロ法の実践的入門 (KS理工学専門書) 

 

実際にプログラミングによる実装と説明がよくまとまっています。

 

やさしいMCMC入門: 有限マルコフ連鎖とアルゴリズム 

 

決してやさしくはなく理論もしっかり説明されています。

 

参考URL

 

参考にした資料に関するURLを下記にまとめます。

参考にしたYouTube講義は以下の講義です。

 

 

 

参考オンデマンド講義

 

Udemyの下記の講義もPythonによる実装と理論を理解するのに有益でした。

 

【PythonとStanで学ぶ】仕組みが分かるベイズ統計学入門

ベイズ統計の文脈で具体的にマルコフ連鎖モンテカルロ法を利用するので、実務で使用したい方必見の講義です。

 

まとめ

 

モンテカルロ法の基礎を簡単にまとめました。

皆さんも理解したら一度立ち止まり、式を追ってみたり、実際に実装してみたりして理解をより深く定着させてください。

モンテカルロ法は、確率統計理論の性質に基づき期待値評価の正当性が評価されます。

そのため、確率統計の性質を理解することで、モンテカルロ法の理解をさらに深めることができます。

もし、確率統計の勉強に興味がある方は下記の参考書を参考にしてください。

 

確実に力がつく統計学の参考書10選【東大生のオススメ】確実に力がつく統計学の参考書10選【東大生のオススメ】 この悩みと疑問を徹底解決していきます。 私が紹介する統計学の参考書を...

 

また、ベイズ統計をはじめ、機械学習の文脈でもマルコフ連鎖モンテカルロ法はよく使用されます。

例えば、データの背後にある確率分布を近似的に表現する生成モデル等でも頻繁に使用されます(下記参照)

 

【入門】生成モデルと統計的機械学習について(機械学習)機械学習の一つである生成モデルの基本的な枠組みについて説明しました。これらの枠組みを最初の段階で理解することで高度なモデルも理解することが可能になります。...

 

 

ABOUT ME
努力のガリレオ
【運営者】 : 東大で理論物理を研究中(経歴)東京大学, TOEIC950点, NASA留学, カナダ滞在経験有り, 最優秀塾講師賞, オンライン英会話講師試験合格, ブログと独自コンテンツで収益6桁達成 【編集者】: イングリッシュアドバイザーとして勤務中(経歴)中学校教諭一種免許取得[英語],カナダ留学経験あり, TOEIC650点