マルコフ連鎖は、「未来の状態は現在の状態のみに依存し、過去の履歴には依存しない」という マルコフ性 を持つ確率過程です。Webページのランキング（PageRank）、自然言語処理、待ち行列理論、統計力学など、非常に幅広い分野で活用されています。
マルコフ連鎖の長期的な振る舞いを理解する上で中心的な役割を果たすのが 定常分布（stationary distribution） です。定常分布は、「十分長い時間が経過した後に各状態にいる確率」を表します。本記事では、遷移確率行列の性質から出発して、Chapman-Kolmogorov方程式、定常分布の存在条件と一意性、エルゴード定理、さらに吸収マルコフ連鎖までを丁寧に導出し、Pythonで実装・可視化します。
本記事の内容
- 遷移確率行列 $\bm{P}$ の性質（確率的行列）
- Chapman-Kolmogorov方程式の導出
- 定常分布 $\bm{\pi} \bm{P} = \bm{\pi}$ の導出と存在条件
- 既約性・非周期性・正再帰性
- エルゴード定理の証明スケッチ
- 定常分布の一意性の証明
- 吸収マルコフ連鎖
- Pythonでの遷移行列の固有値解析・定常分布計算・収束シミュレーション
前提知識
この記事を読む前に、以下の記事を読んでおくと理解が深まります。
遷移確率行列の定義と性質
マルコフ連鎖の定義
状態空間 $S = \{1, 2, \dots, N\}$ を持つ離散時間マルコフ連鎖 $\{X_n\}_{n \geq 0}$ は、マルコフ性
$$ P(X_{n+1} = j \mid X_n = i, X_{n-1} = i_{n-1}, \dots, X_0 = i_0) = P(X_{n+1} = j \mid X_n = i) $$
を満たす確率過程です。
遷移確率と遷移確率行列
1ステップ遷移確率 を
$$ p_{ij} = P(X_{n+1} = j \mid X_n = i) $$
と定義します。時間に依存しない場合（時間斉次（time-homogeneous） な場合）、これを $N \times N$ の行列にまとめたものが 遷移確率行列（transition probability matrix） です。
$$ \bm{P} = \begin{pmatrix} p_{11} & p_{12} & \cdots & p_{1N} \\ p_{21} & p_{22} & \cdots & p_{2N} \\ \vdots & \vdots & \ddots & \vdots \\ p_{N1} & p_{N2} & \cdots & p_{NN} \end{pmatrix} $$
確率的行列の性質
遷移確率行列 $\bm{P}$ は以下の性質を持ちます。
性質1: 非負性
$$ p_{ij} \geq 0 \quad \forall i, j $$
確率は非負であるため。
性質2: 行和が1（確率的行列, stochastic matrix）
$$ \sum_{j=1}^{N} p_{ij} = 1 \quad \forall i $$
状態 $i$ から出発して、必ずどこかの状態に遷移するためです。この性質を持つ行列を 確率的行列（stochastic matrix） または 確率行列 と呼びます。より正確には 行確率的行列（row stochastic matrix） です。
ベクトル表記で書くと、$\bm{1} = (1, 1, \dots, 1)^T$ として、
$$ \bm{P} \bm{1} = \bm{1} $$
つまり、$\bm{1}$ は $\bm{P}$ の固有値1に対応する右固有ベクトルです。
性質3: 固有値の範囲
確率的行列の固有値 $\lambda$ はすべて $|\lambda| \leq 1$ を満たします。
証明: $\bm{P}\bm{v} = \lambda \bm{v}$ とし、$|v_k| = \max_i |v_i|$ とします。$k$ 番目の成分を見ると、
$$ \lambda v_k = \sum_{j=1}^{N} p_{kj} v_j $$
両辺の絶対値を取ると、
$$ \begin{align} |\lambda| |v_k| &= \left|\sum_{j=1}^{N} p_{kj} v_j\right| \\ &\leq \sum_{j=1}^{N} p_{kj} |v_j| \quad (\text{三角不等式}) \\ &\leq \sum_{j=1}^{N} p_{kj} |v_k| \quad (\because |v_j| \leq |v_k|) \\ &= |v_k| \quad (\because \text{行和} = 1) \end{align} $$
$|v_k| > 0$ で割ると $|\lambda| \leq 1$ が得られます。$\square$
Chapman-Kolmogorov方程式
nステップ遷移確率
nステップ遷移確率 を
$$ p_{ij}^{(n)} = P(X_{m+n} = j \mid X_m = i) $$
と定義します。時間斉次性から、これは $m$ に依存しません。
主張: $p_{ij}^{(n)}$ は遷移確率行列の $n$ 乗の $(i, j)$ 成分に等しい。
$$ p_{ij}^{(n)} = (\bm{P}^n)_{ij} $$
Chapman-Kolmogorov方程式の導出
任意の $0 \leq m \leq n$ に対して、全確率の公式を用いると、
$$ \begin{align} p_{ij}^{(n)} &= P(X_n = j \mid X_0 = i) \\ &= \sum_{k \in S} P(X_n = j, X_m = k \mid X_0 = i) \\ &= \sum_{k \in S} P(X_n = j \mid X_m = k, X_0 = i) \cdot P(X_m = k \mid X_0 = i) \\ &= \sum_{k \in S} P(X_n = j \mid X_m = k) \cdot P(X_m = k \mid X_0 = i) \quad (\because \text{マルコフ性}) \\ &= \sum_{k \in S} p_{kj}^{(n-m)} \cdot p_{ik}^{(m)} \end{align} $$
これを Chapman-Kolmogorov方程式 と呼びます。行列の形で書くと、
$$ \bm{P}^{(n)} = \bm{P}^{(m)} \cdot \bm{P}^{(n-m)} $$
特に $m$ ステップと $(n-m)$ ステップの遷移を合成すると $n$ ステップの遷移になるという直感的な意味を持ちます。$m = 1$ の場合は、
$$ \bm{P}^{(n)} = \bm{P} \cdot \bm{P}^{(n-1)} = \bm{P}^n $$
となり、$n$ ステップ遷移確率行列は $\bm{P}$ の $n$ 乗であることが帰納法で示されます。
状態確率ベクトルの時間発展
時刻 $n$ での状態確率ベクトルを
$$ \bm{\pi}^{(n)} = (\pi_1^{(n)}, \pi_2^{(n)}, \dots, \pi_N^{(n)}) $$
ここで $\pi_j^{(n)} = P(X_n = j)$ とします。全確率の公式より、
$$ \pi_j^{(n)} = \sum_{i \in S} \pi_i^{(0)} p_{ij}^{(n)} $$
行列表記では、
$$ \bm{\pi}^{(n)} = \bm{\pi}^{(0)} \bm{P}^n $$
初期分布 $\bm{\pi}^{(0)}$ に遷移確率行列を $n$ 回掛けると、時刻 $n$ での分布が得られます。
定常分布
定常分布の定義
定義: 確率ベクトル $\bm{\pi} = (\pi_1, \pi_2, \dots, \pi_N)$ が 定常分布（stationary distribution） であるとは、以下の条件を満たすことをいいます。
$$ \bm{\pi} \bm{P} = \bm{\pi}, \quad \pi_i \geq 0, \quad \sum_{i=1}^{N} \pi_i = 1 $$
つまり、$\bm{\pi}$ は $\bm{P}$ の 左固有ベクトル であり、固有値1に対応します。
直感的には、状態確率分布が $\bm{\pi}$ であるとき、1ステップ遷移した後も分布が $\bm{\pi}$ のまま変わらないということです。「平衡状態」と解釈できます。
成分ごとの表記
$$ \pi_j = \sum_{i=1}^{N} \pi_i p_{ij} \quad \forall j $$
これは「状態 $j$ にいる確率は、すべての状態 $i$ から $j$ に遷移する確率の加重和に等しい」ことを意味します。入ってくるフロー（流入）と出ていくフロー（流出）が釣り合っています。
定常分布の求め方
$\bm{\pi} \bm{P} = \bm{\pi}$ を変形すると、
$$ \bm{\pi}(\bm{P} – \bm{I}) = \bm{0} $$
これは $(\bm{P}^T – \bm{I})\bm{\pi}^T = \bm{0}$ と等価です。つまり $\bm{\pi}^T$ は $\bm{P}^T$ の固有値1に対応する固有ベクトルです。
これに正規化条件 $\sum_i \pi_i = 1$ を加えた連立方程式を解くことで定常分布が求まります。具体的には、$(\bm{P}^T – \bm{I})$ の1行を $\bm{1}^T$（全成分1のベクトル）に置き換えた連立方程式を解きます。
既約性・非周期性・正再帰性
定常分布の存在と一意性を保証するために、マルコフ連鎖の構造に関する3つの重要な概念を導入します。
既約性（Irreducibility）
定義: マルコフ連鎖が 既約（irreducible） であるとは、任意の状態ペア $(i, j)$ に対して、ある $n \geq 1$ が存在して $p_{ij}^{(n)} > 0$ であることをいいます。
直感的には、どの状態からどの状態へも有限ステップで到達できるということです。状態空間が1つの「連結成分」であることを意味します。
対偶: 既約でない場合、状態空間は複数の 閉じた連絡クラス（communicating class） に分割されます。
非周期性（Aperiodicity）
定義: 状態 $i$ の 周期 $d(i)$ を
$$ d(i) = \gcd\{n \geq 1 : p_{ii}^{(n)} > 0\} $$
と定義します（gcd は最大公約数）。$d(i) = 1$ のとき状態 $i$ は 非周期的（aperiodic） であるといいます。
既約マルコフ連鎖では、すべての状態の周期は同じです。すべての状態が非周期的であるとき、マルコフ連鎖は非周期的であるといいます。
例: 状態 $\{0, 1\}$ で $p_{01} = p_{10} = 1$ のマルコフ連鎖は周期 $d = 2$ です。偶数ステップでのみ元の状態に戻るためです。
正再帰性（Positive Recurrence）
定義: 状態 $i$ に出発して再び $i$ に戻るまでの時間（初回帰時間, first return time）を
$$ T_i = \min\{n \geq 1 : X_n = i \mid X_0 = i\} $$
とします。
- $P(T_i < \infty) = 1$ のとき、状態 $i$ は 再帰的（recurrent）
- $E[T_i] < \infty$ のとき、状態 $i$ は 正再帰的（positive recurrent）
- $E[T_i] = \infty$ のとき、状態 $i$ は 零再帰的（null recurrent）
有限状態空間では、既約マルコフ連鎖は常に正再帰的です。
エルゴード定理
定理の主張
定理（マルコフ連鎖のエルゴード定理）: 既約かつ非周期的かつ正再帰的なマルコフ連鎖に対して、一意の定常分布 $\bm{\pi}$ が存在し、初期分布によらず、
$$ \lim_{n \to \infty} p_{ij}^{(n)} = \pi_j \quad \forall i, j $$
が成り立つ。行列の形で書くと、
$$ \lim_{n \to \infty} \bm{P}^n = \bm{1} \bm{\pi} = \begin{pmatrix} \pi_1 & \pi_2 & \cdots & \pi_N \\ \pi_1 & \pi_2 & \cdots & \pi_N \\ \vdots & \vdots & \ddots & \vdots \\ \pi_1 & \pi_2 & \cdots & \pi_N \end{pmatrix} $$
さらに、定常分布と平均初回帰時間の間に
$$ \pi_i = \frac{1}{E[T_i]} $$
の関係が成り立ちます。
証明スケッチ
ステップ1: 定常分布の存在
有限既約正再帰マルコフ連鎖に対して、各状態 $j$ について、
$$ \pi_j = \frac{1}{E[T_j]} $$
と定義します。ここで $E[T_j]$ は状態 $j$ への平均初回帰時間です。正再帰性から $E[T_j] < \infty$ なので $\pi_j > 0$ です。
再生定理（Renewal Theory） により、状態 $j$ を訪問する長期的な頻度は $1 / E[T_j]$ に等しくなります。具体的には、
$$ \frac{1}{n} \sum_{k=0}^{n-1} \mathbf{1}(X_k = j) \to \frac{1}{E[T_j]} = \pi_j \quad (n \to \infty) $$
が確率1で成り立ちます。
$\bm{\pi}$ が定常分布であることの確認: $\bm{\pi} \bm{P} = \bm{\pi}$ を示す必要があります。
$$ (\bm{\pi} \bm{P})_j = \sum_i \pi_i p_{ij} $$
再生的な議論から、時刻 $n$ で状態 $j$ にいる確率は、時刻 $n – 1$ で各状態 $i$ にいて1ステップで $j$ に遷移する確率の和であり、$n \to \infty$ では定常分布に収束するため、
$$ \pi_j = \sum_i \pi_i p_{ij} $$
が成り立ちます。
ステップ2: 一意性の証明
$\bm{\pi}$ と $\bm{\pi}’$ が共に定常分布であるとします。$\bm{\pi}’ \bm{P}^n = \bm{\pi}’$ がすべての $n$ で成り立つので、
$$ \pi_j’ = \sum_i \pi_i’ p_{ij}^{(n)} $$
$n \to \infty$ で $p_{ij}^{(n)} \to \pi_j$（ステップ1で示した収束）を用いると、
$$ \pi_j’ = \sum_i \pi_i’ \pi_j = \pi_j \sum_i \pi_i’ = \pi_j \cdot 1 = \pi_j $$
よって $\bm{\pi}’ = \bm{\pi}$ であり、定常分布は一意です。$\square$
ステップ3: 非周期性の役割
非周期性の条件を外すと、$p_{ij}^{(n)}$ は $n \to \infty$ で振動し、極限を持たない場合があります（周期的挙動）。しかし、Cesaro平均 の意味では収束します。
$$ \lim_{n \to \infty} \frac{1}{n} \sum_{k=0}^{n-1} p_{ij}^{(k)} = \pi_j $$
詳細釣り合い条件
定常分布を求める上で有用な十分条件として 詳細釣り合い条件（detailed balance condition） があります。
定義: 確率ベクトル $\bm{\pi}$ が 詳細釣り合い条件 を満たすとは、すべての状態ペア $(i, j)$ に対して、
$$ \pi_i p_{ij} = \pi_j p_{ji} $$
が成り立つことをいいます。
命題: 詳細釣り合い条件を満たす $\bm{\pi}$ は定常分布である。
証明:
$$ \begin{align} (\bm{\pi} \bm{P})_j &= \sum_i \pi_i p_{ij} \\ &= \sum_i \pi_j p_{ji} \quad (\because \text{詳細釣り合い条件}) \\ &= \pi_j \sum_i p_{ji} \\ &= \pi_j \cdot 1 \quad (\because \text{行和} = 1) \\ &= \pi_j \end{align} $$
よって $\bm{\pi} \bm{P} = \bm{\pi}$ が成り立ちます。$\square$
詳細釣り合い条件は、MCMCのメトロポリス=ヘイスティングス法の設計原理としても重要です。
吸収マルコフ連鎖
定義
状態 $i$ が 吸収状態（absorbing state） であるとは、$p_{ii} = 1$ であることをいいます。一度入ると二度と出られない状態です。
少なくとも1つの吸収状態を持ち、すべての非吸収状態から少なくとも1つの吸収状態に到達可能なマルコフ連鎖を 吸収マルコフ連鎖（absorbing Markov chain） と呼びます。
標準形
吸収マルコフ連鎖の遷移確率行列は、状態の順序を入れ替えて以下の標準形に書くことができます。
$$ \bm{P} = \begin{pmatrix} \bm{Q} & \bm{R} \\ \bm{0} & \bm{I} \end{pmatrix} $$
ここで、
- $\bm{Q}$: 非吸収状態間の遷移確率（$t \times t$ 行列、$t$ は非吸収状態の数）
- $\bm{R}$: 非吸収状態から吸収状態への遷移確率
- $\bm{I}$: 吸収状態の単位行列
- $\bm{0}$: 零行列（吸収状態から非吸収状態へは遷移しない）
基本行列
基本行列（fundamental matrix） $\bm{N}$ を
$$ \bm{N} = (\bm{I} – \bm{Q})^{-1} = \sum_{k=0}^{\infty} \bm{Q}^k $$
と定義します。$\bm{Q}$ のすべての固有値の絶対値が1未満であるため、このNeumann級数は収束します。
$\bm{N}$ の $(i, j)$ 成分 $n_{ij}$ は、「状態 $i$ から出発して吸収されるまでに状態 $j$ に滞在する期待回数」を表します。
吸収までの期待ステップ数: 状態 $i$ から出発して吸収されるまでの期待ステップ数は、
$$ E[T_{\text{abs}} \mid X_0 = i] = \sum_j n_{ij} = (\bm{N} \bm{1})_i $$
吸収確率: 状態 $i$ から出発して吸収状態 $j$ に吸収される確率は、$(\bm{N}\bm{R})_{ij}$ で与えられます。
Pythonでの実装
定常分布の計算
import numpy as np
import matplotlib.pyplot as plt
from scipy import linalg
# --- 3状態マルコフ連鎖の遷移確率行列 ---
P = np.array([
[0.1, 0.6, 0.3],
[0.4, 0.2, 0.4],
[0.3, 0.3, 0.4]
])
# 行和が1であることの確認
print("行和:", P.sum(axis=1))
# --- 方法1: 固有値分解による定常分布の計算 ---
# P^T の固有値1に対応する固有ベクトルを求める
eigenvalues, eigenvectors = linalg.eig(P.T)
print("\n固有値:", eigenvalues)
# 固有値1に最も近い固有ベクトルを選択
idx = np.argmin(np.abs(eigenvalues - 1.0))
pi_eigen = np.real(eigenvectors[:, idx])
pi_eigen = pi_eigen / pi_eigen.sum() # 正規化
print("定常分布（固有値分解）:", pi_eigen)
# --- 方法2: 連立方程式を解く ---
# πP = π かつ Σπ_i = 1
# (P^T - I)π^T = 0 の1行を正規化条件に置き換え
A = P.T - np.eye(3)
A[-1, :] = 1 # 最後の行を正規化条件に置き換え
b = np.zeros(3)
b[-1] = 1
pi_solve = linalg.solve(A, b)
print("定常分布（連立方程式）:", pi_solve)
# --- 方法3: べき乗法（Power Iteration） ---
pi_power = np.array([1/3, 1/3, 1/3]) # 初期分布
for _ in range(1000):
pi_power = pi_power @ P
print("定常分布（べき乗法）:", pi_power)
# 検証: πP = π
print("\n検証 πP:", pi_eigen @ P)
print("検証 π: ", pi_eigen)
print("差のノルム:", np.linalg.norm(pi_eigen @ P - pi_eigen))
収束シミュレーション
import numpy as np
import matplotlib.pyplot as plt
# 遷移確率行列
P = np.array([
[0.1, 0.6, 0.3],
[0.4, 0.2, 0.4],
[0.3, 0.3, 0.4]
])
# 定常分布の計算
from scipy import linalg
A = P.T - np.eye(3)
A[-1, :] = 1
b = np.zeros(3)
b[-1] = 1
pi_stationary = linalg.solve(A, b)
# --- 異なる初期分布からの収束 ---
initial_distributions = [
np.array([1.0, 0.0, 0.0]), # 状態0から開始
np.array([0.0, 1.0, 0.0]), # 状態1から開始
np.array([0.0, 0.0, 1.0]), # 状態2から開始
np.array([0.5, 0.3, 0.2]), # 混合初期分布
]
labels = ['Start: State 0', 'Start: State 1', 'Start: State 2', 'Start: Mixed']
colors_list = ['tab:blue', 'tab:orange', 'tab:green', 'tab:red']
n_steps = 30
fig, axes = plt.subplots(1, 3, figsize=(16, 5))
state_names = ['State 0', 'State 1', 'State 2']
for state_idx in range(3):
ax = axes[state_idx]
for init_idx, pi0 in enumerate(initial_distributions):
probs = np.zeros(n_steps + 1)
pi_n = pi0.copy()
for n in range(n_steps + 1):
probs[n] = pi_n[state_idx]
if n < n_steps:
pi_n = pi_n @ P
ax.plot(range(n_steps + 1), probs, color=colors_list[init_idx],
linewidth=1.5, label=labels[init_idx])
ax.axhline(y=pi_stationary[state_idx], color='k', linestyle='--',
linewidth=2, label=f'π = {pi_stationary[state_idx]:.4f}')
ax.set_title(f"Convergence to π({state_names[state_idx]})")
ax.set_xlabel("Step n")
ax.set_ylabel("Probability")
ax.legend(fontsize=8)
ax.grid(True, alpha=0.3)
ax.set_ylim(0, 1.05)
plt.tight_layout()
plt.savefig("markov_convergence.png", dpi=150, bbox_inches="tight")
plt.show()
P^n の可視化
遷移確率行列の $n$ 乗が定常分布行列に収束する様子を可視化します。
import numpy as np
import matplotlib.pyplot as plt
P = np.array([
[0.1, 0.6, 0.3],
[0.4, 0.2, 0.4],
[0.3, 0.3, 0.4]
])
# P^n を計算して各成分の推移を描画
n_max = 30
P_n_values = np.zeros((n_max + 1, 3, 3))
P_n = np.eye(3)
for n in range(n_max + 1):
P_n_values[n] = P_n
P_n = P_n @ P
fig, axes = plt.subplots(3, 3, figsize=(14, 12))
for i in range(3):
for j in range(3):
axes[i, j].plot(range(n_max + 1), P_n_values[:, i, j],
'b-', linewidth=1.5)
axes[i, j].axhline(y=pi_stationary[j], color='r', linestyle='--',
linewidth=1.5)
axes[i, j].set_title(f"$P^n$[{i},{j}] → π_{j} = {pi_stationary[j]:.4f}")
axes[i, j].set_xlabel("n")
axes[i, j].set_ylim(0, 1)
axes[i, j].grid(True, alpha=0.3)
plt.suptitle("Convergence of $P^n$ to Stationary Distribution", fontsize=14, y=1.01)
plt.tight_layout()
plt.savefig("Pn_convergence.png", dpi=150, bbox_inches="tight")
plt.show()
固有値解析
import numpy as np
import matplotlib.pyplot as plt
from scipy import linalg
P = np.array([
[0.1, 0.6, 0.3],
[0.4, 0.2, 0.4],
[0.3, 0.3, 0.4]
])
# 固有値と固有ベクトル
eigenvalues, left_eigvecs = linalg.eig(P, left=True, right=False)
print("=== 固有値解析 ===")
print(f"固有値: {eigenvalues}")
print(f"固有値の絶対値: {np.abs(eigenvalues)}")
print()
# 固有値を複素平面上にプロット
fig, axes = plt.subplots(1, 2, figsize=(14, 5))
# 固有値の配置
theta = np.linspace(0, 2 * np.pi, 100)
axes[0].plot(np.cos(theta), np.sin(theta), 'k--', alpha=0.3, label='Unit circle')
axes[0].scatter(np.real(eigenvalues), np.imag(eigenvalues),
s=100, c='red', zorder=5, label='Eigenvalues')
for i, ev in enumerate(eigenvalues):
axes[0].annotate(f'λ_{i+1}={ev:.4f}',
(np.real(ev), np.imag(ev)),
textcoords="offset points", xytext=(10, 10), fontsize=10)
axes[0].set_xlim(-1.3, 1.3)
axes[0].set_ylim(-1.3, 1.3)
axes[0].set_aspect('equal')
axes[0].axhline(y=0, color='k', linewidth=0.5)
axes[0].axvline(x=0, color='k', linewidth=0.5)
axes[0].set_title("Eigenvalues of P in Complex Plane")
axes[0].set_xlabel("Real")
axes[0].set_ylabel("Imaginary")
axes[0].legend()
axes[0].grid(True, alpha=0.3)
# 収束速度: |λ_2|^n
lambda_sorted = sorted(np.abs(eigenvalues), reverse=True)
lambda_2 = lambda_sorted[1] # 第2固有値の絶対値
n_vals = np.arange(0, 31)
convergence_rate = lambda_2 ** n_vals
axes[1].semilogy(n_vals, convergence_rate, 'b-', linewidth=2,
label=f'|λ₂|^n = {lambda_2:.4f}^n')
axes[1].set_title("Convergence Rate: $|\\lambda_2|^n$")
axes[1].set_xlabel("Step n")
axes[1].set_ylabel("$|\\lambda_2|^n$")
axes[1].legend()
axes[1].grid(True, alpha=0.3)
plt.tight_layout()
plt.savefig("eigenvalue_analysis.png", dpi=150, bbox_inches="tight")
plt.show()
print(f"第2固有値の絶対値 |λ₂| = {lambda_2:.6f}")
print(f"収束の半減期: n ≈ {-np.log(2)/np.log(lambda_2):.1f} ステップ")
収束速度は第2固有値の絶対値 $|\lambda_2|$ で決まります。$|\lambda_2|$ が1に近いほど収束が遅く、0に近いほど速やかに定常分布に収束します。具体的には、定常分布との誤差は $O(|\lambda_2|^n)$ で減少します。
吸収マルコフ連鎖
import numpy as np
import matplotlib.pyplot as plt
# ギャンブラーの破産問題
# 所持金が0（破産）または N（目標達成）になったら終了
# 各ステップで確率pで+1、確率1-pで-1
N = 10 # 目標金額
p = 0.4 # 勝つ確率
# 状態空間: 0, 1, 2, ..., N
# 吸収状態: 0（破産）, N（目標達成）
# 非吸収状態: 1, 2, ..., N-1
t = N - 1 # 非吸収状態の数
# Q行列（非吸収状態間の遷移）
Q = np.zeros((t, t))
for i in range(t):
if i > 0:
Q[i, i - 1] = 1 - p # 所持金が1減る
if i < t - 1:
Q[i, i + 1] = p # 所持金が1増える
# R行列（非吸収状態→吸収状態への遷移）
R = np.zeros((t, 2))
R[0, 0] = 1 - p # 状態1→状態0（破産）
R[t - 1, 1] = p # 状態N-1→状態N（目標達成）
# 基本行列
I_t = np.eye(t)
N_fund = np.linalg.inv(I_t - Q)
# 吸収までの期待ステップ数
expected_steps = N_fund @ np.ones(t)
# 各吸収状態への吸収確率
absorption_probs = N_fund @ R
print("=== ギャンブラーの破産問題 ===")
print(f"目標金額: {N}, 勝つ確率: {p}")
print()
print("初期所持金 | 破産確率 | 目標達成確率 | 期待ステップ数")
print("-" * 55)
for i in range(t):
start = i + 1 # 所持金（1-indexed）
print(f" {start:2d} | {absorption_probs[i, 0]:.4f} | {absorption_probs[i, 1]:.4f} | {expected_steps[i]:.1f}")
# 可視化
fig, axes = plt.subplots(1, 3, figsize=(16, 5))
states = np.arange(1, N)
# (1) 吸収確率
axes[0].bar(states - 0.15, absorption_probs[:, 0], width=0.3,
color='tab:red', alpha=0.8, label='Ruin (state 0)')
axes[0].bar(states + 0.15, absorption_probs[:, 1], width=0.3,
color='tab:green', alpha=0.8, label=f'Goal (state {N})')
axes[0].set_title(f"Absorption Probabilities (p={p})")
axes[0].set_xlabel("Initial fortune")
axes[0].set_ylabel("Probability")
axes[0].legend()
axes[0].grid(True, alpha=0.3, axis='y')
# (2) 期待ステップ数
axes[1].bar(states, expected_steps, color='steelblue', alpha=0.8)
axes[1].set_title("Expected Steps to Absorption")
axes[1].set_xlabel("Initial fortune")
axes[1].set_ylabel("Expected steps")
axes[1].grid(True, alpha=0.3, axis='y')
# (3) サンプルパスのシミュレーション
np.random.seed(42)
n_simulations = 10
for sim in range(n_simulations):
fortune = 5 # 初期所持金
path = [fortune]
while 0 < fortune < N:
if np.random.random() < p:
fortune += 1
else:
fortune -= 1
path.append(fortune)
color = 'tab:red' if fortune == 0 else 'tab:green'
axes[2].plot(path, color=color, alpha=0.6, linewidth=0.8)
axes[2].axhline(y=0, color='tab:red', linestyle='--', linewidth=1)
axes[2].axhline(y=N, color='tab:green', linestyle='--', linewidth=1)
axes[2].set_title(f"Sample Paths (start=5, p={p})")
axes[2].set_xlabel("Step")
axes[2].set_ylabel("Fortune")
axes[2].set_ylim(-0.5, N + 0.5)
axes[2].grid(True, alpha=0.3)
plt.tight_layout()
plt.savefig("absorbing_markov_chain.png", dpi=150, bbox_inches="tight")
plt.show()
まとめ
本記事では、マルコフ連鎖の定常分布と遷移確率行列について包括的に解説しました。
- 遷移確率行列 $\bm{P}$ は確率的行列であり、固有値はすべて $|\lambda| \leq 1$ を満たす
- Chapman-Kolmogorov方程式により、$n$ ステップ遷移確率は $\bm{P}^n$ で計算される
- 定常分布 $\bm{\pi}$ は $\bm{\pi}\bm{P} = \bm{\pi}$ を満たす左固有ベクトルである
- 既約性・非周期性・正再帰性のもとでエルゴード定理が成り立ち、定常分布への収束が保証される
- 収束速度は第2固有値 $|\lambda_2|$ で決まり、$|\lambda_2|$ が小さいほど速く収束する
- 詳細釣り合い条件は定常分布であるための十分条件であり、MCMCの設計原理となる
- 吸収マルコフ連鎖では基本行列 $\bm{N} = (\bm{I} – \bm{Q})^{-1}$ を用いて吸収確率や期待ステップ数を解析的に計算できる
次のステップとして、以下の記事も参考にしてください。