マルコフ連鎖 - ヘロログ

明日の天気は「今日が晴れかどうか」に大きく依存するが、「3日前に雨だったかどうか」はあまり関係ない。株価の明日の動きは今日の価格に依存するが、1年前の価格の記憶は薄い。このように「未来が現在の状態だけに依存し、過去の経緯には依存しない」という性質をもつ確率的な過程は、自然界や社会のいたるところに現れる。

この性質を数学的に定式化したのがマルコフ連鎖（Markov chain）である。ロシアの数学者アンドレイ・マルコフ（Andrey A. Markov, 1856–1922）は、1906年に従属する確率変数の列に関する最初の論文を発表した。さらに1913年、プーシキンの韻文小説『エヴゲーニイ・オネーギン』の冒頭2万文字を母音と子音に分類し、文字の出現が前の文字に依存する確率構造を分析した。これがマルコフ連鎖の最初の実証的応用である。

現代ではマルコフ連鎖はきわめて幅広い分野で活用されている。Googleの初期の検索アルゴリズムであるPageRankは、ウェブページ間のリンク構造を巨大なマルコフ連鎖として定式化したものである。ベイズ統計学ではマルコフ連鎖モンテカルロ法（MCMC）が複雑な事後分布からのサンプリングに不可欠な手法となっている。そのほか、遺伝子配列の解析、音声認識、金融工学における株価モデリングなど、応用は多岐にわたる。

マルコフ連鎖の定義

マルコフ性

実数値確率変数の列 $X = (X_n)_{n=0,1,2,\dots}$ と集合 $B \subset \mathbb{R}$ に対して、

P(X_{n+1} \in B \mid X_n, X_{n-1}, \dots, X_1, X_0) = P(X_{n+1} \in B \mid X_n)

が成り立つとき、これを確率変数列 $X$ のマルコフ性（Markov property）という。マルコフ性をもつ $X$ をマルコフ連鎖という。

直感的に述べると、現在の状態 $X_n$ が与えられれば、それ以外の過去の情報は未来 $X_{n+1}$ を予測するのに不要であるということである。たとえば天気の例では「明日の天気を予測するには今日の天気だけがわかれば十分で、昨日以前の天気の履歴は追加の情報を与えない」ことに対応する。この性質は「無記憶性」とも呼ばれる。

左辺では $X_0, X_1, \dots, X_n$ というすべての過去の情報を条件にしているのに対し、右辺では現在の状態 $X_n$ のみを条件にしている。マルコフ性が成り立てば、未来の予測においてこの2つの条件付き確率が等しいことを意味する。

状態空間と齊時性

$n$ を時間の変数とみなすとき、 $X_n$ は時点 $n$ における「状態」を表す。 $X_n$ がとりうる値の集合を状態空間（state space） $S$ という。応用上は $S = \{1, 2, \dots, N\}$ のように離散的で有限な状態空間をもつ場合が重要であり、このとき $X$ を有限マルコフ連鎖（finite Markov chain）という。

マルコフ連鎖 $X$ の $m$ ステップ推移確率を

p_m^{(n)}(x, B) := P(X_{n+m} \in B \mid X_n = x)

と定義する。これは「時点 $n$ に状態 $x$ にいたとき、 $m$ ステップ後に集合 $B$ に入る確率」を表す。この推移確率が時点 $n$ によらないとき、すなわち推移の規則が時間とともに変わらないとき、マルコフ連鎖 $X$ は齊時的（time homogeneous）であるという。以下では齊時的有限マルコフ連鎖に限定して議論する。

推移確率行列

状態空間 $S = \{1, 2, \dots, N\}$ 上の齊時的有限マルコフ連鎖 $X = (X_n)_{n=0,1,2,\dots}$ に対して、状態 $i$ から状態 $j$ への1ステップでの推移確率を

p(i, j) = P(X_{n+1} = j \mid X_n = i), \quad i, j \in S

とおく。これを $(i, j)$ 成分にもつ $N \times N$ 行列

Q = \begin{pmatrix} p(1,1) & p(1,2) & \cdots & p(1,N) \\ p(2,1) & p(2,2) & \cdots & p(2,N) \\ \vdots & \vdots & \ddots & \vdots \\ p(N,1) & p(N,2) & \cdots & p(N,N) \end{pmatrix}

