本稿ではMCMC法の解説のため、MC法による積分の計算方法（モンテカルロ積分）から、MCMCによる手法の概要を見ていく。MCMC法は有名かつ知り尽くされた手法で、多くの良質な説明資料 [1], [2], [3], [4], [5] が存在している。従ってここの説明は読まずに、資料を見てもらった方が理解が早いかもしれない。
一般に MC（Monte-Calro, モンテカルロ）法 は、サンプリング（サンプルを乱数から生成すること）によってシミュレーションや数値計算を行う手法である。特に確率分布が関わる積分値 [6] を近似的に求めるMC法はモンテカルロ積分と呼ばれる。モンテカルロ積分は確率的な推論の一種であり、大数の法則 [7] によって、十分なサンプル数をとれば近似精度をいくらでも良くする事ができる。サンプリングの手間がある為、近似分布をあらかじめ仮定する様な決定論的な推論よりも遥かに推論が遅い。しかし、MCは近似分布が求められないような場合にも適用可能であり、汎用性が高いと言える。
MC法によって原理的には任意の解を求められるが、十分なサンプル数の要求というのが大きな問題を孕んでいる。サンプリングの自由度（範囲及び次元）が大きくなると、解の計算にあまり寄与しない（無駄な）サンプルが増えてしまう。計算を現実的かつ効率的に行うためには、サンプルの選択が重要になる。
そして MCMC（Markov Chain Monte-Calro, マルコフ連鎖モンテカルロ）法 は、新しいサンプルを以前に生成したサンプルに確率的に依存して（サンプルの列がマルコフ連鎖となる様に）生成するMC法である。MCMCでは、新しく生成したサンプルを採択（採用）するか棄却（捨てる）するかも確率的に判断する。この手続きによって、無駄なサンプルを極力減らすようにサンプリングを実行することができる。
MC法による積分 - モンテカルロ積分
確率変数を\(d\)次元の実数値ベクトル [8] \(\boldsymbol{x} = [x_{1},\dots,x_{d}]^{\mathsf{T}} \in X \subset \mathbb{R}^{d}\)とする。ここで\(X\)は全事象 [9] の集合である。\(\boldsymbol{x}\)の確率分布を\(r(\boldsymbol{x})\)とし、関数\(h\)の確率分布\(r\)による平均（期待値）
を求めることを考える。ここで、\(\mathrm{E}_{p}[\cdot]\)は確率分布\(p\)による平均を表す。\(I\)において、関数\(h\)の形に制約を与えておらず積分として様々な値が計算できる。例を挙げると:
- \(h(\boldsymbol{x}) = \boldsymbol{x}\) : この場合は\(\mathrm{E}_{r}[\boldsymbol{x}]\)、即ち\(\boldsymbol{x}\)の平均を求める
- \(h(\boldsymbol{x}) = (\boldsymbol{x} - \mathrm{E_{r}}[\boldsymbol{x}])(\boldsymbol{x} - \mathrm{E_{r}}[\boldsymbol{x}])^{\mathsf{T}}\) : \(\boldsymbol{x}\)の分散を求める
- … その他 [10]
もし\(r(\boldsymbol{x})\)が既知で、分布\(r\)から簡単に独立にサンプリングできる [11] ならば、\(r(\boldsymbol{x})\)からの独立な（他のサンプルに依存して生成しない）\(n\)個のサンプルを\(\boldsymbol{x}_{1}, \boldsymbol{x}_{2}, \dots, \boldsymbol{x}_{n}\)と書くと、\(I\)の標本平均による近似値\(\hat{I}\)は
で計算できる。大数の法則により、サンプル数の極限を取れば標本平均は真の平均に一致する:
この様にして平均を求める方法を モンテカルロ積分(Monte-Carlo Integration) という。一般にモンテカルロ法(Monte-Carlo Method)はサンプリングによってシミュレーションや数値計算を行う事を指す。
重点サンプリング
モンテカルロ積分によって、原理的には\(\hat{I}\)を多くのサンプルで計算する事で\(I\)を精度良く計算できる。しかし実際確率分布\(r(\boldsymbol{x})\)は複雑であることが多く、その場合\(r(\boldsymbol{x})\)から直接サンプリングするのは困難となる。そこで、より簡単でサンプリング可能な確率分布（提案分布 という）\(q(\boldsymbol{x})\)を用意して、そこからサンプリングする事を考える。\(q(\boldsymbol{x})\)を使えば、\(I\)は次の様に変形できる:
モンテカルロ積分の時と同じ様にに考え、次は\(\boldsymbol{x}_{1},\dots,\boldsymbol{x}_{n}\)を\(q(\boldsymbol{x})\)からの独立な\(n\)個のサンプルにすれば、\(I\)の近似値\(\hat{I}_{IS}\)として
が得られる。ここで\(w(\boldsymbol{x}_{i}) = r(\boldsymbol{x}_{i})/q(\boldsymbol{x}_{i})\)はサンプル\(\boldsymbol{x}_{i}\)に対する重みと見ることができる。この様に、重みが付いたサンプルで平均を求める手法を 重点サンプリング(Importance Sampling) という。重点サンプリングにおいても、\(q(\boldsymbol{x})\)がある条件を満たしていれば、大数の法則によって\(\displaystyle\lim_{n \to \infty} \hat{I}_{IS} = I\)となることが保証されている。
MCMC
重点サンプリングの考え方によって、確率分布\(r\)が複雑でも替わりに提案分布\(q\)を用いてサンプリングを行えばモンテカルロ積分が計算できる事が確かめられた。しかし、" \(r\)より簡単でサンプリング可能な\(q\)" を構成する事自体が一般に困難である。特に次元\(d\)が増加すれば\(r\)が複雑になるのはもちろん、全事象\(X\)の自由度が増加し次元の呪い [12] を引き起こす。即ち、\(r\)を\(q\)で良く近似出来てない時に毎回独立にサンプリングを行っていると、空間\(X\)から当てずっぽうなサンプルを取得しているのと同様な状態になる。
そこで、簡単な提案分布\(q\)を用いて、かつ逐次的に以前のサンプルを使用して新しくサンプルを生成する手法が90年代以降使われる様になってきた。この場合、サンプル列はマルコフ連鎖(Markov Chain)をなす。そして、マルコフ連鎖で生成したサンプルによるMC法をMCMC（Markov Chain Monte-Calro）法という。サンプル間の独立性は担保されなくなる為にMC法の基本原理が成立しなくなるが、提案分布（マルコフ連鎖の遷移確率）がある性質を満たせば、十分なサンプル数で確率分布\(r\)からのサンプリングが実現できる。
遷移確率の条件 - 詳細釣り合い条件
概要でも既に述べたが、MCMCは生成したサンプル列がマルコフ連鎖をなすように生成する。今、サンプル列\(\boldsymbol{x}_{0}, \boldsymbol{x}_{1}, \dots\)はマルコフ連鎖をなすので、生成した時刻（ステップ）で実際に観測した状態を\(\boldsymbol{e}_{0}, \boldsymbol{e}_{1}, \dots \ (\boldsymbol{e}_{i} \in X \ i=0,1,\dots)\)と書くと、任意の時刻\(n \geq 0\)で、
が成り立つ（この性質をマルコフ性 [13] という）。即ち、サンプルは直前のサンプルのみに依存して生成する。この様にサンプルを生成する場合、実はマルコフ連鎖が エルゴード的(ergodic) という性質を満たせば、大量のサンプルを用いた時にある分布（定常分布）\(\pi\)からサンプリングしているのと同様になる。
マルコフ連鎖がエルゴード的であるとは、規約性（どの状態からでも任意の状態へ遷移できる）と正再帰性（任意の状態へ何回でも遷移できる）非周期性（任意の状態は一回の遷移で元に戻れる）を全て同時に満たすことを言う [14]。 エルゴード的なマルコフ連鎖と定常分布\(\pi\)の関係は、次の定理で表せる:
マルコフ連鎖の収束
マルコフ連鎖\(\boldsymbol{x}_{0}, \boldsymbol{x}_{1}, \dots\)がエルゴード的であり、その遷移確率行列を\(\boldsymbol{P}\)とおく。\(\pi\)を\(\boldsymbol{P}\)の定常（不変）分布とした時、任意の初期状態から始まるマルコフ連鎖はサンプル数の極限において定常分布\(\pi\)に収束する。
ここで遷移確率行列\(\boldsymbol{P}\)とは、その\((i,j)\)要素\((\boldsymbol{P})\_{ij} = p_{ij}\ (i,j \in X)\)が任意の時刻\(t \geq 0\)で
を満たすような行列である [15]。 また、定常分布とは時刻が経過しようとも不変なマルコフ連鎖（一般に確率過程）の各状態の確率分布である [16]。即ち、十分に長いマルコフ連鎖を観測すれば、どの状態にいる傾向があるのかを定常分布によって知ることができる。
上記の議論により、マルコフ連鎖がエルゴード的であればサンプリングが定常分布に従う事は分かったが、次は遷移確率の設計が問題となる。遷移確率を規約性と正再帰性と非周期性とを満たすように設定するのは案外容易 [17] であるが、それだけでは定常分布の存在のみを保証するので、その定常分布が希望する分布に一致するとは限らない。次に問題となるのは、希望の確率分布\(r\)を定常分布とするように遷移確率を設計することである。その問題は次の 詳細釣り合い条件(detailed balance condition) という条件によって解決できる。
詳細釣り合い条件
希望する確率分布\(r\)と遷移確率\(p\)が次の条件を満たす時、そのマルコフ連鎖の定常分布\(\pi\)は\(r\)に一致する:
ここで\(r_{i} = r(\boldsymbol{x} = i)\)である（証明は補足に示した）。
詳細釣り合い条件を満たす遷移確率を用いさえすれば、十分大きな\(m>0\)を取った時に、マルコフ連鎖\(\boldsymbol{x}_{m}, \boldsymbol{x_{m+1}},\dots\)は\(r\)からのサンプルとなる。 次の節で紹介するアルゴリズムの遷移確率は、いずれも詳細釣り合い条件を満たすように設計されている。
メトロポリス-ヘイスティングス法
メトロポリス-ヘイスティングス法は、サンプルは重点サンプリングの時と同じように提案分布によって生成し、そして新しく生成したサンプルを 採択（採用）するか、もしくは 棄却 （捨てる）のかを 採択確率(acceptance rate) と呼ばれる確率によって決め、採択された場合は新しい状態に遷移し、棄却された場合には遷移は行わずに（状態を変えずに）もう一度サンプリングし直す、という手続きを繰り返す手法である。
メトロポリス-ヘイスティングス法の更新規則を導出してみる。 まず、状態\(i \in X\)から状態\(j \in X\)に遷移する時の提案分布を条件付き確率\(q(\boldsymbol{x}_{n+1}=j|\boldsymbol{x}_{n}=i) = q_{ij}\)と書き、また状態\(i\)にいる時に状態\(j\)を採択する確率（採択確率）を\(\alpha(i \to j)\)と表す。すると、\(i\)から\(j\)への遷移確率\(p_{ij}\)は\(q_{ij}\)と\(\alpha(i \to j)\)の積で表せる:
そして、詳細釣り合い条件から、
となる。採択確率はこの条件を満たす様に設計する。メトロポリス-ヘイスティングス法では特に、
とする [18]。アルゴリズムの実行中には、この式によって採択確率を計算し、\([0,1]\)の範囲の一様乱数を発生させて採択/棄却を判断する。
これでメトロポリス-ヘイスティングス法が実行できるが、その利点を2つ挙げる:
\(r\)が厳密計算出来なくても良い \(r\)は一般に複雑なので直接的な計算は難しいが、上の採択確率の式は確率の比率のみに注目している。従って分布が厳密に計算できなくてもアルゴリズムを実行できる。比率さえ一致すれば良いので、分布\(r\)の近似分布\(\hat{r}\)として
\begin{equation*} \hat{r} = \frac{1}{Z_{r}} r \end{equation*}としても良い事になる(\(Z_{r}\):正規化定数)。特に、近似分布をボルツマン-ギブス分布
\begin{equation*} \hat{r}(\boldsymbol{x}) = \frac{1}{Z_{r}} \exp(-r(\boldsymbol{x})/T) \end{equation*}とする場合が多い。ここで、\(T>0\)は温度パラメタ [19] である。
\(q_{ij} = q_{ji}\)が成り立つ場合には、より簡単にサンプリングできる \(q_{ij} = q_{ji}\)が成立する提案分布で有名なものに酔歩連鎖(random walk chain)がある:
\begin{equation*} q_{ij} = {\cal N}(i, \sigma^{2}\boldsymbol{I}) \end{equation*}即ち平均（中心）を現在状態\(i\)、分散を\(\sigma\) [20] とした正規分布からの乱択でサンプリングを行う [21] 。 正規分布以外でも、\(i\)を平均とした一様分布、多変量\(t\)分布でも実行できる。
ギブスサンプリング
ギブスサンプリング(Gibbs Sampling, 熱浴法とも)は提案分布の変数を1個ずつ更新していく手法である。 主に多次元確率分布 [22] の推定に用いられる事が多い。説明のため、現在の状態を組\(\boldsymbol{x} = (x_{1}, x_{2}, \dots, x_{d})\)と書く。状態の更新の際には、変数を1つ選び出し [23] て\(x_{i} \to x_{i}^{\prime}\)と遷移させる(\(i=1,\dots,d\))。更新後の状態を\(\boldsymbol{x}^{\prime} = (x_{1}, \dots, x_{i-1}, x_{i}^{\prime}, x_{i+1}, \dots, x_{d})\)と書く。ここで、遷移確率\(q(\boldsymbol{x}^{\prime}|\boldsymbol{x})\)は次で定義される:
即ち、選択した変数\(x_{i}\)以外を全て``固定’’した確率分布\(r\)から\(x_{i}^{\prime}\)を新しくサンプリングする。上記右辺が計算できる場合にのみ、ギブスサンプリングは適用可能となる。
この更新規則が詳細釣り合い条件を満たすことは、再びベイズの定理を用いて、
により確認できる。また、メトロポリス-ヘイスティングス法の採択確率の式から、
となり、ギブスサンプリングはメトロポリス-ヘイスティングス法で採択確率を\(1\)（必ず採択）するようにした特別の場合である事が分かる。採択/棄却の手順を踏まくくても良く、しかも遷移確率\(q\)は予め計算できるので、高速な推定ができるようになっている。
MCMCによる最適化
MCMCは関数最適化に用いることもできる。今、サンプリングを行う確率分布をボルツマン-ギブス分布
とした時、定義式により、\(f(\boldsymbol{x})\)が小さな値を与える点ではその確率\(r(\boldsymbol{x})\)は同時に大きくことが即座に観察できる。従って、MCMCによって\(r(\boldsymbol{x})\)からのサンプリングを行えば、\(f(\boldsymbol{x})\)が小さな値をとる点を集中してサンプリングできる事から、\(f(\boldsymbol{x})\)の最小化（最大化の場合は\(-f(\boldsymbol{x})\)の最小化に置き換えれば良い）を考える事ができる。実際、関数\(f\)の最小値を与える点を\(\boldsymbol{x}^{\ast}\)と表せば、サンプル数\(N\)の極限において最小値\(f(\boldsymbol{x}^{\ast})\)が確率1で得られる事:
が示せる。以下、その証明を示す。
（証明） MCMCにおいて、定常分布を\(r\)とする様に（詳細釣り合い条件を満たす様に）サンプリングを行う。この時マルコフ連鎖\(\boldsymbol{x}_{1}, \boldsymbol{x}_{2}, \dots, \boldsymbol{x}_{n},\dots\)は、十分大きな\(n > 1\)においては\(r(\boldsymbol{x})\)からのサンプルとみなせる。関数\(f\)に最小値\(f(\boldsymbol{x}^{\ast})\)が存在すれば、\(\boldsymbol{x}^{\ast}\)をサンプリングする確率\(r(\boldsymbol{x}^{\ast})\)も存在が保証され、分布の中で最大の確率を与えている。従って、\(n\)回目以降のマルコフ連鎖\(\boldsymbol{x}_{n}, \boldsymbol{x}_{n+1},\dots\)において、\(m \geq n\)回目に初めて\(\boldsymbol{x}^{\ast}\)がサンプリングできる確率\(P(\boldsymbol{x}_{m} = \boldsymbol{x}^{\ast})\)は、幾何分布と同じ様に、
によって計算できる。また、初めて\(\boldsymbol{x^{\ast}}\)がサンプリングできるまでの回数が\(N \geq n\)回以内となる確率は、
となる。 ここでサンプル数の極限\(N \to \infty\)をとると、初項\(r(\boldsymbol{x}^{\ast})\)、項比\(1-r(\boldsymbol{x}^{\ast})\)とした等比級数の和の公式より、
が得られる。即ち、サンプリングを無限に繰り返せば\(\boldsymbol{x}^{\ast}\)が確率1で得られることが示された。この結果は、サンプルの関数列\(f(\boldsymbol{x}_{1}), f(\boldsymbol{x}_{2}), \dots\)の中に少なくとも1つ\(f(\boldsymbol{x}^{\ast})\)が存在する事と同値である。
焼きなまし法
以上でMCMCによる最適化が理論的に可能なことが示されたが、最適化の際に特に問題となるのは分布\(r\)の温度パラメタ\(T\)である。\(T\)が大きければ、\(\exp\)内部の\(f(\boldsymbol{x})\)の値に影響されず\(r(\boldsymbol{x})\)は一様分布に近くなり、一様乱数からのサンプリングと殆ど変わらなくなる。逆に\(T\)が\(0\)に近いと\(r(\boldsymbol{x})\)は\(f(\boldsymbol{x})\)の値に大きく影響されるが、サンプリングが特定の場所だけに集中してしまって局所最適値しか得られない場合がある。この様に\(T\)は適切に決定する必要があるが、\(T\)の適切な決定法は存在せず、問題依存となる場合が多い。
そこで、最初は\(T\)（温度）を高い状態から初めてサンプリングの度に少しずつ\(T\)を下げていくやり方があり、これを焼きなまし法（Simulated Annealing, SA）と呼ぶ。この様に\(T\)を変化させると最初は空間全体の中から大雑把な\(f\)の値を取得し、後に最適値の近傍を集中してサンプリングすることができるために効率的な探索が期待できる。証明は省くが、温度パラメタの系列\(T_{1}, T_{2}, \dots\)が次の条件を満たせばサンプリングによって\(\boldsymbol{x}^{\ast}\)が得られる事（収束定理）が示されている:
ここで、\(D\)は問題によって決まる定数である。
エルゴード的なマルコフ連鎖の定常分布
上記の議論で、「マルコフ連鎖がエルゴード的ならば、一意な定常分布が存在する」という事に触れた。この定理についての証明を述べていくが、準備として確率過程についての用語や記法の定義、基本的な定理の証明を行う。大方の証明はここを参照した。なお、状態空間（全事象）\(X\)は有限集合であるとする。
離散時間マルコフ連鎖
確率過程（サンプル列） \(\boldsymbol{x}_{0}, \boldsymbol{x}_{1}, \dots\) が次を満たす時、離散時間マルコフ連鎖という。
またこの性質をマルコフ性という。
遷移確率の斉時性、nステップ遷移確率
任意の状態\(i,j \in X\)と非負整数\(n \geq 0\)に対して
を、状態\(i\)から状態\(j\)への遷移確率という。\(p_{ij}(n)\)が\(n\)と独立で常に\(p_{ij}(n) = p_{ij}(0) = p_{ij}\)となる時、離散時間マルコフ連鎖は斉時であるという。今後、遷移確率は\(p_{ij}\)を用いて表す。また、
は状態\(i\)から始まって\(n\)ステップ後に状態が\(j\)になる確率を表しており、\(n\)ステップ遷移確率と呼ぶ。
チャップマン−コルモゴロフ方程式
任意の状態\(i,j \in X\)に対し、\(n\)ステップ遷移確率\(p_{ij}^{(n)}\)は次を満たす:
（証明）
到達可能、連結
ある状態\(i,j \in X\)に対して\(p_{ij}^{(n)} > 0\)なる非負整数\(n \geq 0\)が存在する時、状態\(j\)は状態\(i\)から到達可能であると言い、\(i\to j\)と表す。 また\(i \to j \land j \to i\)ならば、\(i\)と\(j\)は連結しているといい、\(i \leftrightarrow j\)と表す。
連結関係は、反射性(\(i \leftrightarrow i\))、対称性(\(i \leftrightarrow j \Leftrightarrow j \leftrightarrow i\))、推移性(\(i \leftrightarrow j \land j \leftrightarrow k \Rightarrow i \leftrightarrow k\))が成り立つ。
連結クラス（連結成分）
\(X\)の部分集合\(C \subseteq X\)において、
- \(i \in C \land j \in C \implies i \leftrightarrow j\)
- \(i \in C \land i \leftrightarrow j \implies j \in C\)
が常に成立する時、\(C\)を\(X\)の連結クラス（連結成分）という。定義より、\(C\)の要素は互いに連結している。また、連結クラス\(C\)の任意の状態\(i \in C\)から\(j \notin C\)に到達できない時、\(C\)は閉じていると言う。
規約性
\(X\)内の全ての状態が単一の閉じた連結クラスに属する、即ち\(X\)の全ての要素が互いに連結している時、そのマルコフ連鎖は規約であるという。
周期性
状態\(i \in X\)に対して\(p_{ii}^{(n)} > 0\)となる（\(n\)ステップ後に元の状態に戻る）\(n\)の最大公約数\(d\)を、状態\(i\)の周期と呼ぶ。\(d = 1\)の時は状態\(i\)は非周期的と呼ばれ、\(d \geq 2\)の時は周期的であると呼ばれる。
再帰的、過渡的
確率変数\(T_{j}\)を次で定義する:
即ち、離散時間マルコフ連鎖が初めて状態\(j\)を訪れる時刻を表す。また、\(T_{i}\)を用いて次の値を定義する:
\(f_{i}\)は将来状態\(i\)に戻ってくる確率を表しており、\(f_{i}=1\)ならば確率\(1\)で状態\(i\)を訪れる（無限にしばしば訪れる）ので状態\(i\)は再帰的であるという。\(f_{i} < 1\)ならば状態\(i\)は過渡的であるという。また、\(m_{i}\)は初期状態が\(i\)の時に、再び状態\(i\)に戻るまでの時間の期待値を表しており、\(m_{i} < \infty\)ならば状態\(i\)は正再帰的（有限時間で\(i\)に戻る）であるといい、\(m_{i} = \infty\)ならば状態\(i\)は零再帰的であるという。
エルゴード的な離散時間マルコフ連鎖
離散時間マルコフ連鎖\(\boldsymbol{x}_{0}, \boldsymbol{x}_{1},\dots\)が規約かつ正再帰かつ非周期的であるならば、この離散時間マルコフ連鎖はエルゴード的とも呼ばれる
ここまでで用語の定義は揃ったので、それではエルゴード的なマルコフ連鎖の定常分布の存在についての定理を証明する。
エルゴード的な離散時間マルコフ連鎖の定常分布
離散時間マルコフ連鎖\(\boldsymbol{x}_{0}, \boldsymbol{x}_{1}, \dots\)がエルゴード的ならば、任意の状態\(i, j \in X\)について次が成り立つ:
- \(\displaystyle\lim_{n \to \infty} p_{ij}^{(n)} = \lim_{n \to \infty} p_{jj}^{(n)} = \frac{1}{m_{j}} = \pi_{j}\)
- \(\pi_{j}\)は\(\displaystyle \pi_{j} = \sum_{i \in X} \pi_{i} p_{ij}\)と\(\displaystyle\sum_{j \in X}\pi_{j} = 1\)を満たす解であり、唯一に定まる。
2.を満たす\(\pi_{j}\)を極限分布（定常状態分布）と言う。
一方、初期分布として\(P(\boldsymbol{x}_{0} = j) = \pi_{j}\)を持つ離散時間マルコフ連鎖では、任意の\(n \geq 1\)に対して\(P(\boldsymbol{x}_{n}=j) = \pi_{j}\)が成り立ち、\(\boldsymbol{x}_{n}\)は\(n\)と独立した分布を持つ。この様に、時間に関して不変な分布\(\pi_{j} = P(\boldsymbol{x}_{n} = j)\ n = 0,1,\dots\)を 定常分布 と呼ぶ。
（証明）まず1.から考える。最初に\(i\neq j\)なる状態に対して
を（初期状態が\(i\)で、初めて\(j\)に訪れる時刻が\(k\)となる確率）おく。この時、
の観察により、\(n \geq 1\)なる\(n\)に対して帰納的に
が成立する（最初の\(k\)ステップで状態\(j\)に行き、その後\(n-k\)ステップ後に再び\(j\)に行く）ことが分かる。また、任意の\(i\)と\(j\)は連結している（\(i \leftrightarrow j\)）ので、
（状態\(i\)から始まり、\(j\)へいつかは訪れる確率は\(1\)）が成り立つ。一方\(p_{jj}^{(n)}\)は、
と展開できる。数列\(p_{jj}^{(n)}\)の極限\(\displaystyle\lim_{n \to \infty} p_{jj}^{(n)}\)を求める為、ここでは数列の 母関数を定義し、（片側）Z変換の最終値定理 [24] を用いる。その為、今、\(\displaystyle G(z) = \sum_{n=0}^{\infty}p_{jj}^{(n)}z^{n},\ U(z) = \sum_{n=1}^{\infty}u_{n}z^{n}\)なる母関数を定義し、上式の両辺に\(z^{n}\)を掛けて\(n=1,2,\dots\)についての和を取ると、
ここで、右辺式の最後の式変形には冪級数の積の公式 [25] を用いている。最終値定理を適用する事を考えると、この場合は、
が成立する [26] ので、\(\displaystyle\lim_{n \to \infty} p_{jj}^{(n)}\)の結果として、
が得られる。さて、この結果より、任意の正数\(\epsilon > 0\)に対して\(n \geq N\)なる全ての\(n\)が
を同時に満たすような\(N\)を取ることができる。今、\(n \geq 2N\)に対し、
よって、\(\displaystyle \lim_{n \to \infty} p_{ij}^{(n)} = \pi_{j} = \lim_{n \to \infty} p_{jj}^{(n)}\)。
次に2. の\(\pi_{j}\)の一意性を示す。まず、\(\displaystyle \sum_{j \in X}p_{ij}^{(n)} = 1\)（どこかの状態には確率1で遷移している）より、この式で\(n \to \infty\)ならしめれば、1. により
を得る。また、\(a_{j}(n) = P(\boldsymbol{x}_{n} = j)\)（時刻\(n\)で状態\(j\)を訪れる確率）とおくと、
が成立し、チャップマン−コルモゴロフ方程式により、\(n, m \geq0\)なる整数に対し、
この式の両辺を\(m \to \infty\)ならしめれば、極限と和の交換法則より、
を得る。特に\(n=1\)とすれば、\(\displaystyle \pi_{j} = \sum_{i \in X} \pi_{j} p_{ij}\)が得られる。 次に一意性を示す。今、\(\pi_{i}^{\prime}\ (i \in X)\)が、
を満たすとする。上述の議論により、全ての正整数\(n \geq 0\)に対し、
を得る。\(n \to \infty\)とすると、
となって、一意性が示される。
詳細釣り合い条件の証明
最後に詳細釣り合い条件を示す。今、確率分布\(r\)と遷移確率が
を満たしているとする。この時両辺ともに状態\(i\)について和をとると、
2.により、\(r\)は定常分布の解となっている事が分かる。
脚注
| [1] | 古澄英雄, 「21世紀の統計科学」第Ⅲ巻 第10章 マルコフ連鎖モンテカルロ法入門 |
| [2] | 笠原正治, 確率過程論基礎 |
| [3] | 中川裕志, マルコフ連鎖モンテカルロ法 |
| [4] | 福島孝治, マルコフ連鎖モンテカルロ法の実践 |
| [5] | tera monagi, マルコフ連鎖モンテカルロ法入門-1 |
| [6] | 主に、確率分布の平均（期待値）、分散が対象となる |
| [7] | 十分な回数の独立な試行を行った経験分布は理論的（真の）分布に一致する、という法則。例えばコイン投げをひたすら繰り返せば、表及び裏が出る 頻度の比率 はそれぞれ\(1/2\)に近づいていく。厳密には大数の法則は2種類（強、弱法則）あり、確率の応用において非常に非常に重要な法則であるが、ここでは説明をしない。 |
| [8] | 確率変数のとる値が実数値でなくとも、事象が有限個存在（\(\iff\)全事象が有限集合）する場合（例。サイコロとかコインを投げる試行）は議論で用いている分布を離散確率分布で考えれば良い。 |
| [9] | 起こりえる全ての事象の集合。 |
| [10] | 他の個人的に興味深い例：強化学習において\(X\)を選択した行動列の集合、\(h:X \to \mathbb{R}\)を報酬関数とすれば、\(h(\boldsymbol{x})\)で行動列の報酬が計算でき、\(I\)の計算結果は報酬の期待値となる。報酬の期待値が計算できることはエージェントの行動決定において大変有用である。 |
| [11] | 一様分布や正規分布等のよく知られた分布は、サンプリングアルゴリズムも確立されている。一様分布はメルセンヌ・ツイスタ、正規分布にはボックス-ミューラー法といった具合である。 |
| [12] | 空間の次元が増加すると、その空間の自由度が直感に反して 指数的 に増加すること。例えば、ユークリッド空間で一辺の長さが\(a\)の\(n\)次元超立方体を占める直径\(a\)の超球体の割合を計算してみると\(\frac{\sqrt{(\pi(a/2)^{2})^{n}}}{a^{n} \Gamma(\frac{n}{2}+1)}\)であり、\(n\)を増加させると階乗オーダー（即ち、指数オーダーよりも早く）で減少する事が分かる。従って、一様乱数を用いていると、\(n\)次元空間で超球体の内部にサンプルが入る確率が階乗オーダーで小さくなる。 |
| [13] | 厳密には、直前の1つのサンプルのみに依存するので1階マルコフ性と呼ばれる。 |
| [14] | 詳細は補足で述べる。一般にエルゴード的とは、長時間に渡って観測した状態の平均（長時間平均）と、状態空間の平均（位相平均）が一致するという事を表す概念である。エルゴード理論がある様に、厳密な数学理論が展開されるが、ここではマルコフ連鎖以外については詳しくは説明しない（筆者がついていけてない）。 |
| [15] | 連続な状態空間では、遷移確率行列の代わりに
\begin{equation*}
P(\boldsymbol{x}_{t+1} \in C|\boldsymbol{x}_{t} = \boldsymbol{e}_{t}) = \int_{C} T(\boldsymbol{e}_{t}, \boldsymbol{y}) d \boldsymbol{y} \quad C \subset X, \boldsymbol{e}_{t} \in X
\end{equation*}
となる様な条件付き確率分布\(T(\boldsymbol{x}, \boldsymbol{y})\)（遷移核）を用いれば良い。 |
| [16] | 形式的に書くと、状態\(j \in X\)の定常分布\(\pi_{j}\)は\(\pi_{j} = P(\boldsymbol{x}_{n} = j)\ n=0,1,\dots\)と表される。 |
| [17] | 例えば、現在状態を中心とした正規分布からでの乱択でも3つの性質を満たし、マルコフ連鎖はエルゴード的となる。 |
| [18] | これが詳細釣り合い条件を満たすことは、場合分けにより分かる:
\begin{equation*}
p_{ji} = q_{ji} \alpha(j \to i) = q_{ji} \frac{r_{i}q_{ij}}{r_{j}q_{ji}} = \frac{r_{i}}{r_{j}}q_{ij} = \frac{r_{i}}{r_{j}} p_{ij} \iff r_{i}p_{ij} = r_{j}p_{ji}
\end{equation*}
\begin{equation*}
p_{ij} = q_{ij} \alpha(i \to j) = q_{ij} \frac{r_{j}q_{ji}}{r_{i}q_{ij}} = \frac{r_{j}}{r_{i}}q_{ji} = \frac{r_{j}}{r_{i}} p_{ji} \iff r_{i}p_{ij} = r_{j}p_{ji}
\end{equation*}
|
| [19] | 温度パラメタの調節は非常に難しい事が知られている。実験結果を見て経験的に設定される事がほとんどである。 |
| [20] | 分散パラメタの調節も非常に難しい。分散を大きくすると遷移幅（ステップサイズという）が大きくなって定常分布に落ち着くまでに時間が掛かり、分散を小さくし過ぎると遷移の動きが小さく、探索が十分に行われない危険性がある。一般に分散パラメタと温度パラメタにはトレードオフの関係がある。 |
| [21] | \(q_{ij} = q_{ji}\)が成立する理由は、この場合\(j\)は
\begin{equation*}
j = i + \varepsilon \quad \varepsilon \sim {\cal N}(\boldsymbol{0}, \sigma^{2} \boldsymbol{I})
\end{equation*}
と書けるので、平均が\(\boldsymbol{0}\)かつ正規分布の対称性により、
\begin{equation*}
i = j - \varepsilon = j + \varepsilon
\end{equation*}
よって\(q_{ij} = {\cal N}(i, \sigma^{2}\boldsymbol{I}) = {\cal N}(j, \sigma^{2}\boldsymbol{I}) = q_{ji}\)を満たす |
| [22] | 機械学習では、ベイジアンネットワークやボルツマンマシン（深層学習の一部）等のモデル学習に使われる |
| [23] | 毎回ランダムで選んでも、順番に全変数を1個ずつ選んでも良い |
| [24] | 数列\(a_{n}\)の母関数を\(F(z) = \displaystyle\sum_{n=0}^{\infty}a_{n}z^{n}\)とする。今、複素数\(s \in \mathbb{C}\)を用いて\(z = \exp(-s)\)とおき、\(n\)の和を\(t\)の積分に置き換えると、
\begin{equation*}
F(\exp(-s)) = \int_{0}^{\infty} a_{t}\exp(-st) dt
\end{equation*}
これは数列\(a_{t}\)のラプラス変換に他ならない。従ってラプラス変換の最終値定理を適用できる。離散の場合のラプラス変換を（片側）Z変換と呼ぶ。 |
| [25] | 2つの冪級数を\(\displaystyle\sum_{n=0}^{\infty}a_{n}z^{n}, \sum_{n=0}^{\infty}b_{n}z^{n}\)とし、積の結果を\(\displaystyle\sum_{n=0}^{\infty}c_{n}z^{n}\)とする。等号を立てると、
\begin{align*}
\left(\sum_{n=0}^{\infty}a_{n}z^{n}\right)\left(\sum_{n=0}^{\infty}b_{n}z^{n} \right) &= a_{0}b_{0}z^{0} + (a_{0}b_{1} + a_{1}b_{0})z^{1} + (a_{0}b_{2}+a_{1}b_{1}+a_{2}b_{0})z^{2} + \dots \\
&= \sum_{n=0}^{\infty}c_{n}z^{n} = c_{0}z^{0} + c_{1}z^{1} + c_{2}z^{2} + \dots
\end{align*}
係数比較により、\(c_{0} = a_{0}b_{0},\ c_{1} = a_{0}b_{1} + a_{1}b_{0},\ \dots\)が成立し、よって\(c_{n} = \displaystyle \sum_{k=0}^{n}a_{k}b_{n-k}\)となる。ここの例では、\(a_{k} = u_{k}z^{k},\ b_{k} = p_{jj}^{(k)}z^{k}\)とおけば良い。 |
| [26] | （証明）母関数（Z変換）を\(F(z) = \displaystyle \sum_{n=0}^{\infty}a_{n}z^{n}\)とおくと、
\begin{align*}
\lim_{z \to 1} (1-z) F(z) &= \lim_{z \to 1}(1-z) \sum_{n=0}^{\infty} a_{n}z^{n} = \lim_{z \to 1} \sum_{n=0}^{\infty} a_{n} (z^{n} - z^{n+1}) = \lim_{z \to 1} \lim_{n \to \infty} \sum_{k=0}^{n} a_{k} (z^{k} - z^{k+1}) \\
&= \lim_{z \to 1} \lim_{n \to \infty} \sum_{k=0}^{n} (a_{k} - a_{k-1}) z^{k} \quad (\because a_{-1} = 0, また a_{0}(z^{0}-z^{1}) + a_{1}(z^{1}-z^{2}) +\dots = a_{0}z^{0} + (a_{0}-a_{1})z^{1} + \dots) \\
&= \lim_{n \to \infty} \sum_{k=0}^{n} (a_{k} - a_{k-1}) = \lim_{n \to \infty} a_{n}
\end{align*}
|