この記事を読んでできること
MCMC（Markov Chain Monte Carlo）により、特定の確率分布に従い乱数（確率変数）を発生する方法について学べます。
掲載しているPythonコードを用いて、2状態の単純な例において乱数生成をシミュレートできます。
話の流れ
MCMCは複雑な確率分布に従う乱数を生成できる有効な手法です。その説明では遷移行列や採択確率の式が登場しますがこれらがどんな役割を果たすのかを見ていきます。
この記事では、2状態の単純な例を用いながら、「状態間の遷移件数をつり合わせる」という視点からMCMCの考え方を解説します。
以下の構成で話を進めていきます
- 確率密度関数に基づく分布を得る方法について
- 遷移行列と定常分布
- 詳細つり合い条件という考え方
- 仮の遷移行列と遷移の採択確率
- Pythonコードによるシミュレーション
確率密度関数に基づく分布を得るためには？
ポイント：確率密度関数から直接乱数を生成する方法とは
乱数を用いたシミュレーションでは、特定の確率分布に従う乱数を生成したい場面がよくあります。
例えば、指数分布
$$
p(x)=\lambda e^{-\lambda x}
\qquad (x \geq 0)
$$
に従う乱数を生成したいとします。
一般に、確率密度関数 (Probability Density Function, PDF) そのものは確率を与えるものではなく、ある区間に値が現れる確率は、PDFを積分することで得られます。
そこでまず累積分布関数
$$
F(x)=\int_{-\infty}^{x} p(t),dt
$$
を考えます。
指数分布の場合、
$$
F(x)=1-e^{-\lambda x}
$$
となります。
ここで、一様乱数 $u$ を区間 $[0,1]$ から生成し、
$$
u=F(x)
$$
を満たす $x$ を求めれば、その $x$ こそが「PDFから生成した確率変数」です。
指数分布では逆関数を解析的に求めることができ、
$$
x=F^{-1}(u)
=-\frac{1}{\lambda}\log(1-u)
$$
となるため、指数分布に従う乱数を容易に生成できますね。
しかし、現実に扱う確率分布は必ずしも単純ではないです。累積分布関数を解析的に計算できない場合や、逆関数を求められない場合も多いでしょう。
そのような場合には、上記の方法で直接乱数を生成することが難しいです。
MCMCは、この問題に対する一つの解決策です。累積分布関数やその逆関数を求める代わりに、ある状態から別の状態への遷移を繰り返すことで、目的の確率分布に従うサンプルを得ることができます。
遷移行列と定常分布
ポイント：行列を使って乱数を生成する
MCMCの考え方を理解するために、まずは状態が 2 つしかない単純な系を考えてみます。
状態「0」と状態「1」があり、ある時刻ごとに状態が変化するとします。このとき、状態の遷移確率をまとめたものが遷移行列です。
P=
\begin{pmatrix}
P_{00} & P_{01} \\
P_{10} & P_{11}
\end{pmatrix}
ここで、各成分は次のような意味を持ちます。
- $P_{00}$ : 状態0にいるとき、次の時刻も状態0にとどまる確率
- $P_{01}$ : 状態0にいるとき、次の時刻に状態1へ遷移する確率
- $P_{10}$ : 状態1にいるとき、次の時刻に状態0へ遷移する確率
- $P_{11}$ : 状態1にいるとき、次の時刻も状態1にとどまる確率
確率を表す行列なので、
P_{00}+P_{01}=1,
$$
P_{10}+P_{11}=1
$$
が成立します。
次に、ある時刻において各状態に存在する確率を
\boldsymbol{\pi}
=
\begin{pmatrix}
\pi_0 \\
\pi_1
\end{pmatrix}
というベクトルで表します。
例えば
\boldsymbol{\pi}
=
\begin{pmatrix}
0.3 \\
0.7
\end{pmatrix}
であれば、状態0にいる確率が 30 %、状態1にいる確率が 70 % であることを意味します。
このとき、1ステップ後の確率分布は
\boldsymbol{\pi}' = P^T\boldsymbol{\pi}
によって求められます。($P$の「転置」だからね！ただの$P$じゃないよ)
実際に計算すると、
\begin{pmatrix}
\pi_0' \\
\pi_1'
\end{pmatrix}
=
\begin{pmatrix}
P_{00} & P_{10}\\
P_{01} & P_{11}
\end{pmatrix}
\begin{pmatrix}
\pi_0 \\
\pi_1
\end{pmatrix}
となり、各状態へ流入する確率が自然に計算されます。
ここで、遷移件数に注目してみます。
例えば十分に大きな試行回数 $N$ を考えると、状態0に存在するサンプル数はおよそ
$$
N\pi_0
$$
です。
そのうち状態1へ移動する割合は $P_{01}$ なので、状態0から状態1への遷移件数は
$$
N\pi_0P_{01}
$$
と表せます。
同様に、状態1から状態0への遷移件数は
$$
N\pi_1P_{10}
$$
となります。
つまり、
$$
\pi_0P_{01}
$$
という量は「状態0から状態1への流れの大きさ」を表していると解釈できます。
この考え方は後に登場する詳細つり合い条件の理解において重要になります。
※念押しで確認しておくと、次ステップで状態0をとっている確率$\pi_0'$っていうのは、もともと状態0だったものがそのまま変わらなかった（＝発生確率$\pi_0 P_{00}$）ものと、状態1だったものが状態0に変化した（＝発生確率$\pi_1 P_{10}$）ものの和というわけだ。
$$
\pi_0' = \pi_0 P_{00} +\pi_1 P_{10}
$$
さて、遷移行列を繰り返し作用させていくと、ある確率分布へ収束することがあります。
特に、
P^T\boldsymbol{\pi}
=
\boldsymbol{\pi}
を満たす確率ベクトルが存在する場合、このベクトルを定常分布と呼びます。
定常分布では、状態遷移が起きているにもかかわらず、全体として見た状態の存在確率は変化しません。
言い換えると、状態間でサンプルが行き来していても、流入と流出がつり合っているため、分布の形が保たれている状態です。
MCMCでは、この定常分布を目的の確率分布に一致させることを目指します。そのため、まずは「どのような遷移を行えば望みの定常分布が得られるのか」を考える必要があります。
詳細つり合い条件という考え方
ポイント：状態変化の流れに注目する
前節では、遷移行列を繰り返し作用させることで定常分布が得られることを見ました。
それでは逆に、ある確率分布
\boldsymbol{\pi}
=
\begin{pmatrix}
\pi_0 \\
\pi_1
\end{pmatrix}
を定常分布として持つような遷移行列を作ることはできるのでしょうか。
状態数が少ない場合であれば、そのような遷移行列を構成することはそれほど難しくありません。しかし実際のMCMCでは、状態空間が非常に大きかったり、連続値を扱ったりすることが多く、目的の定常分布を持つ遷移行列を直接求めることは容易ではありません。
そこでMCMCでは、遷移行列そのものではなく、状態間の「流れ」に注目します。
前節で見たように、十分大きな試行回数 $N$ に対して、
$$
N\pi_0P_{01}
$$
は状態0から状態1への遷移件数を表し、
$$
N\pi_1P_{10}
$$
は状態1から状態0への遷移件数を表します。
もし定常分布が実現されているのであれば、状態0から状態1へ移動するサンプル数と、状態1から状態0へ移動するサンプル数は平均的に一致しているはずです。（前節の説明で言うと、「次ステップ」の分布が今の分布と変わらないということ）
すなわち、
N\pi_0P_{01}
=
N\pi_1P_{10}
が成立します。
両辺を $N$ で割ると、
\pi_0P_{01}
=
\pi_1P_{10}
を得ます。
この関係式は、状態0と状態1の間で流入と流出が完全につり合っていることを意味しています。
より一般には、
\pi_iP_{ij}
=
\pi_jP_{ji}
という形で表されます。
この条件を
詳細つり合い条件（Detailed Balance Condition）
と呼びます。
詳細つり合い条件が成立しているとき、各状態間の遷移件数が対ごとにつり合うため、結果として確率分布全体も変化しなくなります。つまり、詳細つり合い条件は定常分布を実現するための十分条件になっています。
MCMCでは、この詳細つり合い条件を満たすように状態遷移を設計します。
言い換えると、「目的の定常分布を持つ遷移行列を直接求める」のではなく、「状態間の流れがつり合うように遷移規則を作る」という方針を採用します。
次節では、この考え方を用いて、任意の確率分布を実現するための状態遷移の設計方法を見ていきます。
仮の遷移行列と遷移の採択確率
ポイント：「採択確率」を導入し状態遷移を制御する
詳細つり合い条件
\pi_i P_{ij}
=
\pi_j P_{ji}
が成立していれば、目的の確率分布 $\boldsymbol{\pi}$ は定常分布になります。
しかし、ここで一つ問題があります。
この条件を満たす遷移行列 $P$ を直接求めることは、一般には容易ではありません。特に状態空間が大きい場合や連続変数を扱う場合には、遷移行列そのものを構成することは現実的ではありません。
そこでMCMCでは、まず状態遷移の候補を与える仮の遷移行列
M=
\begin{pmatrix}
M_{00} & M_{01} \\
M_{10} & M_{11}
\end{pmatrix}
を用意します。
例えば、
M=
\begin{pmatrix}
0.3 & 0.7\\
0.5 & 0.5
\end{pmatrix}
のような行列を考えてみましょう。
この行列に従って状態遷移を行えば、状態0から状態1への遷移件数は
$$
N\pi_0M_{01}
$$
状態1から状態0への遷移件数は
$$
N\pi_1M_{10}
$$
となります。
一般には、
$$
\pi_0M_{01}
\neq
\pi_1M_{10}
$$
であり、詳細つり合い条件は成立しません。
つまり、このままでは目的の分布は定常分布になりません。
そこで導入されるのが採択確率です。
仮の遷移行列によって提案された状態遷移を必ず受け入れるのではなく、一定の確率で却下することを考えます。
例えば、
\pi_0M_{01}
>
\pi_1M_{10}
であったとします。
この場合、状態0から状態1への流れが大きすぎるため、詳細つり合い条件が崩れています。
そこで、状態0から状態1への遷移を確率 $\alpha_{01}$ でのみ受け入れることにします。
すると実際の遷移件数は
$$
N\pi_0M_{01}\alpha_{01}
$$
となります。
一方で逆方向の流れはそのままとすると、
$$
N\pi_1M_{10}
$$
です。
ここで
N\pi_0M_{01}\alpha_{01}
=
N\pi_1M_{10}
となるように $\alpha_{01}$ を選べば、両方向の遷移件数がつり合います。
したがって、
\alpha_{01}
=
\frac{\pi_1M_{10}}
{\pi_0M_{01}}
となります。
ただし採択確率は 1 を超えることができないため、一般には
\alpha_{01}
=
\min
\left(
1,
\frac{\pi_1M_{10}}
{\pi_0M_{01}}
\right)
と定義します。
同様に、任意の状態 $i,j$ に対して
\alpha_{ij}
=
\min
\left(
1,
\frac{\pi_jM_{ji}}
{\pi_iM_{ij}}
\right)
を用いることで、状態間の流れを補正できます。
このようにMCMCでは、仮の遷移行列そのものを変更するのではなく、提案された遷移に採択確率を導入することで詳細つり合い条件を実現します。
言い換えると、
- 仮の遷移行列は「どこへ移動したいか」を決める
- 採択確率は「流れの過不足を補正する」
という役割を担っています。
この仕組みによって、遷移行列を直接構成しなくても、目的とする確率分布を定常分布として持つマルコフ連鎖を作ることができます。
Pythonコードによるシミュレーション
ここまで見てきたように、MCMCでは仮の遷移行列と採択確率を用いることで、状態間の流れを調整し、目的の分布を定常分布として実現します。
例として2状態系の簡単なプログラムを用いて動作を確認してみます。
このプログラムでは、
-
pが目的分布 -
mが仮の遷移行列 -
p_adoptが採択確率
を表しています。
仮の遷移行列によって次の状態候補を決定し、その後で採択確率に基づいて実際に遷移するかどうかを判定しています。
コードを実行すると、状態0と状態1が出現した回数から経験的な確率分布が得られます。これを目的分布と比較することで、採択確率による補正が正しく機能していることを確認できます。
また、目的分布や仮の遷移行列を変更して実験してみると、仮の遷移行列が変わっても最終的には同じ分布へ収束することが分かります。一方で、収束までに必要なステップ数や状態の移動のしやすさは変化します。
この違いは、高次元問題において仮の遷移行列の設計が重要になる理由にもつながっています。
import numpy as np
n = 1000 #サンプリングの試行回数
p = np.array([4,9]) #目的分布（自由に設定する）
z = p.sum()
p = (1/z)*p #目的分布の規格化
# ↓これが「仮の遷移行列」（値を変えて実験してみよう）
m = (1/10)*np.array([[4,6],
[5,5]])
rng = np.random.default_rng(seed=1)
state = 0 #初期状態は｢0｣に設定
x = 1 #状態0を取った回数を記録
y = 0 #状態1を取った回数を記録
for i in range(1,n):
u = rng.random() #遷移先の決定に使う乱数
v = rng.random() #採択に使う乱数
if state == 0:
p_adopt = min(1, p[1]*m[1,0]/(p[0]*m[0,1]))
# ↑採択確率
if (u > m[0,0]) and (v < p_adopt):
y += 1
state = 1
else:
x += 1
else:
p_adopt = min(1, p[0]*m[0,1]/(p[1]*m[1,0]))
if (u < m[1,0] and (v < p_adopt)):
x += 1
state = 0
else:
y += 1
print(x/n, y/n) #サンプリング結果の度数分布を
print(p[0], p[1]) #目的分布と比較してみよう
# 目的分布、仮の遷移号列を変えて実験してみよう
まとめ
この記事では、MCMCによる乱数生成の考え方を2状態系の例を用いて見てきました。
まず、確率密度関数から直接乱数を生成するには、累積分布関数やその逆関数が必要になることを確認しました。しかし、複雑な分布ではそれらを求めることが難しい場合があります。
そこでMCMCでは、目的の分布を定常分布として持つマルコフ連鎖を構成するという考え方を採用します。
また、遷移行列そのものを直接求める代わりに、状態間の遷移件数がつり合うという詳細つり合い条件に注目しました。そして、仮の遷移行列によって提案された遷移を採択確率で補正することで、この条件を実現できることを見ました。
言い換えると、MCMCは仮の遷移行列を更新しながら分布を作る手法ではなく、仮の遷移規則に対して適切な採択確率を導入することで、目的の分布へ収束する状態遷移を実現する手法だと捉えることができます。
今回扱ったのは状態が2つしかない非常に単純な例でしたが、詳細つり合い条件と採択確率という基本的な考え方は、より高次元で複雑な問題に対してもそのまま用いられています。