を推移確率行列（transition matrix）という。行列の第 $i$ 行が「状態 $i$ から出発したときの行き先の確率分布」を表している。たとえば第1行 $(p(1,1), p(1,2), \dots, p(1,N))$ は「状態1にいるとき、次のステップで各状態に移る確率」を並べたものである。

推移確率行列は次の2つの性質を満たす。

確率行列の性質

① 各成分は非負： $p(i, j) \geq 0$ （すべての $i, j$ に対して）
② 各行の成分の和が1： $\displaystyle\sum_{j \in S} p(i, j) = 1$ （すべての $i$ に対して）

この2つの性質を満たす行列を確率行列（probability matrix）という。①はどの推移確率も0以上であることを、②はある状態から出発して次のステップで必ずどこかの状態に移ることを意味する。

状態確率ベクトルと推移の性質

時点 $n$ における各状態の確率を $p_n(k) := P(X_n = k)$ （ $k \in S$ ）とおくとき、

\boldsymbol{\pi}_n = (p_n(1), \; p_n(2), \; \dots, \; p_n(N))

を時点 $n$ における $X$ の状態確率ベクトル（state probability vector）という。特に $\boldsymbol{\pi}_0$ を初期分布（initial distribution）という。成分の総和が1となるベクトルを確率ベクトルという。

漸化式の導出

時点 $n$ で状態 $k$ にいる確率は、「時点 $n-1$ でどの状態にいたか」を全通り考えて足し合わせることで得られる（全確率の法則）。すなわち

p_n(k) = \sum_{i=1}^{N} p_{n-1}(i) \cdot p(i, k)

これは「前の時点で状態 $i$ にいた確率 $p_{n-1}(i)$ 」に「状態 $i$ から状態 $k$ へ移る確率 $p(i,k)$ 」をかけて、すべての出発状態 $i$ について合計したものである。これを行列の積として書くと、次のコンパクトな漸化式になる。

状態確率の漸化式

\boldsymbol{\pi}_n = \boldsymbol{\pi}_{n-1} Q, \quad n = 1, 2, \dots

この漸化式を繰り返し適用すると、

\boldsymbol{\pi}_n = \boldsymbol{\pi}_{n-1} Q = (\boldsymbol{\pi}_{n-2} Q) Q = \boldsymbol{\pi}_{n-2} Q^2 = \cdots = \boldsymbol{\pi}_0 Q^n

が得られる。すなわち、初期分布 $\boldsymbol{\pi}_0$ と推移確率行列の $n$ 乗 $Q^n$ だけから、任意の時点の状態確率が求められる。

チャップマン–コルモゴロフ方程式

状態 $i$ から $m$ ステップで状態 $j$ に到達する確率 $p_m(i,j)$ を $(i,j)$ 成分にもつ行列を $Q(m)$ と書く。このとき、 $m+l$ ステップの推移確率を、途中のどこかの状態 $k$ を経由して分解できる。具体的には、状態 $i$ から出発して $m$ ステップでまず状態 $k$ に到達し、そこからさらに $l$ ステップで状態 $j$ に到達する経路を、中間状態 $k$ のすべての可能性について足し合わせると

p_{m+l}(i, j) = \sum_{k=1}^{N} p_m(i, k) \, p_l(k, j)

が成り立つ。右辺は行列の積 $Q(m) Q(l)$ の $(i,j)$ 成分そのものであるから、行列の形で書くと次のチャップマン–コルモゴロフ方程式が得られる。

チャップマン–コルモゴロフ方程式

Q(m + l) = Q(m) \, Q(l)

齊時的マルコフ連鎖では $Q(m) = Q^m$ であるから、これは $Q^{m+l} = Q^m Q^l$ という行列のべき乗の性質に他ならない。たとえば、3ステップの推移確率行列は $Q(3) = Q(2) \cdot Q(1) = Q^2 Q = Q^3$ として、1ステップの推移確率行列 $Q$ を3回かけることで得られる。

計算例：2状態マルコフ連鎖

天気の推移モデル

ある地域の天気を「晴れ」（状態1）と「雨」（状態2）の2状態で表す。天気の推移が次の確率に従うとする。

・晴れの日の翌日：晴れ 0.7、雨 0.3
・雨の日の翌日：晴れ 0.4、雨 0.6

