マルコフ連鎖が長時間経過した後にどのような分布に落ち着くかを記述するのが定常分布です。定常分布の存在と一意性は、マルコフ連鎖の理論で最も重要な結果の一つであり、PageRankアルゴリズムやMCMC法の理論的基盤でもあります。
本記事では、定常分布の定義、存在条件、計算法を解説し、PageRankとの関係をPythonで実装します。
本記事の内容
- 定常分布 $\bm{\pi}\bm{P} = \bm{\pi}$ の定義
- 存在条件（既約・非周期）
- 定常分布の計算法
- エルゴード定理
- PageRankとの関係
- Pythonでの実装と可視化
前提知識
この記事を読む前に、以下の記事を読んでおくと理解が深まります。
定常分布の定義
確率ベクトル $\bm{\pi} = (\pi_1, \pi_2, \dots, \pi_N)$ が遷移確率行列 $\bm{P}$ の定常分布（stationary distribution）であるとは、以下の条件を満たすことです。
$$ \boxed{\bm{\pi}\bm{P} = \bm{\pi}} $$
成分で書くと：
$$ \pi_j = \sum_{i=1}^{N} \pi_i p_{ij}, \quad j = 1, 2, \dots, N $$
さらに確率分布であるための条件：
$$ \pi_i \geq 0, \quad \sum_{i=1}^{N} \pi_i = 1 $$
定常分布の意味
初期分布が $\bm{\pi}^{(0)} = \bm{\pi}$ のとき：
$$ \bm{\pi}^{(1)} = \bm{\pi}\bm{P} = \bm{\pi} $$
つまり、定常分布を初期分布とすると、すべての時刻で分布が変わりません。これが「定常」と呼ばれる所以です。
固有値との関係
$\bm{\pi}\bm{P} = \bm{\pi}$ は $\bm{P}^T\bm{\pi}^T = \bm{\pi}^T$ と同値です。つまり、$\bm{\pi}^T$ は $\bm{P}^T$ の固有値1に対応する左固有ベクトルです。
確率行列 $\bm{P}$ の固有値は以下の性質を持ちます。
- 最大固有値は1: $\lambda_1 = 1$
- すべての固有値の絶対値は1以下: $|\lambda_i| \leq 1$
存在と一意性の条件
既約性（Irreducibility）
すべての状態が互いに到達可能であるとき、マルコフ連鎖は既約と呼ばれます。
$$ \forall i, j \in \mathcal{S}, \exists n \geq 1 : p_{ij}^{(n)} > 0 $$
非周期性（Aperiodicity）
すべての状態の周期が1であるとき、非周期と呼ばれます。
$$ d_i = \gcd\{n \geq 1 : p_{ii}^{(n)} > 0\} = 1 $$
十分条件: ある状態 $i$ で $p_{ii} > 0$（自己ループあり）ならば、既約なマルコフ連鎖は非周期です。
定常分布の存在定理
有限状態の既約マルコフ連鎖には、一意の定常分布が存在します。
$$ \text{既約} \implies \exists ! \bm{\pi} : \bm{\pi}\bm{P} = \bm{\pi}, \quad \pi_i > 0 $$
エルゴード定理
既約かつ非周期（エルゴード的）な有限状態マルコフ連鎖では、初期分布によらず分布が定常分布に収束します。
$$ \boxed{\lim_{n \to \infty} p_{ij}^{(n)} = \pi_j \quad \text{（すべての } i, j \text{）}} $$
つまり $\lim_{n \to \infty} \bm{P}^n$ の各行がすべて $\bm{\pi}$ に等しくなります。
さらに、時間平均も定常分布に収束します（強エルゴード性）。
$$ \frac{1}{n}\sum_{k=0}^{n-1} \mathbf{1}_{X_k = j} \xrightarrow{a.s.} \pi_j $$
定常分布の計算法
方法1: 連立方程式を解く
$\bm{\pi}\bm{P} = \bm{\pi}$ と $\sum_i \pi_i = 1$ から連立方程式を立てて解きます。
3状態の例（$\bm{\pi}\bm{P} = \bm{\pi}$）：
$$ \begin{cases} \pi_1 p_{11} + \pi_2 p_{21} + \pi_3 p_{31} = \pi_1 \\ \pi_1 p_{12} + \pi_2 p_{22} + \pi_3 p_{32} = \pi_2 \\ \pi_1 p_{13} + \pi_2 p_{23} + \pi_3 p_{33} = \pi_3 \\ \pi_1 + \pi_2 + \pi_3 = 1 \end{cases} $$
3つ目の方程式は最初の2つから導出されるため、3つ目を $\sum_i \pi_i = 1$ で置き換えます。
方法2: 固有値分解
$\bm{P}^T$ の固有値1に対応する固有ベクトルを求め、正規化します。
方法3: $\bm{P}^n$ の極限
$\bm{P}^n$ を大きな $n$ で計算し、各行を読み取ります（エルゴード的な場合）。
方法4: 詳細釣り合い条件
以下の条件を満たす $\bm{\pi}$ は定常分布です（十分条件）。
$$ \pi_i p_{ij} = \pi_j p_{ji} \quad \text{（すべての } i, j \text{）} $$
これを詳細釣り合い条件（detailed balance）と呼びます。この条件を満たすマルコフ連鎖を可逆と呼びます。
検証：両辺を $j$ について和を取ると $\pi_i \sum_j p_{ij} = \sum_j \pi_j p_{ji}$ より $\pi_i = \sum_j \pi_j p_{ji}$、これは $\bm{\pi}\bm{P} = \bm{\pi}$ と同値です。
PageRankとの関係
GoogleのPageRankアルゴリズムは、Webページのリンク構造をマルコフ連鎖としてモデル化し、その定常分布をページの重要度として使います。
モデル
$N$ 個のWebページがあり、ページ $i$ から出ているリンクの数を $L_i$ とします。ランダムサーファーが各リンクを等確率で辿るとすると、遷移確率は：
$$ p_{ij} = \frac{1}{L_i} \quad \text{（ページ $i$ からページ $j$ へのリンクが存在する場合）} $$
ダンピングファクター
実際のWebにはリンクが出ていないページ（dangling node）や孤立したクラスタがあり、既約性が保証されません。そこで、確率 $\alpha$ でリンクを辿り、確率 $1-\alpha$ でランダムなページに飛ぶモデルを使います。
$$ \bm{M} = \alpha \bm{P} + (1-\alpha)\frac{1}{N}\bm{J} $$
ここで $\bm{J}$ は全要素が1の $N \times N$ 行列です。$\alpha$ はダンピングファクター（通常 $\alpha = 0.85$）。
$\bm{M}$ の全要素が正になるため、ペロン-フロベニウスの定理から一意の定常分布が存在し、これがPageRankベクトルです。
Pythonでの実装
import numpy as np
import matplotlib.pyplot as plt
np.random.seed(42)
fig, axes = plt.subplots(2, 2, figsize=(14, 10))
# --- 遷移確率行列（天気モデル） ---
P = np.array([
[0.7, 0.2, 0.1],
[0.3, 0.4, 0.3],
[0.2, 0.3, 0.5]
])
states = ['晴れ', '曇り', '雨']
# --- (1) 方法1: 連立方程式で定常分布を求める ---
# πP = π, Σπ = 1
# (P^T - I)π^T = 0 に正規化条件を追加
A = (P.T - np.eye(3))
# 最後の方程式を正規化条件に置換
A[-1, :] = 1
b = np.zeros(3)
b[-1] = 1
pi_eq = np.linalg.solve(A, b)
print("方法1（連立方程式）: π =", pi_eq.round(6))
# --- (2) 方法2: 固有値分解 ---
eigenvalues, eigenvectors = np.linalg.eig(P.T)
idx = np.argmin(np.abs(eigenvalues - 1.0))
pi_eig = np.real(eigenvectors[:, idx])
pi_eig = pi_eig / pi_eig.sum()
print("方法2（固有値分解）: π =", pi_eig.round(6))
# --- (3) 方法3: P^n の極限 ---
Pn = np.linalg.matrix_power(P, 100)
pi_power = Pn[0]
print("方法3（P^100の行） : π =", pi_power.round(6))
# --- 可視化 (1): 定常分布の棒グラフ ---
x = np.arange(3)
axes[0, 0].bar(x, pi_eq, color=['gold', 'gray', 'steelblue'], alpha=0.8, edgecolor='black')
for i, v in enumerate(pi_eq):
axes[0, 0].text(i, v + 0.01, f'{v:.4f}', ha='center', fontsize=11)
axes[0, 0].set_xticks(x)
axes[0, 0].set_xticklabels(states, fontsize=12)
axes[0, 0].set_ylabel('$\\pi_i$')
axes[0, 0].set_title('Stationary Distribution $\\bm{\\pi}P = \\bm{\\pi}$')
axes[0, 0].grid(True, alpha=0.3, axis='y')
axes[0, 0].set_ylim(0, 0.6)
# --- 可視化 (2): 初期分布からの収束 ---
initial_distributions = [
np.array([1, 0, 0]), # 晴れ確定
np.array([0, 1, 0]), # 曇り確定
np.array([0, 0, 1]), # 雨確定
]
init_labels = ['晴れから', '曇りから', '雨から']
linestyles = ['-', '--', ':']
n_max = 25
for init_dist, label, ls in zip(initial_distributions, init_labels, linestyles):
pi_n = init_dist.copy().astype(float)
history = [pi_n.copy()]
for _ in range(n_max):
pi_n = pi_n @ P
history.append(pi_n.copy())
history = np.array(history)
# 晴れの確率のみプロット
axes[0, 1].plot(range(n_max + 1), history[:, 0], 'r' + ls, linewidth=2,
label=f'{label}: $\\pi_{{晴れ}}$')
axes[0, 1].axhline(y=pi_eq[0], color='k', linestyle='--', alpha=0.5, label=f'定常値 = {pi_eq[0]:.4f}')
axes[0, 1].set_xlabel('Step $n$')
axes[0, 1].set_ylabel('$\\pi_{晴れ}^{(n)}$')
axes[0, 1].set_title('Convergence to Stationary Distribution')
axes[0, 1].legend(fontsize=8)
axes[0, 1].grid(True, alpha=0.3)
# --- 可視化 (3): 固有値スペクトル ---
eigenvalues_all = np.linalg.eigvals(P)
theta = np.linspace(0, 2*np.pi, 100)
axes[1, 0].plot(np.cos(theta), np.sin(theta), 'k--', alpha=0.3, label='単位円')
for i, ev in enumerate(eigenvalues_all):
axes[1, 0].plot(np.real(ev), np.imag(ev), 'ro', markersize=10, zorder=5)
axes[1, 0].annotate(f'$\\lambda_{i+1}={ev:.4f}$',
(np.real(ev), np.imag(ev)),
textcoords="offset points", xytext=(10, 10), fontsize=10)
axes[1, 0].set_xlabel('Real')
axes[1, 0].set_ylabel('Imaginary')
axes[1, 0].set_title('Eigenvalues of $P$ (all $|\\lambda| \\leq 1$)')
axes[1, 0].set_aspect('equal')
axes[1, 0].grid(True, alpha=0.3)
axes[1, 0].legend()
axes[1, 0].set_xlim(-1.5, 1.5)
axes[1, 0].set_ylim(-1.5, 1.5)
# --- 可視化 (4): PageRank の例 ---
# 4ページのWebグラフ
# Page 0 → 1, 2
# Page 1 → 2
# Page 2 → 0
# Page 3 → 0, 2
n_pages = 4
links = {0: [1, 2], 1: [2], 2: [0], 3: [0, 2]}
# 遷移確率行列
P_web = np.zeros((n_pages, n_pages))
for i, targets in links.items():
for j in targets:
P_web[i, j] = 1.0 / len(targets)
# ダンピングファクター
alpha = 0.85
M = alpha * P_web + (1 - alpha) / n_pages * np.ones((n_pages, n_pages))
# 定常分布（PageRank）を固有値分解で計算
eigenvalues_pr, eigenvectors_pr = np.linalg.eig(M.T)
idx_pr = np.argmin(np.abs(eigenvalues_pr - 1.0))
pagerank = np.real(eigenvectors_pr[:, idx_pr])
pagerank = pagerank / pagerank.sum()
# べき乗法でも計算
pi_pr = np.ones(n_pages) / n_pages
for _ in range(100):
pi_pr = pi_pr @ M
print(f"\nPageRank（固有値分解）: {pagerank.round(4)}")
print(f"PageRank（べき乗法） : {pi_pr.round(4)}")
# 棒グラフ
page_labels = [f'Page {i}' for i in range(n_pages)]
colors_pr = plt.cm.Set2(np.linspace(0, 1, n_pages))
bars = axes[1, 1].bar(range(n_pages), pagerank, color=colors_pr, edgecolor='black', alpha=0.8)
for i, v in enumerate(pagerank):
axes[1, 1].text(i, v + 0.005, f'{v:.4f}', ha='center', fontsize=10)
axes[1, 1].set_xticks(range(n_pages))
axes[1, 1].set_xticklabels(page_labels)
axes[1, 1].set_ylabel('PageRank')
axes[1, 1].set_title(f'PageRank ($\\alpha = {alpha}$)')
axes[1, 1].grid(True, alpha=0.3, axis='y')
plt.tight_layout()
plt.savefig('stationary_distribution.png', dpi=150, bbox_inches='tight')
plt.show()
# 詳細釣り合い条件の検証
print("\n=== 詳細釣り合い条件の検証（天気モデル） ===")
is_reversible = True
for i in range(3):
for j in range(i+1, 3):
lhs = pi_eq[i] * P[i, j]
rhs = pi_eq[j] * P[j, i]
match = np.isclose(lhs, rhs)
if not match:
is_reversible = False
print(f"π({states[i]})×P({states[i]}→{states[j]}) = {lhs:.6f}, "
f"π({states[j]})×P({states[j]}→{states[i]}) = {rhs:.6f}, "
f"{'✓' if match else '✗'}")
print(f"このマルコフ連鎖は{'可逆' if is_reversible else '非可逆'}です")
# 収束速度
lambda_2 = sorted(np.abs(eigenvalues_all))[-2]
print(f"\n=== 収束速度 ===")
print(f"第2固有値の絶対値: |λ₂| = {lambda_2:.4f}")
print(f"収束速度 ∝ |λ₂|^n（小さいほど速い）")
print(f"混合時間の目安: 1/|log|λ₂|| ≈ {-1/np.log(lambda_2):.1f} ステップ")
まとめ
本記事では、マルコフ連鎖の定常分布について解説しました。
- 定常分布: $\bm{\pi}\bm{P} = \bm{\pi}$, $\sum_i \pi_i = 1$ を満たす確率ベクトル
- 存在条件: 有限状態の既約マルコフ連鎖には一意の定常分布が存在
- エルゴード定理: 既約かつ非周期なら、$\lim_{n \to \infty} p_{ij}^{(n)} = \pi_j$（初期分布によらず収束）
- 計算法: 連立方程式、固有値分解、$\bm{P}^n$ の極限、べき乗法
- 詳細釣り合い: $\pi_i p_{ij} = \pi_j p_{ji}$ は定常分布の十分条件（可逆マルコフ連鎖）
- PageRank: ダンピングファクター付きのマルコフ連鎖の定常分布
次のステップとして、以下の記事も参考にしてください。