マルコフ連鎖(離散時間)の基礎

マルコフ連鎖は確率過程の中で最も基本的かつ応用範囲の広いモデルです。「現在の状態さえ分かれば、未来の予測に過去の情報は不要」というマルコフ性により、解析が格段に容易になります。

本記事では、離散時間マルコフ連鎖の定義、遷移確率行列、チャップマン-コルモゴロフ方程式、状態の分類を解説し、Pythonでシミュレーションを行います。

本記事の内容

  • マルコフ性の定義
  • 遷移確率行列 $\bm{P}$
  • チャップマン-コルモゴロフ方程式
  • 状態の分類(再帰・一過)
  • Pythonでのシミュレーションと可視化

前提知識

この記事を読む前に、以下の記事を読んでおくと理解が深まります。

マルコフ連鎖の定義

マルコフ性

離散時間の確率過程 $\{X_n : n = 0, 1, 2, \dots\}$ がマルコフ性を持つとは、任意の $n$ と状態 $i_0, i_1, \dots, i_{n+1}$ に対して:

$$ P(X_{n+1} = i_{n+1} \mid X_n = i_n, X_{n-1} = i_{n-1}, \dots, X_0 = i_0) = P(X_{n+1} = i_{n+1} \mid X_n = i_n) $$

が成立することです。

時間斉次マルコフ連鎖

さらに遷移確率が時間 $n$ に依存しないとき、時間斉次(homogeneous)と呼びます。

$$ P(X_{n+1} = j \mid X_n = i) = p_{ij} \quad \text{($n$ に依存しない)} $$

以下、特に断りがない限り時間斉次を仮定します。

遷移確率行列

定義

状態空間 $\mathcal{S} = \{1, 2, \dots, N\}$ のマルコフ連鎖の遷移確率行列は:

$$ \bm{P} = (p_{ij})_{N \times N}, \quad p_{ij} = P(X_{n+1} = j \mid X_n = i) $$

確率行列の性質

$\bm{P}$ は以下の性質を満たします。

  1. 非負性: $p_{ij} \geq 0$ (すべての $i, j$)
  2. 行和が1: $\sum_{j=1}^{N} p_{ij} = 1$ (すべての $i$)

このような行列を確率行列(stochastic matrix)と呼びます。

$n$ ステップ遷移確率

$n$ ステップ後の遷移確率は:

$$ p_{ij}^{(n)} = P(X_{m+n} = j \mid X_m = i) $$

行列で表すと:

$$ \bm{P}^{(n)} = \bm{P}^n $$

つまり、$n$ ステップ遷移確率行列は $\bm{P}$ の $n$ 乗です。

状態ベクトルの時間発展

時刻 $n$ での状態の確率分布を行ベクトル $\bm{\pi}^{(n)} = (\pi_1^{(n)}, \pi_2^{(n)}, \dots, \pi_N^{(n)})$ で表すと:

$$ \bm{\pi}^{(n+1)} = \bm{\pi}^{(n)} \bm{P} $$

$$ \bm{\pi}^{(n)} = \bm{\pi}^{(0)} \bm{P}^n $$

初期分布 $\bm{\pi}^{(0)}$ と遷移確率行列 $\bm{P}$ がマルコフ連鎖を完全に特定します。

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

導出

$n$ ステップ遷移確率を $m$ ステップと $(n-m)$ ステップに分解できます。

$$ \begin{align} p_{ij}^{(n)} &= P(X_n = j \mid X_0 = i) \\ &= \sum_{k \in \mathcal{S}} P(X_n = j \mid X_m = k) \cdot P(X_m = k \mid X_0 = i) \\ &= \sum_{k \in \mathcal{S}} p_{ik}^{(m)} \cdot p_{kj}^{(n-m)} \end{align} $$

2行目では全確率の公式を用い、3行目ではマルコフ性を用いました。

行列表現では:

$$ \boxed{\bm{P}^{(n)} = \bm{P}^{(m)} \bm{P}^{(n-m)}} $$

これは行列の積 $\bm{P}^n = \bm{P}^m \bm{P}^{n-m}$ そのものです。

意味

「$i$ から $n$ ステップで $j$ に到達する確率」は、「途中の任意の時点 $m$ でどの状態 $k$ を経由するか」のすべての可能性を足し合わせたものです。

状態の分類

到達可能性

状態 $j$ が状態 $i$ から到達可能(accessible)であるとは、ある $n \geq 1$ が存在して $p_{ij}^{(n)} > 0$ が成立することです。$i \to j$ と書きます。