この推移確率行列は

Q = \begin{pmatrix} 0.7 & 0.3 \\ 0.4 & 0.6 \end{pmatrix}

である。各行の和が $0.7 + 0.3 = 1$ 、 $0.4 + 0.6 = 1$ となり、確率行列の条件を満たしている。

状態確率ベクトルの計算

今日が晴れだとわかっている場合、初期分布は $\boldsymbol{\pi}_0 = (1, \; 0)$ である。1日後の状態確率ベクトルは

\boldsymbol{\pi}_1 = \boldsymbol{\pi}_0 Q = (1, \; 0) \begin{pmatrix} 0.7 & 0.3 \\ 0.4 & 0.6 \end{pmatrix} = (0.7, \; 0.3)

すなわち、明日は70%の確率で晴れ、30%の確率で雨である。さらに2日後は

\begin{aligned} \boldsymbol{\pi}_2 &= \boldsymbol{\pi}_1 Q = (0.7, \; 0.3) \begin{pmatrix} 0.7 & 0.3 \\ 0.4 & 0.6 \end{pmatrix} \\ &= (0.7 \times 0.7 + 0.3 \times 0.4, \;\; 0.7 \times 0.3 + 0.3 \times 0.6) \\ &= (0.61, \; 0.39) \end{aligned}

3日後は $\boldsymbol{\pi}_3 = (0.583, \; 0.417)$ 、4日後は $\boldsymbol{\pi}_4 = (0.575, \; 0.425)$ と推移し、次第にある一定の値に収束していく様子が見てとれる。この収束先が次に述べる定常分布である。

定常分布

定義

状態確率ベクトル $\boldsymbol{\pi}_n$ が $n \to \infty$ で極限に収束する場合を考えよう。

\boldsymbol{\pi} := \lim_{n \to \infty} \boldsymbol{\pi}_n

漸化式 $\boldsymbol{\pi}_n = \boldsymbol{\pi}_{n-1} Q$ の両辺で $n \to \infty$ とすると、次の関係が得られる。

定常分布の条件

\boldsymbol{\pi} = \boldsymbol{\pi} Q

すなわち、 $\boldsymbol{\pi}$ は行列 $Q$ の固有値1に対する左固有ベクトルである。この条件の意味を考えてみよう。 $\boldsymbol{\pi} = \boldsymbol{\pi} Q$ とは「 $\boldsymbol{\pi}$ に推移確率行列 $Q$ をかけても自分自身と同じになる」ということである。したがって、もしマルコフ連鎖の初期分布がちょうど $\boldsymbol{\pi}$ であれば、1回の推移によっても状態確率は変化しない。さらに長時間経過後に状態確率が $\boldsymbol{\pi}$ に収束してしまうと、その後は何回推移しても状態確率は $\boldsymbol{\pi}$ のまま変わらない。この $\boldsymbol{\pi}$ を $X$ の定常分布（stationary distribution）という。

定常分布の求め方

定常分布 $\boldsymbol{\pi}$ を求めるには、 $\boldsymbol{\pi} = \boldsymbol{\pi} Q$ を連立方程式として解けばよい。手順は次の通りである。

① $\boldsymbol{\pi} = \boldsymbol{\pi} Q$ を変形して $\boldsymbol{\pi}(Q - I) = \mathbf{0}$ とする（ $I$ は単位行列）
② これは $N$ 個の未知数に対して $N$ 本の方程式だが、1本は他の方程式から導ける（線形従属）
③ そこで1本を確率ベクトルの条件 $\pi_1 + \pi_2 + \cdots + \pi_N = 1$ に置き換えて解く

2状態の天気モデルの定常分布

先ほどの天気モデルに戻ろう。定常分布を $\boldsymbol{\pi} = (a, \; b)$ とおくと、 $\boldsymbol{\pi} = \boldsymbol{\pi} Q$ より

\begin{cases} a = 0.7a + 0.4b \\ b = 0.3a + 0.6b \end{cases}

第1式を変形すると $0.3a = 0.4b$ 、すなわち $a = \dfrac{4b}{3}$ を得る。これと $a + b = 1$ をあわせて

