サイコロを振り続けるとき、次の目は過去の履歴と無関係に決まります。一方、人間の感情や天気は「直前の状態」に強く引きずられるけれども、もっと昔のことは案外関係なかったりします。「未来は現在だけに依存する」というこの単純化を確率モデルに落とし込んだものがマルコフ連鎖です。20世紀初頭にロシアの数学者 A. A. Markov が「プーシキンの詩のなかの子音と母音の並び」を解析するために導入したこの仕組みは、いまや確率論・統計学・機械学習・物理学・経済学・情報理論まで、信じられないほど広い分野の共通言語になっています。
身近な応用例だけでも枚挙にいとまがありません。Google を世界一の検索エンジンに押し上げた PageRank はWebリンクから作る巨大なマルコフ連鎖の定常分布計算ですし、ベイズ統計でほぼ唯一の高次元事後分布計算の現実解である MCMC（マルコフ連鎖モンテカルロ） は、サンプルそのものをマルコフ連鎖として作ります。音声認識や形態素解析のバックボーンである HMM（隠れマルコフモデル）、強化学習の理論的支柱となる マルコフ決定過程（MDP）、自然言語処理における N-gram言語モデル、さらに統計物理の Ising モデルや材料科学のキネティックモンテカルロまで、いずれもマルコフ連鎖が心臓部です。
ところが、初学者向けの解説は「定義式と天気の例」で終わってしまうことが多く、いざ応用に進もうとすると「定常分布はなぜ存在するのか？」「いつ収束するのか？」「Chapman–Kolmogorov は何が嬉しいのか？」といった本質的な疑問でつまずきがちです。本記事ではマルコフ連鎖を応用に耐える形で理解できるよう、定義・代数的構造・解析的性質・極限定理・実装・応用を一気通貫で扱います。
本記事の内容
- マルコフ性の厳密な定義と直感的意味
- 遷移行列・初期分布・状態分布の代数
- Chapman–Kolmogorov 等式と $n$ ステップ遷移
- 既約性・再帰性・周期性による状態の分類
- 定常分布の存在・一意性・収束定理
- エルゴード定理（時間平均＝空間平均）
- Pythonでのシミュレーション・固有値計算・収束可視化
- 応用: PageRank, MCMC, HMM, マルコフ連鎖言語モデル
前提知識
この記事を読む前に、以下の概念を理解しておくとスムーズです。
- 条件付き確率と全確率の公式
- 確率変数と確率分布
- 行列の積・固有値・固有ベクトル
- （任意）測度論的な確率の言葉づかい
マルコフ性とは何か — 直感と定義
直感的理解 — 「過去は現在に集約されている」
シャッフルされたトランプの山から1枚ずつめくっていく場面を考えてみてください。残りの山にどんな札が何枚あるか「現在の山の状態」を全部覚えていれば、次にめくられる札の確率分布は完全に決まります。逆に、最初から今までにめくられた順番を細かく記憶していたとしても、現在の山の状態が同じなら、得られる情報量は変わりません。これがマルコフ性の核心です。「未来予測に必要な情報は、現在の状態のなかにすべて畳み込まれている」。
別の例: 酔った人が街路をふらふら歩く「ランダムウォーク」。次にどっちに行くかは「今いる交差点」だけで決まり、どんな経路を経てそこに来たかは関係ありません。ボードゲームの双六、エレベーターが停まる階のシーケンス、株価レジームの切り替わり——いずれも「現在の状態 → 次の状態」の遷移確率さえ書ければモデル化できます。
厳密な定義
時間添字 $n = 0, 1, 2, \ldots$ をもつ確率変数の列 $\{X_n\}$ が、可算な状態空間 $\mathcal{S}$ に値を取るとします。$\{X_n\}$ がマルコフ性を持つとは、任意の $n$ と任意の状態列 $i_0, \ldots, i_{n-1}, i, j \in \mathcal{S}$ に対して
$$ P(X_{n+1} = j \mid X_n = i, X_{n-1} = i_{n-1}, \ldots, X_0 = i_0) = P(X_{n+1} = j \mid X_n = i) $$
が成り立つことです。左辺は「過去全履歴を条件にした次状態の確率」、右辺は「現在状態だけを条件にした次状態の確率」。これらが等しいというのが、マルコフ性の主張です。
さらに右辺が $n$ に依存しないとき、すなわち
$$ P(X_{n+1} = j \mid X_n = i) = p_{ij} \quad (\text{$n$ によらず}) $$
を満たすとき、$\{X_n\}$ を時間的同次（time-homogeneous）なマルコフ連鎖と呼びます。以下、本記事では特に断らない限り時間的同次を仮定します。
マルコフ性が「強い仮定」であることを忘れない
マルコフ性は便利ですが、現実のデータが必ずしも満たすわけではありません。たとえば「人の食欲」は「直近3日間の食事内容」に依存するかもしれず、その場合は1次マルコフでは不足です。しかし、状態空間を拡張して「過去3日分の食事の組」を1つの状態とみなせば、再びマルコフ性を回復できます。これを「状態拡張によるマルコフ化」と呼び、HMMやLSTMの設計思想にも通じます。
別の見方をすると、マルコフ性は確率モデルの記憶長を1ステップに制限する仮定です。記憶長を増やしたければ状態空間を大きくする——このトレードオフは現代の時系列モデリングでも本質的で、Transformer の注意機構が「直近 $k$ トークンの組み合わせ」を状態として扱っていると見ることもできます。シンプルな構造を維持しながら表現力をどこまで広げられるか、というのが応用上のキモです。
ここまでで「現在の状態を起点に次の状態が決まる」という抽象的な構造を把握しました。この構造を具体的に扱うためには、遷移確率を行列で表現するのが圧倒的に便利です。
遷移行列 — マルコフ連鎖の代数的表現
遷移確率の定義
状態 $i$ から状態 $j$ への1ステップ遷移確率を
$$ p_{ij} = P(X_{n+1} = j \mid X_n = i) $$
と書きます。$|\mathcal{S}| = N$（有限状態）と仮定し、$i, j \in \{1, 2, \ldots, N\}$ とすると、遷移確率は $N \times N$ の行列にまとめられます。
$$ \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} $$
これを遷移行列（transition matrix）と呼びます。$\bm{P}$ は次の2つの性質を満たします。
- 非負性: すべての $i, j$ で $p_{ij} \geq 0$
- 行確率性: 各行の和が $\sum_{j=1}^{N} p_{ij} = 1$
このような行列を （行）確率行列 または stochastic matrix と呼びます。直感的には、$i$ 行目は「現在状態が $i$ のときの次状態の確率分布」を表します。だから行の和が1になるのです。
状態分布の時間発展
時刻 $n$ における状態の確率分布を、行ベクトル $\bm{\pi}^{(n)} = (\pi_1^{(n)}, \ldots, \pi_N^{(n)})$ で表します。ここで $\pi_i^{(n)} = P(X_n = i)$ です。マルコフ性と全確率の公式から、
$$ \pi_j^{(n+1)} = \sum_{i=1}^{N} P(X_{n+1} = j \mid X_n = i) P(X_n = i) = \sum_{i=1}^{N} p_{ij} \pi_i^{(n)} $$
が得られます。行列とベクトルで書けば
$$ \bm{\pi}^{(n+1)} = \bm{\pi}^{(n)} \bm{P} $$
です。マルコフ連鎖の時間発展は、行ベクトルへの行列の右乗算として表現できる——これがマルコフ連鎖を線形代数の道具で解析できる理由です。初期分布 $\bm{\pi}^{(0)}$ さえ与えれば、
$$ \bm{\pi}^{(n)} = \bm{\pi}^{(0)} \bm{P}^n $$
で任意時刻の分布が計算できます。
具体例: 天気モデル
「晴れ」と「雨」の2状態モデルを取り上げます。
$$ \bm{P} = \begin{pmatrix} 0.8 & 0.2 \\ 0.4 & 0.6 \end{pmatrix} $$
ここで第1行・第2行は「今日晴れの場合」「今日雨の場合」の翌日の確率を表します。各行の和が確かに 1 になっていることを確認してください。今日晴れる確率が 1（$\bm{\pi}^{(0)} = (1, 0)$）の場合、明日の分布は
$$ \bm{\pi}^{(1)} = (1, 0) \bm{P} = (0.8, 0.2) $$
となり、明日も晴れの確率は0.8、雨の確率は0.2です。明後日は
$$ \bm{\pi}^{(2)} = \bm{\pi}^{(1)} \bm{P} = (0.8, 0.2) \begin{pmatrix} 0.8 & 0.2 \\ 0.4 & 0.6 \end{pmatrix} = (0.72, 0.28) $$
行列のべき乗を計算する方が早いので $\bm{P}^2$ を直接計算してもよく、
$$ \bm{P}^2 = \begin{pmatrix} 0.8 & 0.2 \\ 0.4 & 0.6 \end{pmatrix} \begin{pmatrix} 0.8 & 0.2 \\ 0.4 & 0.6 \end{pmatrix} = \begin{pmatrix} 0.72 & 0.28 \\ 0.56 & 0.44 \end{pmatrix} $$
$\bm{P}^2$ の $(i, j)$ 成分は「$i$ から出発して2ステップ後に $j$ にいる確率」になっています。
次に、この「$n$ ステップ遷移」を一般化する Chapman–Kolmogorov の等式に進みます。
Chapman–Kolmogorov 等式 — $n$ ステップ遷移の代数
$n$ ステップ遷移確率
状態 $i$ から $n$ ステップ後に状態 $j$ にいる確率を
$$ p_{ij}^{(n)} = P(X_{m+n} = j \mid X_m = i) $$
と定義します（時間的同次なので $m$ にはよりません）。これらを並べた行列を $\bm{P}^{(n)}$ と書きます。
Chapman–Kolmogorov 等式の主張
任意の $0 \leq k \leq n$ に対し
$$ p_{ij}^{(n)} = \sum_{l \in \mathcal{S}} p_{il}^{(k)} p_{lj}^{(n-k)} $$
が成り立ち、行列でまとめると
$$ \bm{P}^{(n)} = \bm{P}^{(k)} \bm{P}^{(n-k)} = \bm{P}^n $$
です。証明は全確率の公式とマルコフ性を機械的に使うだけです。
導出
$X_m = i$ を所与にして $X_{m+n} = j$ の確率を計算します。途中時刻 $X_{m+k}$ が取る値 $l$ を場合分けすることで、
$$ P(X_{m+n} = j \mid X_m = i) = \sum_l P(X_{m+n} = j, X_{m+k} = l \mid X_m = i) $$
ここで右辺の各項に乗法定理を適用します。$P(A, B \mid C) = P(A \mid B, C) P(B \mid C)$ より
$$ P(X_{m+n} = j, X_{m+k} = l \mid X_m = i) = P(X_{m+n} = j \mid X_{m+k} = l, X_m = i) \cdot P(X_{m+k} = l \mid X_m = i) $$
マルコフ性により、第1項の条件から $X_m = i$ は落とせます（過去の情報は現在に集約済み）。すると
$$ P(X_{m+n} = j \mid X_{m+k} = l, X_m = i) = P(X_{m+n} = j \mid X_{m+k} = l) = p_{lj}^{(n-k)} $$
であり、第2項は $p_{il}^{(k)}$ そのものです。合わせて
$$ p_{ij}^{(n)} = \sum_l p_{il}^{(k)} p_{lj}^{(n-k)} $$
を得ます。これは行列積の定義そのものですから、$\bm{P}^{(n)} = \bm{P}^{(k)} \bm{P}^{(n-k)}$ という行列恒等式に他なりません。
何が嬉しいのか
Chapman–Kolmogorov を一見「自明な恒等式」と感じるかもしれませんが、これは強力な計算手段です。$\bm{P}^{1024}$ を求めたいとき、$\bm{P}$ を1024回掛けるのではなく、二乗を10回繰り返せば $\bm{P}^2, \bm{P}^4, \bm{P}^8, \ldots, \bm{P}^{1024}$ が得られます。これは O(log n) 計算量で、シミュレーションを介さずに正確な確率分布を求める王道です。
$\bm{P}^n$ の各行が $n$ が大きいときに「同じ分布」に収束する様子が見えてきます。この「同じ分布」こそが定常分布ですが、その存在と一意性を保証するには、まず状態の構造を分類しておく必要があります。
状態の分類 — 既約性・再帰性・周期性
到達可能性と連絡
状態 $j$ が状態 $i$ から到達可能（accessible）であるとは、ある $n \geq 0$ が存在して $p_{ij}^{(n)} > 0$ となること、つまり「$i$ から出発して $j$ にたどり着く可能性がある」ことを言います。$i \to j$ と書きます。
$i \to j$ かつ $j \to i$ のとき、$i$ と $j$ は連絡する（communicate）と言い、$i \leftrightarrow j$ と書きます。連絡関係は同値関係になり、状態空間を連絡類（communicating class）に分割します。
既約性
すべての状態ペアが互いに連絡するとき、マルコフ連鎖は既約（irreducible）であると言います。既約とは要するに「どの状態からでも、どの状態にも有限ステップで到達できる」ことです。グラフ理論の言葉で言えば、有向グラフが強連結であることに対応します。
周期性
状態 $i$ の周期 $d(i)$ を
$$ d(i) = \gcd\{n \geq 1 : p_{ii}^{(n)} > 0\} $$
で定義します。$d(i) = 1$ のとき $i$ は非周期的（aperiodic）であると言い、$d(i) > 1$ のときは「$i$ に戻るには必ず $d(i)$ の倍数ステップが必要」という強い構造があります。同じ連絡類に属する状態は同じ周期を持つので、既約な連鎖については単に「連鎖の周期」と呼べます。
再帰性と過渡性
状態 $i$ から出発して、いつかは $i$ に戻る確率を $f_i$ と書きます。
- $f_i = 1$ のとき $i$ を再帰的（recurrent）と呼ぶ
- $f_i < 1$ のとき $i$ を過渡的（transient）と呼ぶ
再帰的な状態にはさらに区別があります。「戻ってくるまでの期待時間」を $E[T_i]$ とすると
- $E[T_i] < \infty$ のとき正再帰（positive recurrent）
- $E[T_i] = \infty$ のとき零再帰（null recurrent）
有限状態空間で既約な連鎖はすべての状態が自動的に正再帰になるので、本記事の実用例ではほぼ「正再帰」と思って差し支えありません。零再帰は無限状態（典型的には2次元以上のランダムウォーク）に固有の現象です。
吸収状態
特殊な再帰状態として、$p_{ii} = 1$ となる状態を吸収状態（absorbing state）と呼びます。一度入ると出られない「終着駅」です。たとえばギャンブラーが破産する状態、SIRモデルの回復済み状態、デバイスの故障状態などが該当します。
ここまでで「どの状態がどのような構造で結ばれているか」が見えるようになりました。次は、十分に時間が経ったときに連鎖が落ち着く先——定常分布——を考えます。
定常分布と収束定理
定常分布の定義
確率ベクトル $\bm{\pi} = (\pi_1, \ldots, \pi_N)$（$\pi_i \geq 0$, $\sum_i \pi_i = 1$）が定常分布（stationary distribution）であるとは
$$ \bm{\pi} \bm{P} = \bm{\pi} $$
すなわち、行ベクトル $\bm{\pi}$ が遷移行列 $\bm{P}$ の固有値1に対応する左固有ベクトルであり、かつ確率分布として正規化されていることを言います。意味合いとしては「$\bm{\pi}$ に従って状態を選んだ後で1ステップ遷移しても、再び $\bm{\pi}$ に従う」というバランス条件です。
存在と一意性 — 有限状態のケース
有限状態空間 $|\mathcal{S}| < \infty$ の場合、既約で非周期的なマルコフ連鎖（= エルゴード的）には、定常分布 $\bm{\pi}$ がただ一つ存在し、しかも任意の初期分布 $\bm{\pi}^{(0)}$ に対して
$$ \bm{\pi}^{(0)} \bm{P}^n \xrightarrow[n \to \infty]{} \bm{\pi} $$
となります。これがマルコフ連鎖の極限定理（fundamental limit theorem）です。
直感的には次のように理解できます。既約性は「どこから出ても全状態を巡れる」ことを、非周期性は「巡る歩数のパターンに最大公約数のずれが生じない」ことを保証し、両方そろうことで初期条件への依存が時間とともに完全に拭われるのです。
スペクトル分解の視点
確率行列 $\bm{P}$ の固有値を $1 = \lambda_1 > |\lambda_2| \geq \cdots \geq |\lambda_N|$ とソートします（既約・非周期なら $|\lambda_2| < 1$ が成り立ちます）。各固有値に対応する右固有ベクトル $\bm{r}_k$ と左固有ベクトル $\bm{l}_k$ を取れば、
$$ \bm{P}^n = \sum_{k=1}^{N} \lambda_k^n \, \bm{r}_k \bm{l}_k^\top $$
と分解できます。$n \to \infty$ で $\lambda_k^n \to 0$（$k \geq 2$）となるため、定常成分 $\bm{r}_1 \bm{l}_1^\top$（= 全行が $\bm{\pi}$ になる行列）だけが残り、まさに収束が起こります。
第2固有値の大きさ $|\lambda_2|$ は収束の速さを支配する量で、誤差は $O(|\lambda_2|^n)$ で減衰します。スペクトルギャップ $1 – |\lambda_2|$ が大きいほど混合（mixing）が速く、$|\lambda_2|$ が1に近いと「ほぼ吸収的」な状況になり収束が極めて遅い、というイメージです。MCMCの設計でも「いかに $|\lambda_2|$ を小さくするか」が実用的な品質を決めるポイントになります。
計算: 天気モデルの定常分布
$\bm{\pi} \bm{P} = \bm{\pi}$ かつ $\pi_S + \pi_R = 1$ を解きましょう。具体的に書き下すと
$$ \begin{cases} 0.8 \pi_S + 0.4 \pi_R = \pi_S \\ 0.2 \pi_S + 0.6 \pi_R = \pi_R \\ \pi_S + \pi_R = 1 \end{cases} $$
第1式を整理すると $0.4 \pi_R = 0.2 \pi_S$、つまり $\pi_S = 2 \pi_R$ です。これを $\pi_S + \pi_R = 1$ に代入して $2\pi_R + \pi_R = 1$、よって
$$ \pi_S = \frac{2}{3}, \quad \pi_R = \frac{1}{3} $$
長期的には「3日のうち2日は晴れ、1日は雨」になるわけです。
詳細つり合いと可逆連鎖
定常分布の十分条件として知られるのが詳細つり合い（detailed balance）です。
$$ \pi_i p_{ij} = \pi_j p_{ji} \quad (\forall i, j) $$
これを満たす $\bm{\pi}$ は定常分布になります。なぜか？両辺を $i$ について和を取ると
$$ \sum_i \pi_i p_{ij} = \sum_i \pi_j p_{ji} = \pi_j \sum_i p_{ji} = \pi_j $$
となり、確かに $\bm{\pi}\bm{P} = \bm{\pi}$ を満たします。詳細つり合いを満たすマルコフ連鎖は可逆（reversible）と呼ばれ、物理学・統計学で重要な役割を果たします。MCMCのMetropolis–Hastings法は、詳細つり合いを自動的に成り立たせるよう提案受理を設計することで目標分布を定常分布にしています。
ここまでで「分布」の収束は理解できました。次に問いたいのは「ひとつの軌跡（サンプル経路）の時間平均」はどうなるか？という、より強い主張——エルゴード定理です。
エルゴード定理 — 時間平均は空間平均に等しい
主張
エルゴード的マルコフ連鎖において、任意の有界関数 $f: \mathcal{S} \to \mathbb{R}$ に対し、サンプル経路 $X_0, X_1, X_2, \ldots$ の時間平均は、定常分布のもとでの期待値（空間平均）に確率1で収束します。
$$ \frac{1}{n} \sum_{t=0}^{n-1} f(X_t) \xrightarrow[n \to \infty]{\text{a.s.}} \sum_{i \in \mathcal{S}} \pi_i f(i) = \mathbb{E}_{\bm{\pi}}[f(X)] $$
ここで a.s. は almost surely（ほとんど確実に）の略です。「1本の長い軌跡を観測すれば、全状態にわたる空間平均が分かる」という主張で、これがマルコフ連鎖モンテカルロ（MCMC）の理論的根拠です。
特別な場合: 再帰時間との関係
$f(i) = \mathbb{1}[i = j]$（状態 $j$ にいるかどうかの指示関数）と取ると、時間平均は「状態 $j$ にいた割合」、空間平均は $\pi_j$ になります。したがってエルゴード定理は
$$ \frac{1}{n} \sum_{t=0}^{n-1} \mathbb{1}[X_t = j] \xrightarrow{\text{a.s.}} \pi_j $$
を意味します。さらにこれを「初めて $j$ に戻るまでの時間」の期待値 $E[T_j]$ と結びつけると
$$ \pi_j = \frac{1}{E[T_j]} $$
という美しい関係式が得られます。定常分布は平均回帰時間の逆数——これは長期挙動を直感的に理解する強力な公式です。
このエルゴード性は、後に出てくる MCMC で「サンプルが目標分布から取れていることが、どう保証されるのか？」の答えにもなっています。
Pythonでの実装と可視化
ここからは理論を Python コードで体感していきましょう。まずは天気モデルの軌跡シミュレーションです。
import numpy as np
import matplotlib.pyplot as plt
def simulate_markov_chain(P, initial_state, n_steps, rng=None):
"""マルコフ連鎖の軌跡を1本シミュレーションする"""
if rng is None:
rng = np.random.default_rng()
n_states = P.shape[0]
states = np.empty(n_steps + 1, dtype=int)
states[0] = initial_state
for t in range(n_steps):
states[t + 1] = rng.choice(n_states, p=P[states[t]])
return states
# 天気モデル
P = np.array([[0.8, 0.2],
[0.4, 0.6]])
state_names = ['Sunny', 'Rainy']
rng = np.random.default_rng(42)
n_steps = 500
states = simulate_markov_chain(P, initial_state=0, n_steps=n_steps, rng=rng)
fig, axes = plt.subplots(2, 1, figsize=(12, 6))
# 上段: 軌跡
ax = axes[0]
ax.step(range(n_steps + 1), states, where='post', linewidth=0.9)
ax.set_yticks([0, 1])
ax.set_yticklabels(state_names)
ax.set_xlabel('time step n')
ax.set_title('Markov chain trajectory (Weather model)')
ax.grid(alpha=0.3)
# 下段: Sunny の経験頻度の収束
ax = axes[1]
sunny_ratio = np.cumsum(states == 0) / np.arange(1, n_steps + 2)
ax.plot(range(n_steps + 1), sunny_ratio, linewidth=1.3)
ax.axhline(2/3, ls='--', label=r'$\pi_{\mathrm{Sunny}} = 2/3$')
ax.set_xlabel('time step n')
ax.set_ylabel('cumulative fraction of Sunny')
ax.set_title('Convergence to stationary distribution')
ax.legend()
ax.grid(alpha=0.3)
plt.tight_layout()
plt.show()
上段は実際の状態の遷移、下段は「時刻 $n$ までに Sunny だった割合」です。下段のカーブは初期は揺らぎが大きいものの、$n$ が大きくなると確実に理論値 $2/3 \approx 0.667$ に収束していきます。これはエルゴード定理の経験的検証になっています。
次に、$\bm{P}^n$ を直接計算して定常分布への収束を確認します。
import numpy as np
import matplotlib.pyplot as plt
P = np.array([[0.8, 0.2],
[0.4, 0.6]])
n_max = 30
Pn_history = np.zeros((n_max + 1, 2, 2))
Pn_history[0] = np.eye(2)
for n in range(1, n_max + 1):
Pn_history[n] = Pn_history[n - 1] @ P
fig, axes = plt.subplots(1, 2, figsize=(12, 4.5))
for row, title in enumerate(['Starting from Sunny', 'Starting from Rainy']):
ax = axes[row]
ax.plot(range(n_max + 1), Pn_history[:, row, 0], 'o-', ms=3, label=r'$P^n_{,\to S}$')
ax.plot(range(n_max + 1), Pn_history[:, row, 1], 's-', ms=3, label=r'$P^n_{,\to R}$')
ax.axhline(2/3, ls='--', alpha=0.5)
ax.axhline(1/3, ls='--', alpha=0.5)
ax.set_xlabel('n')
ax.set_ylabel('probability')
ax.set_title(title)
ax.legend()
ax.grid(alpha=0.3)
plt.tight_layout()
plt.show()
初期状態が Sunny でも Rainy でも、$n \to \infty$ で同じ定常分布 $(2/3, 1/3)$ に収束していくのが見えます。これは「定常分布の一意性 + 任意初期分布からの収束」というエルゴード性の表明そのものです。
定常分布を固有値問題として解く
ループを回さずに、線形代数の道具一発でも定常分布は求められます。
import numpy as np
P = np.array([[0.8, 0.2],
[0.4, 0.6]])
# 定常分布は P^T の固有値1に対応する左固有ベクトル
# (πP = π の両辺を転置すると P^T π^T = π^T)
eigvals, eigvecs = np.linalg.eig(P.T)
idx = np.argmin(np.abs(eigvals - 1.0))
pi = np.real(eigvecs[:, idx])
pi = pi / pi.sum()
print(f"固有値: {eigvals}")
print(f"定常分布 π: {pi}")
print(f"検証 πP - π: {pi @ P - pi}")
固有値の片方が厳密に 1、もう一方が $0.4$ となります。理論上、確率行列は必ず固有値1を持ち、他の固有値は $|\lambda| \leq 1$ を満たします。第2固有値の絶対値 $|\lambda_2|$ は収束の速さを決める量で、$1 – |\lambda_2|$ をスペクトルギャップと呼びます。
3状態への拡張と Chapman–Kolmogorov の検証
import numpy as np
P = np.array([
[0.6, 0.3, 0.1],
[0.2, 0.5, 0.3],
[0.1, 0.2, 0.7],
])
# 行列のべき乗で n ステップ遷移
P_5 = np.linalg.matrix_power(P, 5)
P_20 = np.linalg.matrix_power(P, 20)
print("P^5 =")
print(P_5)
print("\nP^20 =")
print(P_20)
# Chapman-Kolmogorov: P^20 = P^7 * P^13
P_7 = np.linalg.matrix_power(P, 7)
P_13 = np.linalg.matrix_power(P, 13)
diff = np.max(np.abs(P_20 - P_7 @ P_13))
print(f"\n|P^20 - P^7 @ P^13|_max = {diff:.2e}")
P_20 の各行はほぼ同じ値の3つ組になっており、定常分布に収束していることが見て取れます。また Chapman–Kolmogorov の検証は $10^{-15}$ オーダーの誤差で一致——浮動小数点演算の丸め誤差レベルです。
エルゴード性: 時間平均と空間平均の一致
import numpy as np
P = np.array([
[0.7, 0.2, 0.1],
[0.2, 0.5, 0.3],
[0.1, 0.3, 0.6],
])
# 真の定常分布
eigvals, eigvecs = np.linalg.eig(P.T)
pi = np.real(eigvecs[:, np.argmin(np.abs(eigvals - 1.0))])
pi /= pi.sum()
rng = np.random.default_rng(13)
N = 20000
traj = np.zeros(N, dtype=int)
for t in range(1, N):
traj[t] = rng.choice(3, p=P[traj[t - 1]])
time_avg = np.array([(traj == s).mean() for s in range(3)])
print("真の定常分布 π:", pi)
print("時間平均（軌跡1本）:", time_avg)
print("最大誤差:", np.max(np.abs(time_avg - pi)))
1本の長い軌跡から得た「各状態にいた割合」は、理論定常分布と $10^{-3}$ 精度で一致します。これがエルゴード定理の力です。多数本のサンプルを集める代わりに、1本の長いシミュレーションで分布を推定できる——これは MCMC の基本原理そのものです。
理論と実装で基礎を固めたので、ここからは応用編に進みます。マルコフ連鎖の威力を実感できる、4つの代表的な応用を駆け足で見ていきましょう。
応用1: PageRank — Webのマルコフ連鎖
Google を世界最大の企業に押し上げた PageRank アルゴリズムは、Webリンク構造から作るマルコフ連鎖の定常分布計算に過ぎません。
モデル化
Webページを状態、ハイパーリンクを遷移とみなします。ページ $i$ から出ているリンクが $d_i$ 個あり、それぞれ等確率でクリックすると仮定すると、遷移行列 $M$ は
$$ M_{ij} = \begin{cases} 1/d_i & \text{$i$ から $j$ へのリンクがある} \\ 0 & \text{なし} \end{cases} $$
となります。ただしリンクのないページ（dangling node）の扱いと、サーファーがリンクをたどらず突然別のページにジャンプする可能性を入れるため、ダンピング係数 $\alpha$（典型値 0.85）を導入して
$$ G = \alpha M + (1 – \alpha) \frac{1}{N} \bm{1} \bm{1}^\top $$
とします。この $G$ を遷移行列とするマルコフ連鎖の定常分布が PageRank です。
なぜマルコフ連鎖なのか
「無限ランダムサーファー」を考えます。確率 $\alpha$ でリンクをたどり、確率 $1-\alpha$ で全ページから一様にランダムジャンプ。十分長く歩き回ったあとサーファーがいるページの分布が PageRank です。エルゴード定理が成り立つように $G$ を構成しているため、Web全体が強連結でなくても収束が保証されます。
import numpy as np
# 5ページのトイ Web (A,B,C,D,E)
nodes = list("ABCDE")
edges = [("A", "B"), ("A", "C"), ("B", "C"), ("B", "D"),
("C", "A"), ("D", "A"), ("D", "E"), ("E", "A"), ("E", "B")]
N = len(nodes); idx = {n: i for i, n in enumerate(nodes)}
M = np.zeros((N, N))
for src, dst in edges:
M[idx[src], idx[dst]] += 1
M = M / M.sum(axis=1, keepdims=True)
alpha = 0.85
G = alpha * M + (1 - alpha) / N * np.ones((N, N))
# 主固有ベクトル (左) = PageRank
eigvals, eigvecs = np.linalg.eig(G.T)
pr = np.real(eigvecs[:, np.argmin(np.abs(eigvals - 1.0))])
pr /= pr.sum()
for n, v in zip(nodes, pr):
print(f"{n}: PR={v:.4f}")
このように、リンクを多く受けるページほど高い PR を得る一方、「多くのリンクを集めているページからリンクされている」ページも間接的に評価されます——再帰的な重要度測定がマルコフ連鎖の定常分布として自然に表現されているのです。
応用2: MCMC — 目標分布を定常分布に持つ連鎖
ベイズ統計で事後分布 $p(\theta \mid \mathcal{D}) \propto p(\mathcal{D} \mid \theta) p(\theta)$ を求めたいとき、高次元では解析的に積分できません。そこで「サンプルだけ」あれば期待値・分位点・周辺分布が推定できることを利用し、目標分布 $\pi(x)$ からサンプリングするマルコフ連鎖を逆算で設計するのが MCMC のアイデアです。
Metropolis–Hastings 法
提案分布 $q(x’ \mid x)$ で次候補を出し、受理確率
$$ \alpha(x, x’) = \min\left\{1, \frac{\pi(x’) q(x \mid x’)}{\pi(x) q(x’ \mid x)}\right\} $$
で受理する。提案が対称 ($q(x’|x) = q(x|x’)$) なら単に $\min\{1, \pi(x’)/\pi(x)\}$。この設計は詳細つり合い
$$ \pi(x) p(x \to x’) = \pi(x’) p(x’ \to x) $$
を自動的に成立させ、したがって $\pi$ が定常分布になります。エルゴード性が成り立てば、初期値に関係なく $\pi$ からのサンプルが得られるわけです。
import numpy as np
# 目標分布: 2成分混合正規（解析的に不規格でも比だけ分かればOK）
def target_logp(x):
p1 = np.exp(-0.5 * ((x - (-3)) / 1.0) ** 2) / 1.0
p2 = np.exp(-0.5 * ((x - (+3)) / 1.5) ** 2) / 1.5
return np.log(0.4 * p1 + 0.6 * p2 + 1e-300)
rng = np.random.default_rng(7)
n_iter = 5000
sigma_prop = 1.5
samples = np.zeros(n_iter); samples[0] = -8.0
acc = 0
for i in range(1, n_iter):
cur = samples[i - 1]
prop = cur + rng.normal(0, sigma_prop)
log_alpha = target_logp(prop) - target_logp(cur)
if np.log(rng.uniform()) < log_alpha:
samples[i] = prop; acc += 1
else:
samples[i] = cur
print(f"受理率: {acc / n_iter:.3f}")
print(f"サンプル平均: {samples[500:].mean():.3f}")
print(f"サンプル分散: {samples[500:].var():.3f}")
burn-in（最初の数百サンプル）を捨ててから集計すれば、2モード混合正規の平均（$0.4 \cdot (-3) + 0.6 \cdot 3 = 0.6$）に近い値が得られます。マルコフ連鎖が目標分布をきちんと再現していることが見て取れます。
応用3: HMM — マルコフ連鎖 + 観測モデル
隠れマルコフモデル（Hidden Markov Model） は、内部にマルコフ連鎖を持ちながら、状態そのものは直接観測できず、各時刻に「状態に依存する観測 $Y_t$」だけが得られるモデルです。
- 隠れ状態 $Z_t$ が遷移行列 $\bm{A}$ に従うマルコフ連鎖
- 観測 $Y_t$ が出力分布 $B_{Z_t}(\cdot)$ から生成される
音声認識（音素が隠れ状態、音響特徴が観測）、品詞タグ付け（品詞が隠れ状態、単語が観測）、株価レジーム検出（強気/弱気が隠れ状態、リターンが観測）など、応用は無数にあります。
推定問題は3つに整理されます。
- 評価: 観測列の尤度 $p(Y_{1:T} \mid \bm{A}, \bm{B}, \bm{\pi}_0)$ → Forward アルゴリズム
- 復号: 最尤状態系列 $\arg\max_{Z_{1:T}} p(Z_{1:T} \mid Y_{1:T})$ → Viterbi アルゴリズム
- 学習: パラメータ推定 $\hat{\bm{A}}, \hat{\bm{B}}$ → Baum–Welch（EM）アルゴリズム
これらは全て、内部のマルコフ連鎖がもつ「現在の状態だけで未来が決まる」性質を最大限に活用した動的計画法です。状態空間 $|\mathcal{S}| = K$、系列長 $T$ で $O(K^2 T)$ で計算できる——マルコフ性がもたらす計算上の恩恵がよくわかります。
隠れマルコフモデル(HMM)の理論と実装 で詳細を解説しています。
応用4: マルコフ連鎖言語モデル
文字（あるいは単語）を状態とみなし、コーパスから遷移確率 $p(c_{t+1} \mid c_t)$ を最尤推定したものが、最も素朴な統計的言語モデルです。Shannon が1948年の論文で英語の「2次マルコフ近似」「3次マルコフ近似」を作って、近似次数を上げるほど英語らしくなる様子を示したのは有名な話です。
import numpy as np
text = (
"the quick brown fox jumps over the lazy dog. "
"to be or not to be that is the question. " * 10
).lower()
chars = sorted(set(text))
n = len(chars); idx = {c: i for i, c in enumerate(chars)}
M = np.zeros((n, n))
for a, b in zip(text[:-1], text[1:]):
M[idx[a], idx[b]] += 1
M = M / np.maximum(M.sum(axis=1, keepdims=True), 1)
rng = np.random.default_rng(11)
cur = "t"; out = [cur]
for _ in range(200):
nxt_i = rng.choice(n, p=M[idx[cur]])
nxt = chars[nxt_i]; out.append(nxt); cur = nxt
print("".join(out))
生成結果は「英語っぽい音の並び」になりますが、文法構造は持ちません。これは1次マルコフ性の限界です。次数を上げる（N-gram）か、状態を埋め込みベクトルに置き換える（RNN/Transformer）ことで現代の言語モデルにつながります。マルコフ連鎖は、現代のLLMの最も基本的な前身だと言えます。
まとめ
本記事では、マルコフ連鎖の理論と応用を体系的に解説しました。
- マルコフ性: 未来は現在の状態のみに依存し、過去の履歴に直接は依存しない
- 遷移行列: 行確率性 $\sum_j p_{ij} = 1$ を満たす確率行列。状態分布の時間発展は $\bm{\pi}^{(n+1)} = \bm{\pi}^{(n)} \bm{P}$ で表される
- Chapman–Kolmogorov: $\bm{P}^{(n)} = \bm{P}^n$。$n$ ステップ遷移は単純な行列のべき乗で得られる
- 状態の分類: 既約性・周期性・再帰性が長期挙動を決める
- 定常分布: $\bm{\pi}\bm{P} = \bm{\pi}$ の解。エルゴード的（既約 + 非周期）なら一意に存在し、任意初期分布が収束する
- エルゴード定理: 時間平均が空間平均（$\mathbb{E}_{\bm{\pi}}[f]$）に収束。MCMCの理論的支柱
- 応用: PageRank（Web行列の主固有ベクトル）、MCMC（詳細つり合いで目標分布を定常化）、HMM（観測モデル付き）、N-gram言語モデル
マルコフ連鎖は確率過程の入り口でありながら、現代データ科学のあらゆる場面で姿を変えて登場します。
次のステップとして、以下の発展トピックを学んでみてください。