連通(communicating)

$i \to j$ かつ $j \to i$ が成立するとき、$i$ と $j$ は連通すると言い、$i \leftrightarrow j$ と書きます。

連通関係は同値関係であり、状態空間を連通クラスに分割します。すべての状態が一つの連通クラスに属するマルコフ連鎖を既約(irreducible)と呼びます。

再帰状態と一過状態

状態 $i$ から出発して、再び $i$ に戻る確率を $f_i$ とします。

$$ f_i = P(\exists n \geq 1 : X_n = i \mid X_0 = i) $$

  • 再帰状態(recurrent): $f_i = 1$(必ず戻ってくる)
  • 一過状態(transient): $f_i < 1$(戻ってこない確率がある)

再帰性の判定

以下の同値条件が知られています。

$$ \text{状態 } i \text{ が再帰} \iff \sum_{n=0}^{\infty} p_{ii}^{(n)} = \infty $$

$$ \text{状態 } i \text{ が一過} \iff \sum_{n=0}^{\infty} p_{ii}^{(n)} < \infty $$

一過状態 $i$ に対しては、$i$ を訪れる期待回数は有限です。

周期

状態 $i$ の周期 $d_i$ は:

$$ d_i = \gcd\{n \geq 1 : p_{ii}^{(n)} > 0\} $$

$d_i = 1$ のとき、状態 $i$ は非周期(aperiodic)と呼ばれます。$d_i \geq 2$ のとき周期的(periodic)です。

既約なマルコフ連鎖では、すべての状態の周期は同じです。

具体例: 天気モデル

3状態(晴れ、曇り、雨)の天気モデルを考えます。

$$ \bm{P} = \begin{pmatrix} 0.7 & 0.2 & 0.1 \\ 0.3 & 0.4 & 0.3 \\ 0.2 & 0.3 & 0.5 \end{pmatrix} $$

すべての要素が正なので、このマルコフ連鎖は既約かつ非周期です。

3ステップ後の遷移確率行列は:

$$ \bm{P}^3 = \bm{P} \cdot \bm{P} \cdot \bm{P} $$

Pythonでの実装

import numpy as np
import matplotlib.pyplot as plt

np.random.seed(42)

# 遷移確率行列(天気モデル)
P = np.array([
    [0.7, 0.2, 0.1],
    [0.3, 0.4, 0.3],
    [0.2, 0.3, 0.5]
])
states = ['晴れ', '曇り', '雨']

fig, axes = plt.subplots(2, 2, figsize=(14, 10))

# --- (1) マルコフ連鎖のシミュレーション ---
n_steps = 100
n_paths = 5

for path in range(n_paths):
    state = 0  # 晴れから開始
    trajectory = [state]
    for _ in range(n_steps - 1):
        state = np.random.choice(3, p=P[state])
        trajectory.append(state)
    axes[0, 0].step(range(n_steps), trajectory, alpha=0.7, linewidth=1.2)

axes[0, 0].set_yticks([0, 1, 2])
axes[0, 0].set_yticklabels(states)
axes[0, 0].set_xlabel('Time step $n$')
axes[0, 0].set_ylabel('State')
axes[0, 0].set_title('Markov Chain Sample Paths')
axes[0, 0].grid(True, alpha=0.3)

# --- (2) P^n の収束(チャップマン-コルモゴロフ) ---
max_n = 30
p_ij_history = np.zeros((max_n, 3, 3))

Pn = np.eye(3)
for n in range(max_n):
    Pn = Pn @ P
    p_ij_history[n] = Pn

# p_{ij}^{(n)} の時間変化を表示(i=0: 晴れから出発)
colors = ['gold', 'gray', 'blue']
for j in range(3):
    axes[0, 1].plot(range(1, max_n + 1), p_ij_history[:, 0, j],
                   color=colors[j], linewidth=2, label=f'$p_{{晴れ,{states[j]}}}^{{(n)}}$')

# 定常分布(固有値分解で計算)
eigenvalues, eigenvectors = np.linalg.eig(P.T)
idx = np.argmin(np.abs(eigenvalues - 1.0))
pi_stationary = np.real(eigenvectors[:, idx])
pi_stationary = pi_stationary / pi_stationary.sum()

for j in range(3):
    axes[0, 1].axhline(y=pi_stationary[j], color=colors[j], linestyle='--', alpha=0.5)