\dfrac{4b}{3} + b = 1 \quad \Rightarrow \quad \dfrac{7b}{3} = 1 \quad \Rightarrow \quad b = \dfrac{3}{7}

よって $a = \dfrac{4}{7}$ であり、定常分布は

\boldsymbol{\pi} = \left(\dfrac{4}{7}, \; \dfrac{3}{7}\right) \approx (0.571, \; 0.429)

となる。この地域では長期的にみると約57%の日が晴れ、約43%の日が雨であると解釈できる。先ほどの計算で確認したように、初期分布 $\boldsymbol{\pi}_0 = (1, 0)$ から出発した状態確率ベクトルが $(0.575, 0.425)$ と、4日目にはほぼこの定常分布に近づいていた。

吸収状態に注意

推移確率行列に $p(i, i) = 1$ となる行がある場合、状態 $i$ は吸収状態（absorbing state）と呼ばれ、一度その状態に入ると二度と離れることができない。吸収状態をもつマルコフ連鎖では、定常分布が初期分布に依存する場合がある。

計算例：3状態マルコフ連鎖

製品品質の推移モデル

ある工場の製品品質を「良品」（状態1）、「標準」（状態2）、「不良」（状態3）の3段階で評価する。月ごとの品質の推移が次の推移確率行列に従うとする。

Q = \begin{pmatrix} 0.5 & 0.4 & 0.1 \\ 0.2 & 0.6 & 0.2 \\ 0.1 & 0.4 & 0.5 \end{pmatrix}

各行の和は $0.5+0.4+0.1=1$ 、 $0.2+0.6+0.2=1$ 、 $0.1+0.4+0.5=1$ で確率行列の条件を満たしている。

定常分布の計算

定常分布を $\boldsymbol{\pi} = (a, \; b, \; c)$ とおき、 $\boldsymbol{\pi} = \boldsymbol{\pi} Q$ と $a + b + c = 1$ から連立方程式を立てる。

\begin{cases} a = 0.5a + 0.2b + 0.1c \\ b = 0.4a + 0.6b + 0.4c \\ c = 0.1a + 0.2b + 0.5c \\ a + b + c = 1 \end{cases}

第1式を変形すると $0.5a = 0.2b + 0.1c$ 、すなわち $5a = 2b + c$ …(I) を得る。第3式を変形すると $0.5c = 0.1a + 0.2b$ 、すなわち $5c = a + 2b$ …(II) を得る。

(I) と (II) を辺々引くと $5a - 5c = c - a$ 、すなわち $6a = 6c$ より $a = c$ となる。(I) に代入して $5a = 2b + a$ 、すなわち $b = 2a$ を得る。 $a + b + c = 1$ に代入すると

a + 2a + a = 1 \quad \Rightarrow \quad 4a = 1 \quad \Rightarrow \quad a = \dfrac{1}{4}

よって定常分布は

\boldsymbol{\pi} = \left(\dfrac{1}{4}, \; \dfrac{1}{2}, \; \dfrac{1}{4}\right)

となる。長期的にみると、品質が「良品」の確率は25%、「標準」の確率は50%、「不良」の確率は25%に落ち着く。

パラメータ推定

推移確率行列 $Q$ が未知パラメータ $\theta$ に依存しているとする。マルコフ連鎖の観測データ $x_0, x_1, \dots, x_n$ が与えられたとき、そのデータの実現確率は

P(X_0 = x_0, X_1 = x_1, \dots, X_n = x_n) = p_0(x_0) \prod_{j=1}^{n} p_\theta(x_{j-1}, x_j)

と書ける。ただし、初期分布 $p_0$ はパラメータに依存しないとする。この式の意味は明快で、マルコフ性により「 $x_0$ から出発して $x_1$ へ、 $x_1$ から $x_2$ へ、…」という各ステップの推移確率を順にかけ合わせたものがデータ全体の同時確率になる。対数をとってパラメータに依存しない項を無視すれば、対数尤度関数は

マルコフ連鎖の対数尤度関数

\ell_n(\theta) = \sum_{j=1}^{n} \log p_\theta(x_{j-1}, x_j)

$\theta$ の最尤推定値 $\hat{\theta}$ は、尤度方程式

\dfrac{\partial}{\partial \theta} \ell_n(\hat{\theta}) = 0

の解として求められる。

推移回数による推定

