すごろくで、いまコマが 3 のマスにあるとします。「5 手後にコマが 8 のマスにいる確率」を知りたいとき、あなたはどうやって計算しますか。もっとも素朴なやり方は、5 手ぶんのサイコロの目の組み合わせをすべて数え上げることです。しかしマス目が多く手数が増えると、この組み合わせは天文学的な数になり、とても手では追えません。
ここで効いてくるのが「途中の 1 点で場合分けする」という発想です。5 手後にマス 8 にいるためには、たとえば 2 手後に必ずどこかのマスにいるはずです。そこで「2 手後にマス $k$ にいて、そこからさらに 3 手でマス 8 に着く」という確率を、すべての中継マス $k$ について足し合わせれば答えになります。この「途中で切って、中継地点で足し合わせる」という等式こそが、本記事の主役であるチャップマン・コルモゴロフ方程式です。
この方程式は、単なる計算テクニックではありません。ここから「$n$ ステップの推移確率は推移行列の $n$ 乗で求まる」という離散時間マルコフ連鎖の基本公式が出てきますし、連続時間に一般化すると「推移確率行列は生成行列の行列指数 $\exp(Qt)$ で書ける」という美しい結論に到達します。応用先も広く、たとえば
- 待ち行列理論・信頼性工学: 窓口の混雑状態や機器の故障・復旧状態が時間とともにどう遷移するかを、生成行列 $Q$ と $\exp(Qt)$ で予測できます。
- 通信路・音声認識・バイオインフォマティクス: 隠れマルコフモデル (HMM) の状態遷移や、DNA 配列の進化モデル (連続時間マルコフ) の計算は、まさにチャップマン・コルモゴロフ方程式の上に成り立っています。
本記事では、この方程式を全確率の法則とマルコフ性から一切ごまかさずに導出し、離散時間の推移行列べき乗、連続時間の前進・後退方程式、そして解 $\exp(Qt)$ までを一本の道筋でつなぎます。
上の図が本記事全体の見取り図です。始点 $i$ から終点 $j$ へ至る道は、途中の時刻 $m$ で必ずどこか一つの中継状態 $k$ を通過します。左半分の「$i$ から $k$ へ $m$ ステップ」と右半分の「$k$ から $j$ へ $n$ ステップ」の確率を掛けて、すべての中継 $k$ について足し合わせる——この足し算こそがチャップマン・コルモゴロフ方程式です。
本記事の内容
- チャップマン・コルモゴロフ方程式の直感とマルコフ性からの導出
- 離散時間での $\bm{P}^{(m+n)} = \bm{P}^{(m)}\bm{P}^{(n)}$ と $\bm{P}^{(n)} = \bm{P}^n$ の証明
- 連続時間での生成行列 $\bm{Q}$、前進方程式・後退方程式、解 $\bm{P}(t)=\exp(\bm{Q}t)$ の導出
- Python による推移行列のべき乗と
scipy.linalg.expmを使った半群性 $\bm{P}(s+t)=\bm{P}(s)\bm{P}(t)$ の数値確認
前提知識
この記事を読む前に、以下の記事を読んでおくと理解が深まります。
- マルコフ連鎖とは？ — 状態・推移確率・推移行列の基本
- 連続時間マルコフ連鎖 — 状態が連続時間で遷移するモデル
- 定常分布 — 長時間後に落ち着く分布
チャップマン・コルモゴロフ方程式とは
まず言葉のイメージから入りましょう。地図上の A 地点から C 地点まで行く経路を数えたいとします。A から C へ直接向かうルートを一つずつ数えるのは大変ですが、「必ず途中で B ライン（中継都市の集まり）を通る」とわかっていれば、話は単純になります。「A からある中継都市 $k$ まで行き、その $k$ から C まで行く」という組み合わせを、すべての中継都市 $k$ について足し合わせればよいのです。
確率の世界でもまったく同じことが起きます。マルコフ連鎖では、時刻 $0$ に状態 $i$ にいたコマが時刻 $m+n$ に状態 $j$ にいる確率を考えます。この道のりは、時刻 $m$ に「どこか」の状態 $k$ を必ず通過します。したがって「$m$ ステップで $i$ から $k$ へ、続く $n$ ステップで $k$ から $j$ へ」という確率を、すべての中継状態 $k$ について足し上げれば、$m+n$ ステップの推移確率が得られます。これがチャップマン・コルモゴロフ方程式です。
記号で書くと、$n$ ステップ推移確率を
$$ p_{ij}^{(n)} = P(X_{m+n} = j \mid X_m = i) $$
と定義したとき、チャップマン・コルモゴロフ方程式は
$$ \begin{equation} p_{ij}^{(m+n)} = \sum_{k} p_{ik}^{(m)}\, p_{kj}^{(n)} \end{equation} $$
という形になります。右辺は「中継状態 $k$ を経由する確率の和」で、左辺は「トータルの推移確率」です。この式が成り立つ理由は、後で全確率の法則とマルコフ性から丁寧に示します。
大切なのは、この足し算がちょうど行列の積の形をしていることです。$(i,k)$ 成分と $(k,j)$ 成分をかけて $k$ で和をとる操作は、行列の掛け算そのものです。したがって推移確率をならべた行列 $\bm{P}^{(n)}$ を使えば、方程式は
$$ \bm{P}^{(m+n)} = \bm{P}^{(m)}\,\bm{P}^{(n)} $$
とすっきり書けます。この行列版が、離散時間では「べき乗」、連続時間では「行列指数」へと姿を変えていきます。まずはその土台として、マルコフ性と全確率の法則を確認しましょう。
マルコフ性と全確率の法則
チャップマン・コルモゴロフ方程式を支える二本の柱が、マルコフ性と全確率の法則です。順に見ていきます。
マルコフ性 — 未来は現在だけで決まる
マルコフ連鎖の心臓部は「未来は、いまの状態だけで決まり、そこに至るまでの過去の経路には依らない」という性質です。すごろくで言えば、「いまマス 5 にいる」という情報さえあれば、そこまでどんな出目でたどり着いたかは、次にどこへ進むかの確率にまったく影響しない、ということです。数式では、
$$ \begin{equation} 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) \end{equation} $$
と書きます。条件のところに過去の履歴 $X_{n-1},\dots,X_0$ をいくら並べても、右辺のように「直前の状態 $X_n = i$ だけ」に潰れてよい、というのがマルコフ性です。
図の破線で囲った「いまの状態 $X_n$」だけが次の状態 $X_{n+1}$ を決め、その手前に伸びる過去の経路（点線の矢印）はいっさい効きません。この「過去を忘れてよい」性質が、次に導出でみるように、推移確率をきれいな積の和にまとめる鍵になります。
さらに本記事では、推移確率が時刻に依らない時間的に一様（時間斉次）なマルコフ連鎖を扱います。つまり
$$ P(X_{n+1} = j \mid X_n = i) = p_{ij} $$
が $n$ に依らず一定であるとします。この一段の推移確率 $p_{ij}$ をならべた行列を推移行列 $\bm{P} = (p_{ij})$ と呼びます。各行は「状態 $i$ からどこかへ必ず遷移する」ので、$\sum_j p_{ij} = 1$（行和が 1）を満たします。このような行列を確率行列（stochastic matrix）といいます。
全確率の法則 — 中継地点で場合分けする
もう一本の柱が全確率の法則です。ある事象 $B$ が起こる確率を、途中で起こりうる互いに排反な事象 $\{A_k\}$（全部合わせると全体になる）で場合分けして、
$$ \begin{equation} P(B) = \sum_k P(B \mid A_k)\, P(A_k) \end{equation} $$
と足し上げてよい、という定理です。マルコフ連鎖では、「時刻 $m$ にどの状態にいるか」がちょうどこの排反な場合分け $\{A_k\} = \{X_m = k\}$ を与えます。コマは時刻 $m$ に必ずどれか一つの状態にいるので、これらは排反かつ網羅的です。
図のように標本空間 $\Omega$ を隙間なく（網羅的に）かつ重なりなく（排反に）分割する事象 $A_k$ を用意すれば、知りたい確率 $P(B)$ を各ピースごとの寄与 $P(B\mid A_k)P(A_k)$ の和に分解できます。マルコフ連鎖ではこの分割を「時刻 $m$ の状態」でとるのが要点です。
この二本の柱、「未来は現在だけで決まる（マルコフ性）」と「途中の 1 点で場合分けしてよい（全確率の法則）」を組み合わせると、チャップマン・コルモゴロフ方程式が自然に姿を現します。次のセクションで、その導出を一行ずつたどりましょう。
チャップマン・コルモゴロフ方程式の導出
ゴールは、式 (1) の
$$ p_{ij}^{(m+n)} = \sum_{k} p_{ik}^{(m)}\, p_{kj}^{(n)} $$
を、全確率の法則とマルコフ性から導くことです。表記を簡単にするため、時間斉次性を使って基準時刻を $0$ にとり、$p_{ij}^{(n)} = P(X_n = j \mid X_0 = i)$ とします。
まず、$m+n$ ステップの推移確率を、時刻 $m$ の状態で場合分けします。全確率の法則（式 (4)）を条件 $X_0 = i$ のもとで適用すると、時刻 $m$ にどの状態 $k$ を通るかで分解できます。
$$ p_{ij}^{(m+n)} = P(X_{m+n} = j \mid X_0 = i) = \sum_k P(X_{m+n} = j,\, X_m = k \mid X_0 = i) $$
ここで、同時確率 $P(X_{m+n}=j, X_m=k \mid X_0=i)$ を、条件付き確率の定義（乗法定理）で二つに分けます。「$X_m=k$ になる確率」と「その上で $X_{m+n}=j$ になる確率」の積です。
$$ p_{ij}^{(m+n)} = \sum_k P(X_{m+n} = j \mid X_m = k,\, X_0 = i)\, P(X_m = k \mid X_0 = i) $$
ここでマルコフ性を使います。第一因子 $P(X_{m+n}=j \mid X_m=k, X_0=i)$ は、条件に「現在 $X_m=k$」と「過去 $X_0=i$」の両方が入っていますが、マルコフ性（式 (2)）により過去 $X_0=i$ の情報は落とせます。つまり、
$$ P(X_{m+n} = j \mid X_m = k,\, X_0 = i) = P(X_{m+n} = j \mid X_m = k) $$
が成り立ちます。さらに時間斉次性から、これは基準時刻を $0$ にずらして $P(X_n = j \mid X_0 = k) = p_{kj}^{(n)}$ と書けます。一方、第二因子はそのまま $P(X_m = k \mid X_0 = i) = p_{ik}^{(m)}$ です。これらを代入すると、
$$ p_{ij}^{(m+n)} = \sum_k p_{kj}^{(n)}\, p_{ik}^{(m)} = \sum_k p_{ik}^{(m)}\, p_{kj}^{(n)} $$
となり、目標のチャップマン・コルモゴロフ方程式 (1) が導けました。導出で本質的に効いているのは、「過去を忘れてよい」というマルコフ性のおかげで、第一因子が中継状態 $k$ だけの関数 $p_{kj}^{(n)}$ に単純化できたことです。もしマルコフ性がなければ、第一因子は $i$ にも依存してしまい、きれいな積の和にはなりません。
この和は行列積の定義そのものなので、推移確率行列 $\bm{P}^{(n)} = (p_{ij}^{(n)})$ を用いて
$$ \begin{equation} \bm{P}^{(m+n)} = \bm{P}^{(m)}\,\bm{P}^{(n)} \end{equation} $$
と行列形式で書けます。ここまでは離散・連続を問わず（適切に和を積分に置き換えれば）成り立つ一般的な関係です。
図が示すとおり、$\sum_k p_{ik}^{(m)} p_{kj}^{(n)}$ は「左行列 $\bm{P}^{(m)}$ の $i$ 行」と「右行列 $\bm{P}^{(n)}$ の $j$ 列」の内積にほかなりません。中継状態 $k$ での和が、そのまま行列積の定義に一致するのです。だからこそチャップマン・コルモゴロフ方程式は $\bm{P}^{(m+n)} = \bm{P}^{(m)}\bm{P}^{(n)}$ という簡潔な行列式に化けます。次は、この行列版から離散時間の具体的な計算公式を引き出しましょう。
離散時間 — nステップ推移確率は推移行列のべき乗
式 (5) の $\bm{P}^{(m+n)} = \bm{P}^{(m)}\bm{P}^{(n)}$ に、具体的な数を入れてみます。ここから、$n$ ステップ推移確率行列が一段の推移行列 $\bm{P}$ の $n$ 乗になる、という気持ちのよい結論が出ます。
まず、1 ステップの推移確率行列は定義そのもので $\bm{P}^{(1)} = \bm{P}$ です。また、$0$ ステップでは状態は動かないので、$p_{ij}^{(0)} = \delta_{ij}$（同じ状態なら 1、違えば 0）となり、$\bm{P}^{(0)} = \bm{I}$（単位行列）です。
ここで式 (5) で $m=n=1$ とおくと、
$$ \bm{P}^{(2)} = \bm{P}^{(1)}\bm{P}^{(1)} = \bm{P}\bm{P} = \bm{P}^2 $$
が得られます。同様に $m=2, n=1$ とすれば $\bm{P}^{(3)} = \bm{P}^{(2)}\bm{P}^{(1)} = \bm{P}^2 \bm{P} = \bm{P}^3$ です。この操作を繰り返すと、数学的帰納法により
$$ \begin{equation} \bm{P}^{(n)} = \bm{P}^n \end{equation} $$
が任意の $n \ge 0$ で成り立ちます。つまり「$n$ ステップ後にどの状態からどの状態へ移るか」の確率は、一段の推移行列 $\bm{P}$ をただ $n$ 回かけ算するだけで全部わかる、ということです。冒頭のすごろくの例で「5 手後」を計算するのに組み合わせを数え上げる必要はなく、$5 \times 5$ 行列を 5 乗すれば済むわけです。
帰納法の中身を明示しておきましょう。$\bm{P}^{(n)} = \bm{P}^n$ が成り立つと仮定します（帰納法の仮定）。このとき式 (5) で $m=n, n=1$ とおくと、
$$ \bm{P}^{(n+1)} = \bm{P}^{(n)}\bm{P}^{(1)} = \bm{P}^n \cdot \bm{P} = \bm{P}^{n+1} $$
となり、$n+1$ でも成立します。基底 $n=1$ は $\bm{P}^{(1)}=\bm{P}$ で成立しているので、すべての $n$ で成立します。
状態分布の時間発展
推移確率行列がわかると、状態の確率分布の時間発展も一発で書けます。時刻 $0$ の分布を行ベクトル $\bm{\pi}_0 = (\pi_0(i))$（$\pi_0(i) = P(X_0 = i)$）で表すと、時刻 $n$ の分布は
$$ \begin{equation} \bm{\pi}_n = \bm{\pi}_0\, \bm{P}^n \end{equation} $$
で求まります。成分で書けば $\pi_n(j) = \sum_i \pi_0(i)\, p_{ij}^{(n)}$ で、これも「初期状態 $i$ で場合分けして足す」という全確率の法則の形です。ここで $n \to \infty$ の極限で $\bm{\pi}_n$ が落ち着く先が、前提記事で扱った定常分布です。定常分布 $\bm{\pi}$ は $\bm{\pi}\bm{P} = \bm{\pi}$ を満たし、まさにチャップマン・コルモゴロフ方程式が保証する「行列べき乗の極限」として現れます。
離散時間では、時間が $0,1,2,\dots$ と飛び飛びに進むので、べき乗という離散的な操作で十分でした。では、時間が連続的に流れる場合はどうなるでしょうか。「$\bm{P}$ の連続回のべき乗」とは何を意味するのか——この問いが、次の連続時間マルコフ連鎖へと導いてくれます。
連続時間 — 生成行列と行列指数
連続時間マルコフ連鎖では、状態が任意の実数時刻 $t \ge 0$ で遷移します。推移確率も時刻の連続関数になり、
$$ p_{ij}(t) = P(X(s+t) = j \mid X(s) = i) $$
を $(i,j)$ 成分とする推移確率行列 $\bm{P}(t)$ を考えます（時間斉次なので基準時刻 $s$ には依りません）。
半群性（連続版チャップマン・コルモゴロフ方程式）
離散のときと同じ導出をたどると（和はそのまま、時刻を実数にするだけ）、連続時間でもチャップマン・コルモゴロフ方程式が成り立ちます。
$$ \begin{equation} \bm{P}(s+t) = \bm{P}(s)\,\bm{P}(t) \end{equation} $$
これは「時間 $s+t$ ぶんの遷移は、まず $s$ ぶん遷移してから $t$ ぶん遷移するのと同じ」という自然な主張です。$\bm{P}(0) = \bm{I}$ とあわせて、この性質を数学では半群性 (semigroup property) と呼びます。指数関数が $e^{s+t} = e^s e^t$ を満たすのとまったく同じ構造で、この類似が後で $\exp(\bm{Q}t)$ という解に直結します。
生成行列 Q の定義
離散時間の主役が一段の推移行列 $\bm{P}$ だったのに対し、連続時間の主役は生成行列（生成作用素、$Q$ 行列とも）$\bm{Q}$ です。連続時間では「一段」という最小の時間刻みが存在しないので、代わりに「ごく短い時間 $h$ のあいだにどれくらいの割合（レート）で遷移するか」を考えます。
$t=0$ で $\bm{P}(0) = \bm{I}$ ですから、微小時間 $h$ での推移行列 $\bm{P}(h)$ は単位行列からのわずかなずれです。このずれの一次の傾きとして生成行列を定義します。
$$ \begin{equation} \bm{Q} = \left.\frac{d\bm{P}(t)}{dt}\right|_{t=0} = \lim_{h \to 0}\frac{\bm{P}(h) – \bm{I}}{h} \end{equation} $$
言い換えると、微小時間 $h$ に対して
$$ \begin{equation} \bm{P}(h) = \bm{I} + \bm{Q}h + o(h) \end{equation} $$
と一次近似できる、ということです。生成行列の各成分 $q_{ij}$ には明快な意味があります。
- 非対角成分 $q_{ij}\ (i \ne j)$: 状態 $i$ から状態 $j$ への単位時間あたりの遷移レート。式 (11) より $p_{ij}(h) \approx q_{ij}h$ なので、短い時間 $h$ に $i$ から $j$ へ移る確率は $q_{ij}h$ に比例します。当然 $q_{ij} \ge 0$ です。
- 対角成分 $q_{ii}$: 状態 $i$ から出ていく総レートの符号を反転したもの。各行の確率の和は常に 1（$\sum_j p_{ij}(h) = 1$）なので、両辺を $h$ で微分して $h\to 0$ とすると $\sum_j q_{ij} = 0$、すなわち
$$ q_{ii} = -\sum_{j \ne i} q_{ij} $$
となります。したがって生成行列の各行の和は 0 です。これは推移行列の「行和が 1」に対応する、連続時間版の保存則です。
図は後述の天気の例（$q_{01}=0.5$）で、遷移確率 $p_{01}(h)$ の厳密値（曲線）と一次近似 $q_{01}h$（破線）を重ねたものです。原点 $h=0$ での接線の傾きがちょうど生成レート $q_{01}$ になっており、微小時間では両者がぴったり寄り添います。「$\bm{P}(h)=\bm{I}+\bm{Q}h+o(h)$」という展開が、$\bm{Q}$ を単位行列からのずれの傾きとして定義していることが視覚的に確認できます。
前進方程式と後退方程式の導出
さて、半群性（式 (9)）と生成行列（式 (11)）を組み合わせると、$\bm{P}(t)$ が満たす微分方程式が二通り出てきます。これがコルモゴロフの前進方程式と後退方程式です。
まず半群性で $s \to t$、$t \to h$ と読み替えて、$\bm{P}(t+h) = \bm{P}(t)\bm{P}(h)$ とします。ここに $\bm{P}(h) = \bm{I} + \bm{Q}h + o(h)$ を代入すると、
$$ \bm{P}(t+h) = \bm{P}(t)\big(\bm{I} + \bm{Q}h + o(h)\big) = \bm{P}(t) + \bm{P}(t)\bm{Q}\,h + o(h) $$
となります。ここで $\bm{P}(t)$ を左辺から引いて $h$ で割り、$h \to 0$ の極限をとります。左辺は微分の定義そのものなので、
$$ \frac{d\bm{P}(t)}{dt} = \lim_{h\to 0}\frac{\bm{P}(t+h) – \bm{P}(t)}{h} = \bm{P}(t)\bm{Q} $$
が得られます。これが前進方程式 (forward equation) です。
$$ \begin{equation} \frac{d\bm{P}(t)}{dt} = \bm{P}(t)\,\bm{Q} \end{equation} $$
一方、半群性を逆順に $\bm{P}(t+h) = \bm{P}(h)\bm{P}(t)$ と書いて、こんどは左側の $\bm{P}(h)$ に近似 $\bm{I}+\bm{Q}h$ を代入すると、
$$ \bm{P}(t+h) = \big(\bm{I} + \bm{Q}h + o(h)\big)\bm{P}(t) = \bm{P}(t) + \bm{Q}\bm{P}(t)\,h + o(h) $$
となり、同じ極限操作で後退方程式 (backward equation) が得られます。
$$ \begin{equation} \frac{d\bm{P}(t)}{dt} = \bm{Q}\,\bm{P}(t) \end{equation} $$
前進方程式では $\bm{Q}$ が右から、後退方程式では $\bm{Q}$ が左から掛かる点だけが違います。名前の由来は、前進が「終わりの時刻」側で状態を分解する（未来へ進む）視点、後退が「始めの時刻」側で分解する（過去へ戻る）視点に対応することにあります。行列積は一般に非可換なので二つは別の式ですが、以下で見るように解は共通です。
図で並べたとおり、両者の違いは生成行列 $\bm{Q}$ を右から掛けるか左から掛けるかだけです。半群性の $\bm{P}(t+h)=\bm{P}(t)\bm{P}(h)$ と $\bm{P}(t+h)=\bm{P}(h)\bm{P}(t)$ という二通りの分解に対応しており、次に見るように、どちらも同じ行列指数 $\exp(\bm{Q}t)$ に行き着きます。
解は行列指数 exp(Qt)
式 (12)・(13) はどちらも「微分すると自分自身に $\bm{Q}$ が掛かる」という形で、スカラーの微分方程式 $\frac{dp}{dt} = pq$ の行列版です。スカラーなら解は $p(t) = p(0)e^{qt}$ ですから、行列でも指数関数、すなわち行列指数で解けると予想できます。実際、初期条件 $\bm{P}(0) = \bm{I}$ のもとで、
$$ \begin{equation} \bm{P}(t) = \exp(\bm{Q}t) = \sum_{k=0}^{\infty}\frac{(\bm{Q}t)^k}{k!} = \bm{I} + \bm{Q}t + \frac{(\bm{Q}t)^2}{2!} + \frac{(\bm{Q}t)^3}{3!} + \cdots \end{equation} $$
が唯一の解です。これが本記事の到達点です。確認のため、この $\bm{P}(t)$ が前進方程式を満たすことを見ておきます。式 (14) の級数を項ごとに $t$ で微分すると、
$$ \frac{d}{dt}\exp(\bm{Q}t) = \bm{Q} + \bm{Q}^2 t + \frac{\bm{Q}^3 t^2}{2!} + \cdots = \bm{Q}\left(\bm{I} + \bm{Q}t + \frac{(\bm{Q}t)^2}{2!} + \cdots\right) = \bm{Q}\exp(\bm{Q}t) $$
となります。$\bm{Q}$ は $\exp(\bm{Q}t)$ と可換（同じ $\bm{Q}$ の級数なので）ですから、これは $\bm{Q}\exp(\bm{Q}t) = \exp(\bm{Q}t)\bm{Q}$ とも書け、前進・後退の両方を同時に満たします。また $\exp(\bm{Q}\cdot 0) = \bm{I}$ で初期条件も満たします。
さらに、行列指数はそのまま半群性も満たします。$\bm{Q}s$ と $\bm{Q}t$ は可換なので、スカラー指数と同じ計算で
$$ \exp(\bm{Q}s)\exp(\bm{Q}t) = \exp\big(\bm{Q}(s+t)\big) = \bm{P}(s+t) $$
が成り立ち、連続版チャップマン・コルモゴロフ方程式 (9) と完全に整合します。こうして「離散のべき乗 $\bm{P}^n$」と「連続の行列指数 $\exp(\bm{Q}t)$」が、同じチャップマン・コルモゴロフ方程式から生まれた双子の兄弟だとわかります。実際、$\exp(\bm{Q}t)$ を「$\bm{Q}$ が生成する連続なべき乗」と読めば、両者は完全に対応しています。
ここまでで理論の骨格は完成しました。抽象的な式が本当に成り立つのか、次は小さな具体例と Python で手を動かして確かめましょう。
具体例 — 2状態の天気モデル
もっとも簡単な 2 状態の例で、数式の意味を体感します。天気が「晴れ (状態 0)」と「雨 (状態 1)」の 2 状態を、1 日ごとに次のルールで遷移するとします。
- 晴れの翌日は、$0.8$ の確率で晴れ、$0.2$ の確率で雨
- 雨の翌日は、$0.6$ の確率で晴れ、$0.4$ の確率で雨
推移行列は
$$ \bm{P} = \begin{pmatrix} 0.8 & 0.2 \\ 0.6 & 0.4 \end{pmatrix} $$
です。各行の和が 1 になっていることを確認してください。
図は同じ推移行列を状態遷移図で描いたものです。各状態（晴れ・雨）から出る矢印の確率をすべて足すと、自己ループを含めて必ず 1 になります（晴れ: $0.8+0.2$、雨: $0.6+0.4$）。この「出ていく確率の総和が 1」が確率行列の行和 1 の正体です。では「今日晴れているとき、2 日後に雨である確率」を求めましょう。式 (6) より $\bm{P}^{(2)} = \bm{P}^2$ を計算します。$(0,1)$ 成分（晴れ→雨）は、チャップマン・コルモゴロフ方程式で 1 日目の天気 $k$ で場合分けして、
$$ p_{01}^{(2)} = \sum_k p_{0k}\,p_{k1} = p_{00}p_{01} + p_{01}p_{11} = 0.8 \times 0.2 + 0.2 \times 0.4 = 0.16 + 0.08 = 0.24 $$
となります。「今日晴れ → 明日も晴れ → 明後日雨」の経路 ($0.8 \times 0.2$) と、「今日晴れ → 明日雨 → 明後日雨」の経路 ($0.2 \times 0.4$) を足し合わせた形になっているのがわかります。まさに「中継地点（明日の天気）で場合分けして和をとる」というチャップマン・コルモゴロフの精神そのものです。
行列全体を計算すると
$$ \bm{P}^2 = \begin{pmatrix} 0.8 & 0.2 \\ 0.6 & 0.4 \end{pmatrix}\begin{pmatrix} 0.8 & 0.2 \\ 0.6 & 0.4 \end{pmatrix} = \begin{pmatrix} 0.76 & 0.24 \\ 0.72 & 0.28 \end{pmatrix} $$
となり、たしかに $(0,1)$ 成分は $0.24$ です。手計算とチャップマン・コルモゴロフ方程式の答えが一致しました。この調子で $\bm{P}^n$ を計算していくと、$n$ が大きくなるにつれて各行が同じ値に近づき、定常分布に収束していきます。次は、これを Python で一気に確認しましょう。
Pythonでの実装
理論を数値実験で裏づけます。ここで確認したいのは次の 3 点です。(1) 離散時間で $\bm{P}^{(n)} = \bm{P}^n$ が成り立ち、$n$ が大きいと定常分布に収束すること。(2) 連続時間で生成行列 $\bm{Q}$ から scipy.linalg.expm で $\bm{P}(t) = \exp(\bm{Q}t)$ が計算できること。(3) 半群性 $\bm{P}(s+t) = \bm{P}(s)\bm{P}(t)$ が数値的に成り立つこと。
離散時間：推移行列のべき乗
まず 2 状態の天気モデルで、推移行列を繰り返し掛けていったときの分布の推移を追います。
import numpy as np
import matplotlib, matplotlib.pyplot as plt
# 日本語フォント設定
for cand in ["Hiragino Sans", "Yu Gothic", "Noto Sans CJK JP", "IPAexGothic", "Meiryo"]:
if any(cand == f.name for f in matplotlib.font_manager.fontManager.ttflist):
plt.rcParams["font.family"] = cand
break
plt.rcParams["axes.unicode_minus"] = False
# 推移行列（0:晴れ, 1:雨）
P = np.array([[0.8, 0.2],
[0.6, 0.4]])
# n ステップ推移行列 P^n を n=0..12 で計算
ns = range(0, 13)
Pn_list = [np.linalg.matrix_power(P, n) for n in ns]
# 今日晴れ(状態0)スタートの分布の推移 pi_n = pi_0 P^n
pi0 = np.array([1.0, 0.0])
pi_hist = np.array([pi0 @ Pn for Pn in Pn_list])
print("P^2 =\n", np.linalg.matrix_power(P, 2))
print("P^10 =\n", np.linalg.matrix_power(P, 10))
上の出力で P^2 の $(0,1)$ 成分が手計算どおり $0.24$ になっていること、そして P^10 では 2 つの行がほぼ同じ値 $(0.75, 0.25)$ に揃っていることが確認できます。行が同じ値に揃うのは、初期状態を忘れて定常分布へ収束していく様子を表しています。この $(0.75, 0.25)$ は定常分布 $\bm{\pi} = \bm{\pi}\bm{P}$ の解と一致します。
$\bm{P}^n$ を $n=1,2,4,10$ と並べたヒートマップです。$n$ が小さいうちは 2 つの行（晴れ発・雨発）の値が違いますが、べき乗を重ねるにつれて両行が同じ $(0.75, 0.25)$ へ近づいていきます。$\bm{P}^{10}$ ではもう行の区別がつかず、「どの初期状態から始めても同じ定常分布に落ち着く」ことがひと目でわかります。
次に、分布 $\bm{\pi}_n$ の収束の様子をグラフにします。
# 分布の時間発展を可視化
plt.figure(figsize=(8, 5))
plt.plot(list(ns), pi_hist[:, 0], "o-", label="晴れの確率")
plt.plot(list(ns), pi_hist[:, 1], "s-", label="雨の確率")
plt.axhline(0.75, color="gray", ls="--", alpha=0.7, label="定常分布(晴れ=0.75)")
plt.axhline(0.25, color="gray", ls=":", alpha=0.7, label="定常分布(雨=0.25)")
plt.xlabel("ステップ数 n")
plt.ylabel("確率")
plt.title("P^n による分布の時間発展（今日晴れスタート）")
plt.legend()
plt.grid(True, alpha=0.3)
plt.tight_layout()
plt.savefig("ck_discrete.png", dpi=150)
plt.show()
このグラフから、初日は「晴れ 1.0・雨 0.0」だった分布が、ステップを重ねるごとに滑らかに定常分布 $(0.75, 0.25)$ の点線へ収束していく様子が読み取れます。$n=6$ あたりでほぼ収束しており、初期状態の影響が急速に薄れることがわかります。これはチャップマン・コルモゴロフ方程式が保証する「べき乗の極限」の具体的な現れです。
連続時間：生成行列と行列指数
続いて連続時間版です。晴れ→雨のレート $0.5$、雨→晴れのレート $1.0$（単位: 1/日）とする生成行列
$$ \bm{Q} = \begin{pmatrix} -0.5 & 0.5 \\ 1.0 & -1.0 \end{pmatrix} $$
を考えます。各行の和が 0 になっていることに注意してください。scipy.linalg.expm で $\bm{P}(t) = \exp(\bm{Q}t)$ を計算します。
import numpy as np
from scipy.linalg import expm
import matplotlib.pyplot as plt
# 生成行列（各行の和は0）
Q = np.array([[-0.5, 0.5],
[ 1.0, -1.0]])
# P(t) = exp(Q t) を時間ごとに計算
ts = np.linspace(0, 6, 200)
Pt = np.array([expm(Q * t) for t in ts])
# 初期状態=晴れ の分布 pi(t) = pi0 exp(Qt)
pi0 = np.array([1.0, 0.0])
pit = np.array([pi0 @ expm(Q * t) for t in ts])
plt.figure(figsize=(8, 5))
plt.plot(ts, pit[:, 0], label="晴れの確率 P(t)")
plt.plot(ts, pit[:, 1], label="雨の確率 P(t)")
plt.axhline(2/3, color="gray", ls="--", alpha=0.7, label="定常(晴れ=2/3)")
plt.axhline(1/3, color="gray", ls=":", alpha=0.7, label="定常(雨=1/3)")
plt.xlabel("時刻 t [日]")
plt.ylabel("確率")
plt.title("連続時間: P(t)=exp(Qt) による分布の時間発展")
plt.legend()
plt.grid(True, alpha=0.3)
plt.tight_layout()
plt.savefig("ck_continuous.png", dpi=150)
plt.show()
このグラフでは、離散時間の階段状の発展とは違い、確率が時刻 $t$ の連続関数として滑らかに定常分布 $(2/3, 1/3)$ へ近づいていきます。生成行列 $\bm{Q}$ の各行和が 0 だったおかげで、各時刻で確率の和が常に 1 に保たれている点も重要です。離散のべき乗と連続の行列指数が、同じ「初期状態を忘れて定常へ向かう」振る舞いを示すことが見て取れます。
半群性の数値検証
最後に、連続版チャップマン・コルモゴロフ方程式 $\bm{P}(s+t) = \bm{P}(s)\bm{P}(t)$ が数値的に成り立つことを確認します。
import numpy as np
from scipy.linalg import expm
Q = np.array([[-0.5, 0.5],
[ 1.0, -1.0]])
s, t = 1.3, 2.1
lhs = expm(Q * (s + t)) # P(s+t)
rhs = expm(Q * s) @ expm(Q * t) # P(s)P(t)
print("P(s+t) =\n", lhs)
print("P(s)P(t) =\n", rhs)
print("最大絶対誤差:", np.max(np.abs(lhs - rhs)))
左の $\bm{P}(s+t)$ と中央の $\bm{P}(s)\bm{P}(t)$ は、成分がすべて同じ値になっています。右の絶対誤差パネルはどの成分も $10^{-16}$ 台で、色がほぼ真っ暗（=ゼロ）です。抽象的に導いた半群性が、実際の数値計算でも丸め誤差レベルの精度で成立していることが見て取れます。
出力される最大絶対誤差は $10^{-15}$ 程度（浮動小数点の丸め誤差レベル）になり、$\bm{P}(s+t)$ と $\bm{P}(s)\bm{P}(t)$ が機械精度で一致します。これは、$\bm{Q}s$ と $\bm{Q}t$ が可換なため $\exp(\bm{Q}s)\exp(\bm{Q}t) = \exp(\bm{Q}(s+t))$ が厳密に成り立つことの数値的な裏づけです。抽象的に導いた半群性が、実際の計算でもきちんと成立していることが確かめられました。
べき乗と行列指数の対応
離散のべき乗と連続の行列指数の関係も数値で覗いてみましょう。$\bm{P}(1) = \exp(\bm{Q})$ を「1 日ぶんの推移行列」とみなすと、$\bm{P}(n) = \exp(\bm{Q}n) = \big(\exp(\bm{Q})\big)^n$ となり、連続時間を整数時刻で刻めば離散のべき乗に一致します。
import numpy as np
from scipy.linalg import expm
Q = np.array([[-0.5, 0.5],
[ 1.0, -1.0]])
P1 = expm(Q) # 1日ぶんの推移行列
n = 4
lhs = expm(Q * n) # P(4) = exp(4Q)
rhs = np.linalg.matrix_power(P1, n) # (exp Q)^4
print("exp(4Q) =\n", lhs)
print("(expQ)^4 =\n", rhs)
print("最大絶対誤差:", np.max(np.abs(lhs - rhs)))
図は連続時間の雨確率 $[\exp(\bm{Q}t)]$（曲線）と、整数時刻でサンプリングした $[(\exp\bm{Q})^n]$（四角い点）を重ねたものです。四角い点はぴったり連続曲線の上に乗っており、連続時間を 1 日刻みで区切ると推移行列 $\bm{P}(1)=\exp(\bm{Q})$ をもつ離散マルコフ連鎖になることがわかります。離散のべき乗と連続の行列指数がチャップマン・コルモゴロフ方程式を通じて地続きだと、数値の上でも確認できます。
このコードでも両者は機械精度で一致し、$\exp(\bm{Q}n) = \big(\exp(\bm{Q})\big)^n$ が確認できます。連続時間マルコフ連鎖を整数時刻でサンプリングすると、推移行列 $\bm{P}(1) = \exp(\bm{Q})$ をもつ離散時間マルコフ連鎖になる、という事実がこの実験に対応しています。離散と連続がチャップマン・コルモゴロフ方程式を通じて完全に地続きであることが、数値の上でも見えてきます。
まとめ
本記事では、チャップマン・コルモゴロフ方程式を出発点に、離散時間の推移行列べき乗から連続時間の行列指数までを一本の道筋でつなぎました。
- チャップマン・コルモゴロフ方程式: 「途中の 1 点で場合分けして中継状態で和をとる」という等式 $p_{ij}^{(m+n)} = \sum_k p_{ik}^{(m)} p_{kj}^{(n)}$。全確率の法則とマルコフ性から導かれ、行列形式では $\bm{P}^{(m+n)} = \bm{P}^{(m)}\bm{P}^{(n)}$。
- 離散時間: 帰納法により $n$ ステップ推移確率は推移行列のべき乗 $\bm{P}^{(n)} = \bm{P}^n$。分布は $\bm{\pi}_n = \bm{\pi}_0 \bm{P}^n$ で発展し、極限で定常分布へ収束する。
- 連続時間: 生成行列 $\bm{Q}$（非対角=遷移レート、行和 0）を定義し、微小時間展開 $\bm{P}(h) = \bm{I} + \bm{Q}h + o(h)$ と半群性から前進方程式 $\dot{\bm{P}} = \bm{P}\bm{Q}$・後退方程式 $\dot{\bm{P}} = \bm{Q}\bm{P}$ を導出。解は行列指数 $\bm{P}(t) = \exp(\bm{Q}t)$。
- Python:
matrix_powerとscipy.linalg.expmを使い、$\bm{P}^{(n)} = \bm{P}^n$、定常分布への収束、半群性 $\bm{P}(s+t) = \bm{P}(s)\bm{P}(t)$ を機械精度で確認した。
チャップマン・コルモゴロフ方程式は、マルコフ連鎖の理論を貫く背骨のような存在です。ここで得た「べき乗と行列指数」の視点は、定常分布・エルゴード性・待ち行列・隠れマルコフモデルといった発展的な話題すべての土台になります。
次のステップとして、以下の記事も参考にしてください。
- マルコフ連鎖とは？ — 推移行列と状態遷移の基礎を再確認する
- 連続時間マルコフ連鎖 — 生成行列と滞在時間の指数分布を深掘りする
- 定常分布 — $\bm{P}^n$ の極限として現れる分布を詳しく学ぶ