axes[0, 1].set_xlabel('Step $n$')
axes[0, 1].set_ylabel('$p_{ij}^{(n)}$')
axes[0, 1].set_title('Convergence of $n$-step Transition Probabilities')
axes[0, 1].legend()
axes[0, 1].grid(True, alpha=0.3)

# --- (3) 状態の頻度(エルゴード性の確認) ---
n_steps = 10000
state = 0
trajectory = [state]
for _ in range(n_steps - 1):
    state = np.random.choice(3, p=P[state])
    trajectory.append(state)
trajectory = np.array(trajectory)

# 累積頻度
freq = np.zeros((n_steps, 3))
for n in range(n_steps):
    for j in range(3):
        freq[n, j] = np.mean(trajectory[:n+1] == j)

for j in range(3):
    axes[1, 0].plot(range(n_steps), freq[:, j], color=colors[j], linewidth=1.5,
                   label=f'{states[j]}: 実測')
    axes[1, 0].axhline(y=pi_stationary[j], color=colors[j], linestyle='--', alpha=0.7)

axes[1, 0].set_xlabel('Step $n$')
axes[1, 0].set_ylabel('Empirical Frequency')
axes[1, 0].set_title('Ergodicity: Frequency → Stationary Distribution')
axes[1, 0].legend(fontsize=9)
axes[1, 0].grid(True, alpha=0.3)
axes[1, 0].set_xlim(0, n_steps)

# --- (4) 再帰状態と一過状態の例 ---
# 吸収マルコフ連鎖
P_absorb = np.array([
    [0.5, 0.5, 0.0, 0.0],
    [0.0, 0.0, 0.8, 0.2],
    [0.0, 0.0, 1.0, 0.0],  # 吸収状態
    [0.0, 0.0, 0.0, 1.0],  # 吸収状態
])
states_abs = ['1 (一過)', '2 (一過)', '3 (吸収)', '4 (吸収)']

n_steps = 50
n_paths = 20

# 各パスの最終状態をカウント
final_states = []
for path in range(n_paths):
    state = 0  # 状態1から開始
    trajectory = [state]
    for _ in range(n_steps - 1):
        state = np.random.choice(4, p=P_absorb[state])
        trajectory.append(state)
    final_states.append(trajectory[-1])
    axes[1, 1].step(range(n_steps), trajectory, alpha=0.3, linewidth=0.8)

axes[1, 1].set_yticks([0, 1, 2, 3])
axes[1, 1].set_yticklabels(states_abs, fontsize=8)
axes[1, 1].set_xlabel('Time step $n$')
axes[1, 1].set_ylabel('State')
axes[1, 1].set_title('Absorbing Markov Chain (Transient → Recurrent)')
axes[1, 1].axhline(y=2, color='red', linestyle='--', alpha=0.5, label='吸収状態')
axes[1, 1].axhline(y=3, color='red', linestyle='--', alpha=0.5)
axes[1, 1].legend(fontsize=9)
axes[1, 1].grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig('markov_chain_discrete.png', dpi=150, bbox_inches='tight')
plt.show()

# 数値結果
print("=== 遷移確率行列 P ===")
print(P)

print("\n=== P^10 ===")
print(np.linalg.matrix_power(P, 10).round(4))

print(f"\n=== 定常分布 π ===")
for j in range(3):
    print(f"π({states[j]}) = {pi_stationary[j]:.4f}")

print(f"\n=== チャップマン-コルモゴロフの検証 ===")
P3 = np.linalg.matrix_power(P, 3)
P2 = np.linalg.matrix_power(P, 2)
P1 = P
print(f"P^3 = P^2 × P^1:")
print(f"  直接計算: {P3[0]}")
print(f"  分解計算: {(P2 @ P1)[0]}")
print(f"  一致: {np.allclose(P3, P2 @ P1)}")

まとめ

本記事では、離散時間マルコフ連鎖の基礎について解説しました。

  • マルコフ性: $P(X_{n+1}=j \mid X_n=i, \dots) = P(X_{n+1}=j \mid X_n=i)$(過去によらない)
  • 遷移確率行列 $\bm{P}$: $p_{ij} \geq 0$, $\sum_j p_{ij} = 1$。$n$ ステップ遷移は $\bm{P}^n$
  • チャップマン-コルモゴロフ方程式: $\bm{P}^{(n)} = \bm{P}^{(m)}\bm{P}^{(n-m)}$
  • 状態の分類: 再帰状態($f_i = 1$)と一過状態($f_i < 1$)
  • 既約・非周期のマルコフ連鎖は定常分布に収束する

次のステップとして、以下の記事も参考にしてください。