推移確率行列 $Q$ のパラメータが具体的な値（ $\theta$ に依存しない）の場合、観測データから直接推移確率を推定できる。状態 $i$ から状態 $j$ への推移が観測された回数を $n_{ij}$ 、状態 $i$ からの推移の総回数を $n_i = \sum_j n_{ij}$ とすると、推移確率の最尤推定量は

\hat{p}(i, j) = \dfrac{n_{ij}}{n_i}

である。これは直感的にも自然な推定量であり、「状態 $i$ からの推移のうち、状態 $j$ に移った割合」を推移確率の推定値とするものである。

推移確率の推定

2状態のマルコフ連鎖の観測データ列 $(x_0, x_1, \dots, x_8) = (1, 1, 2, 2, 1, 2, 1, 1, 2)$ から推移確率行列を推定する。

データから推移回数を数える。

推移	回数
$1 \to 1$	2
$1 \to 2$	3
$2 \to 1$	2
$2 \to 2$	1

状態1からの推移は合計5回、状態2からの推移は合計3回であるから、推移確率の推定値は

\hat{Q} = \begin{pmatrix} \dfrac{2}{5} & \dfrac{3}{5} \\[6pt] \dfrac{2}{3} & \dfrac{1}{3} \end{pmatrix}

となる。

練習問題

問1. 状態空間

S = \{1, 2\}

のマルコフ連鎖が次の推移確率行列に従うとする。

Q = \begin{pmatrix} 0.8 & 0.2 \\ 0.5 & 0.5 \end{pmatrix}

初期分布を

\boldsymbol{\pi}_0 = (0, \; 1)

とするとき、以下を求めよ。
[1] 2ステップ後の状態確率ベクトル

\boldsymbol{\pi}_2

[2] 定常分布

\boldsymbol{\pi}

[1]

\begin{aligned} \boldsymbol{\pi}_1 &= (0, \; 1) \begin{pmatrix} 0.8 & 0.2 \\ 0.5 & 0.5 \end{pmatrix} = (0.5, \; 0.5) \\[6pt] \boldsymbol{\pi}_2 &= (0.5, \; 0.5) \begin{pmatrix} 0.8 & 0.2 \\ 0.5 & 0.5 \end{pmatrix} = (0.65, \; 0.35) \end{aligned}

[2] $\boldsymbol{\pi} = (a, b)$ とおいて $\boldsymbol{\pi} = \boldsymbol{\pi} Q$ より $a = 0.8a + 0.5b$ 。変形して $0.2a = 0.5b$ 、すなわち $a = \dfrac{5b}{2}$ 。 $a + b = 1$ に代入して $\dfrac{7b}{2} = 1$ より $b = \dfrac{2}{7}$ 、 $a = \dfrac{5}{7}$ 。

\boldsymbol{\pi} = \left(\dfrac{5}{7}, \; \dfrac{2}{7}\right) \approx (0.714, \; 0.286)

問2. 2状態マルコフ連鎖の観測データが

(x_0, x_1, \dots, x_9) = (2, 1, 1, 2, 1, 2, 2, 1, 2, 2)

であるとき、以下を求めよ。
[1] 推移確率行列の最尤推定量

\hat{Q}

[2] 推定した推移確率行列にもとづく定常分布

[1] 推移回数を数える。

推移	回数
$1 \to 1$	1
$1 \to 2$	3
$2 \to 1$	3
$2 \to 2$	2

状態1からの推移は計4回、状態2からの推移は計5回であるから

\hat{Q} = \begin{pmatrix} \dfrac{1}{4} & \dfrac{3}{4} \\[6pt] \dfrac{3}{5} & \dfrac{2}{5} \end{pmatrix}

[2] 定常分布を $(a, b)$ とおくと $\dfrac{3}{4}a = \dfrac{3}{5}b$ より $a = \dfrac{4b}{5}$ 。 $a + b = 1$ に代入して $\dfrac{9b}{5} = 1$ より $b = \dfrac{5}{9}$ 、 $a = \dfrac{4}{9}$ 。

\boldsymbol{\pi} = \left(\dfrac{4}{9}, \; \dfrac{5}{9}\right) \approx (0.444, \; 0.556)

まとめ

項目	内容
マルコフ性	$P(X_{n+1} \in B \mid X_n, \dots, X_0) = P(X_{n+1} \in B \mid X_n)$
推移確率行列 $Q$	$p(i,j) \geq 0$ 、各行の和が1の確率行列
状態確率ベクトル	$\boldsymbol{\pi}_n = \boldsymbol{\pi}_0 Q^n$
定常分布	$\boldsymbol{\pi} = \boldsymbol{\pi} Q$ を満たす確率ベクトル
定常分布の求め方	$\boldsymbol{\pi}(Q - I) = \mathbf{0}$ と $\sum \pi_i = 1$ を連立して解く
推移確率の推定	$\hat{p}(i,j) = n_{ij} / n_i$ （推移回数の比）
対数尤度関数	$\ell_n(\theta) = \sum_{j=1}^{n} \log p_\theta(x_{j-1}, x_j)$

Python実装

マルコフ連鎖の状態確率ベクトルの計算、定常分布の導出、観測データからの推移確率推定をPythonで実装する。

markov_chain.py

import numpy as np

def stationary_distribution(Q):
    """推移確率行列Qの定常分布を求める"""
    N = Q.shape[0]
    # pi(Q - I) = 0 と sum(pi) = 1 を連立して解く
    A = (Q.T - np.eye(N))
    A[-1, :] = 1  # 最後の方程式を sum(pi) = 1 に置換
    b = np.zeros(N)
    b[-1] = 1
    pi = np.linalg.solve(A, b)
    return pi

def estimate_transition_matrix(data, num_states):
    """観測データから推移確率行列を推定する"""
    n_ij = np.zeros((num_states, num_states))
    for t in range(len(data) - 1):
        i, j = data[t] - 1, data[t+1] - 1
        n_ij[i, j] += 1
    row_sums = n_ij.sum(axis=1, keepdims=True)
    Q_hat = n_ij / row_sums
    return Q_hat, n_ij

markov_chain_examples.py

# ========== 2状態の天気モデル ==========
print("【2状態マルコフ連鎖：天気モデル】")
Q = np.array([[0.7, 0.3],
              [0.4, 0.6]])

pi_0 = np.array([1.0, 0.0])
print(f"初期分布: {pi_0}")
for n in range(1, 6):
    pi_n = pi_0 @ np.linalg.matrix_power(Q, n)
    print(f"  π_{n} = {np.round(pi_n, 4)}")

pi = stationary_distribution(Q)
print(f"定常分布: {np.round(pi, 4)}")

# ========== 3状態の品質モデル ==========
print("\n【3状態マルコフ連鎖：品質モデル】")
Q3 = np.array([[0.5, 0.4, 0.1],
               [0.2, 0.6, 0.2],
               [0.1, 0.4, 0.5]])

pi3 = stationary_distribution(Q3)
print(f"定常分布: {pi3}")

# Q^n の収束を確認
for n in [1, 5, 10, 30]:
    Q_n = np.linalg.matrix_power(Q3, n)
    print(f"  Q^{n:2d} 第1行 = {np.round(Q_n[0], 4)}")

# ========== パラメータ推定 ==========
print("\n【パラメータ推定】")
data = [1, 1, 2, 2, 1, 2, 1, 1, 2]
Q_hat, n_ij = estimate_transition_matrix(data, num_states=2)
print(f"推移回数:\n{n_ij.astype(int)}")
print(f"推移確率行列の推定値:\n{np.round(Q_hat, 4)}")

【2状態マルコフ連鎖：天気モデル】初期分布: [1. 0.] π_1 = [0.7 0.3 ] π_2 = [0.61 0.39] π_3 = [0.583 0.417 ] π_4 = [0.5749 0.4251] π_5 = [0.5725 0.4275] 定常分布: [0.5714 0.4286] 【3状態マルコフ連鎖：品質モデル】定常分布: [0.25 0.5 0.25] Q^ 1 第1行 = [0.5 0.4 0.1] Q^ 5 第1行 = [0.2552 0.4998 0.245 ] Q^10 第1行 = [0.2501 0.5 0.2499] Q^30 第1行 = [0.25 0.5 0.25] 【パラメータ推定】推移回数: [[2 3] [2 1]] 推移確率行列の推定値: [[0.4 0.6 ] [0.6667 0.3333]]