マルコフ連鎖：未来は過去を忘れる

昨日だけを覚える習慣

第8章では、すでに知っていることをグラフとして描く方法を学びました。変数をノード、「依存する」という関係を矢印で表し、マルコフ因数分解によって結合分布を図から直接読み取ることができました。そのマルコフという言葉が、今度は新しい形で登場します。これまでのグラフはスナップショットでした。固定された変数の集合を一度に扱うものでした。この章では、あらゆる物語がいつか必要とする一つの要素を加えます。時間です。

Chibanyのランチに新しいひねりを加えましょう。最初の章では、学生がChibanyに弁当を持参し、中身を推測しようとしました。今回は学生が毎日二つの弁当を持ってきます。一つはとんかつ、もう一つはハンバーガーで、Chibanyはどちらを食べるかを選びます。隠れた中身を推測するのではなく、実際の日々の選択の連続を観察できるようになります。

T, T, H, T, T, T, H, T, …

そして、この選択には習慣があります。Chibanyはとんかつが大好きなので、とんかつの日の後は大抵また食べたくなります。ただし、時々気分が変わることもあります。そしてハンバーガーの日の後は、ほぼ必ずとんかつに戻ります。

Chibanyのラボメートのジャマルがコーヒー片手に近づいてきます。

ジャマル：「Chibanyが昨日とんかつを食べたと教えたら、今日を予測できる？」

Chibany：「たぶんまたとんかつかな——だいたい繰り返すから。」

ジャマル：「なるほど。では——先週の火曜日に何を食べたかも知っていると、役に立つ？」

Chibany：(少し間を置いて)「……いや。昨日のことがわかれば、先週の火曜日は何も追加しない。」

別のラボメートのアリッサが顔を上げます。

アリッサ：「今、とても有名な仮定を説明したわよ。」

Chibanyの直感——今日を予測するには昨日だけわかれば十分——がこの章のすべてです。それに名前をつけましょう。

時刻$t$における状態を$X_t$と書きます。これは$t$日目に真であることを表します。Chibanyの場合、$X_t$は二つの値のどちらかで、$X_t \in \{\text{T}, \text{H}\}$（とんかつまたはハンバーガー）です。

命名：マルコフ性

Chibanyの直感には名前があります。未来が現在を与えられた条件のもとで過去から独立しているとき、そのプロセスはマルコフ性を持つと言います。

左から右に読んでください。ゼロ日目まで遡る全履歴を与えられたときの明日の状態の確率は、今日だけを与えられたときの明日の状態の確率と等しい。今日（$X_t$）がわかれば、それ以外の歴史——先週の火曜日、先月——は明日について何も加えません。この性質を持つ状態の系列をマルコフ連鎖と呼びます。

これは第8章と同じ独立性の考え方を、単に時刻でインデックス付けしたものです。あそこでは矢印$A \to B$が「$B$は$A$に依存する」を意味していました。ここでは依存性が時間の矢に沿って走ります。$X_{t} \to X_{t+1}$、そして一歩だけです。明日の親はちょうど一つ——今日だけです。

同じ連鎖の二つの顔

マルコフ連鎖には二つの顔があり、どちらも同じ対象を異なる方法で描いたものです。

図。 各移動の確率をラベル付きの矢印で示す、二つの状態：

graph LR
    T(("T<br/>tonkatsu")) -->|0.65| T
    T -->|0.35| H(("H<br/>hamburger"))
    H -->|0.82| T
    H -->|0.18| H
    classDef node fill:none,stroke:#9bbcff,stroke-width:2px,color:#fff
    class T,H node
    linkStyle default stroke:#9bbcff,stroke-width:2px,color:#fff

Tのループ（0.65）は「とんかつの後は65%の確率でとんかつに留まる」を意味します。矢印T → H（0.35）は時々の気分転換です。HからTへの矢印は強い0.82——ハンバーガーはほとんど繰り返されません。

行列。 同じ四つの数を格子に並べます。推移行列$P$はすべての一歩の確率をまとめており、行$i$、列$j$の要素は次のように定義されます。

状態$i$にいるときに状態$j$へ移動する確率です。Chibanyの場合：

列を下に読むと明日の選択がどこから来るかがわかり、行を横に読むと今日がどこへ向かうかがわかります。行1は「今日がTなら、明日は(T, H)が確率$(0.65, 0.35)$」を意味します。行2は「今日がHなら」を意味します。この表を行列$P$として書きます。行はから、列はへを表し、どちらも(T, H)の順に並べます。

各行の和が1になることに注目してください。$0.65 + 0.35 = 1$、$0.82 + 0.18 = 1$です。これは必然です——どの状態から始めても何かが明日起こるので、可能な次の状態すべての確率の和は1でなければなりません。すべての行が確率分布（非負で和が1）になっている行列を行確率的と言います。この小さな事実は後で多くの静かな仕事をします。$P$の累乗が有効な確率表のままである理由であり、（後で見るように）$P$が常に1に等しい特別な固有値を持つ理由でもあります。

図と行列は同一の情報を持ちます。図は人間にとって見やすく、行列はコンピュータにとって扱いやすく——そして今から見るように、サンプリングにも扱いやすいです。

行列はサンプラーである

連鎖は実際に一歩をどのように踏むのでしょうか？連続的なダイスを転がします。

今日が**とんかつ（T）**だとすると、行1を読みます。$(0.65, 0.35)$です。明日を選ぶには：

1. 0から1の間に一様に分布する乱数$u$を引きます——Uniform(0, 1)のサンプルです。例えば$u = 0.42$としましょう。
2. $u < 0.65$？ はい（$0.42 < 0.65$） → Tに留まる。（$u$が0.65より上に落ちていたらHに切り替え。）

それだけです。一様乱数を一つ引き、一度比較し、一歩進む。明日も着地した行から同じことを繰り返し、その翌日も——行列と乱数の流れが系列全体を生成します。このアイデアを覚えておいてください。次の二章のすべての種であり、第7週のモンテカルロの種でもあります。推移行列は単なる確率の表ではなく、サンプルを引くレシピです。

連鎖を走らせる

Chibanyの連鎖を実際に走らせましょう。以下は、とんかつ（T）から始めて20日間前進させた5回の独立した実行です（これらは本章末のコードから直接出力されます）。

T H H T T H T H H T T T T T T H H T H T T
T H T T T T T H T T T T T H H T T T T H T
T T T H T T T T H T H T H T T H T H T T T
T T T T T H H T T H T T T T T H H T T T H
T T T T T T H T T T H T T T H T H T T H H

二つのことが目に飛び込んできます。まず、連鎖にはちょうど一歩の記憶があります。とんかつは小さな連続の中で持続する傾向があり、ハンバーガーが現れると次の日にはほぼ必ずTに戻ります——0.82の矢印が約束した通りです。次に、そのような局所的な揺れにもかかわらず、安定した全体的な混合が浮かび上がります。文字を数えると、どの実行でも大まかに70%がT、30%がHです。

この長期的な70/30は偶然ではなく、手で設定したものでもありません。連鎖が落ち着く性質です。それを正確に表現することがこの章の核心です。

定常分布

Chibanyより前に、トランプのデッキをシャッフルすることを考えてみましょう——これもマルコフ連鎖です。状態はデッキの現在の並び順で、推移は「ランダムなカードを取って一番上に移す」などの操作を繰り返すことです。

ジャマル：「ゲームの前にシャッフルするとき、実際に何を達成しようとしているの？」

特定の並び順ではありません——特定のゴール配置に向けてシャッフルしているわけではありません。望むのはすべての並び順が等しく可能になることで、誰もデッキを予測できなくなることです。そして微妙な点があります。「すべての並び順が等しく可能」に達した後は、もう一度シャッフルしてもそこに留まります。よくシャッフルされたデッキを「さらにシャッフル」することはできません。全並び順上の一様分布は自己持続的です。プロセスはそれを変えなくなります。

この自己持続的で、もはや変化しなくなった分布には名前があります。まずは記号から。任意のマルコフ連鎖を非常に長い時間走らせ、問いかけます。各状態に滞在する時間の割合は？ その割合のリストを$\pi$（ギリシャ文字の「パイ」）と呼びます——状態ごとに一つの要素を持ち、和が1になる行ベクトルです。Chibanyの場合、$\pi = (\pi_T, \pi_H)$で、$\pi_T$はとんかつの日の長期的な割合、$\pi_H$はハンバーガーの日の割合です。

この長期的な分布$\pi$は、一歩進めても変化しない唯一の分布である定常分布と一致します。「分布を一歩前進させること」を$P$による乗算として書くと（次のセクションで正確にします）、その性質は

言葉にすると、今日の状態についての信念がすでに$\pi$であれば、もう一歩後の明日の状態についての信念は依然として$\pi$である——永遠に。Chibanyの場合、$\pi$は非常に具体的な問いに答えます。一学期全体で、Chibanyのランチのうちとんかつはどのくらいの割合？ サンプル実行はすでに答えを囁いていました——約70%です。

二状態なら線形方程式を手で解けます。52!通りのカードの並び順ではお手上げです。スケールする方法が必要です——そしてすでに持っています。ただ連鎖を走らせればよいのです。

πを求める：ただ走らせればよい（べき乗反復）

ここに仕掛けがあります。任意の現在の状態についての信念から始めます——確率の行ベクトル$\mathbf{v} = (v_T, v_H)$として書きます——そして$P$を乗算して一歩前進させます。

なぜ$P$を乗算することで「一日前進する」のでしょうか？積を書き出せばわかります。明日のとんかつの確率は、今日Tにいて留まる確率と、今日Hにいて戻ってくる確率の和です。

これは昨日の状態にわたる全確率の法則そのものです——手で行う一歩の予測と同じです。したがって$\mathbf{v}P$は一日後の信念、$\mathbf{v}P^2 = (\mathbf{v}P)P$は二日後、$\mathbf{v}P^3$は三日後、などとなります（指数は日数を数えるだけです）。

乗算を続けると、系列はどこから始めても定常分布$\pi$に収束します。（収束は証明しませんが、今まさに目撃しようとしています。）この手続き——行列を繰り返し乗算すること——をべき乗反復と呼びます。

Chibanyの連鎖でこれを見てみましょう。二つの方法で実行します。とんかつの絶対確信から$\mathbf{v} = (1, 0)$、そしてハンバーガーの絶対確信から$\mathbf{v} = (0, 1)$。

ステップ0では二つのスタートがこれ以上ないほど異なっています。一歩後にはすでに差が縮まっています。ステップ5では同一です——どちらも$(0.701, 0.299)$——そこに留まります。連鎖はどこから始まったかを忘れました。

スタートが関係ない理由：エルゴード性

この「忘却」は魔法ではありません。条件があります。連鎖が混ざる——すべての出発点から同じ$\pi$に収束する——のは、エルゴード的なときです。つまり、どの状態からどの状態へも（複数ステップを経ても）到達でき、連鎖が固定サイクルの中でトラップされていないときです。Chibanyの連鎖はエルゴード的です——TからHへ、またその逆へ到達でき、厳密な周期がありません——したがって混ざり、定常分布は始まりとは無関係に連鎖の性質です。（エルゴード的でない連鎖——例えば二状態間に移動手段がない——は、スタート地点によって異なる長期分布に留まることがあります。その条件が排除するのはそのような例外的なケースです。）

πとは何か：固有値1の固有ベクトル

べき乗反復は$\pi$に収束します。しかし$\pi$は数学的対象として正確に何なのでしょうか？

定義方程式$\pi P = \pi$をもう一度見てください。これは、行ベクトル$\pi$を行列$P$で乗算すると同じベクトルが変わらずに戻ってくると言っています。行列が同じ方向を向いたままにするベクトルには名前があります。固有ベクトルです。行列の固有ベクトルとは、その行列を乗算したとき方向が変わらないベクトルです——長さだけが固有値と呼ばれる数でスケールされます。$\pi P = \pi$は$\pi$を完全に変えないので（正確に1でスケールされる）、$\pi$は固有値1の$P$の固有ベクトルです。（左側で乗算するので、$\pi P = \pi$は技術的には左固有ベクトルです——細部ですが、考え方は同じです。）

すべての行確率的行列はそのような固有ベクトルを持ちます——常に1に等しい固有値が存在します。したがって$\pi$を求める等価な二つの方法があります。連鎖を走らせる（べき乗反復、上で行ったこと）か、線形代数ルーチンを使って固有値方程式を直接解くかです。先ほど見た収束の表はべき乗反復がその固有ベクトルに忍び寄っていたものです。$\times P$のもとで数値が変化しなくなったとき、それを見つけたことになります。

線形代数についてはここまでで十分です。固有値1の固有ベクトルが定常分布であり、この章で行う他のすべてのことは、単に連鎖を走らせることで実行できます。

少し難しいもの：三状態ウォーク

これらのどれかが状態が二つしかないことに依存しているのでしょうか？まったくそうではありません——方法は状態数に関係なく同じです。以下はより複雑な推移行列を持つ三状態連鎖です（状態1、2、3）。

graph LR
    S1(("1")) -->|0.1| S2(("2"))
    S1 -->|0.9| S3(("3"))
    S2 -->|0.5| S1
    S2 -->|0.5| S3
    S3 -->|0.8| S1
    S3 -->|0.2| S2
    classDef node fill:none,stroke:#9bbcff,stroke-width:2px,color:#fff
    class S1,S2,S3 node
    linkStyle default stroke:#9bbcff,stroke-width:2px,color:#fff

行を前と同じように読みます。行1は「状態1から、確率0.1で状態2へ、確率0.9で状態3へ行く」を意味します。（対角はすべてゼロ——この連鎖には自己ループがなく、常に移動します。）

計算する前に、推測してみましょう。このウォークを長時間走らせると——どの状態が最も多く訪問されるでしょうか？矢印をたどります。状態1はトラフィックの0.9を状態3に送り、状態2は1と3に分割し、状態3は0.8を1に戻します。したがって状態1と3は人気の目的地ですが、状態2は矢印のターゲットになることが少ない（状態1から0.1、状態3から0.2のみ）。状態2が最も孤独なことが予想されます。

べき乗反復を実行すると——同じ手順を二つの異なるスタートから——収束先は

まさに推測通りです。状態1と3がほとんどの訪問を吸収し（$0.42$と$0.45$）、状態2はめったに訪問されません（$0.13$）。そしてChibanyの場合と同様、答えはどこから始めるかに依存しません——状態1か状態2から始めて十分に走らせると、同じ$\pi$に着地します。この連鎖もエルゴード的だからです。

二状態、三状態、または52!通りのカードの並び順——同じ方法、同じ種類の答え。連鎖を走らせ（または固有値1の固有ベクトルを解き）、長期的な居場所を読み取ります。

GenJAXとJAXの実装

この章は二種類の計算に明確に分かれており、どのツールがどちらに合っているかについて正直であることが重要です。

• 機構の大部分——$P$の構築、べき乗反復、固有ベクトル——は単純な線形代数です。ベクトルを行列で繰り返し乗算します。jax.numpyがまさに適切なツールであり、生成モデルを持ち出すと却って分かりにくくなります。
• 真に生成的な部分——選択の系列をサンプリングすること、「行列はサンプラーである」——は、第8章〜10章のベイズネットと同様に、GenJAXの@genモデルが真価を発揮するところです。

それぞれに合った場所で使います。

べき乗反復（純粋なJAX）

まず機構から。$P$を構築し、小さなpower_iterateを作って$P$を繰り返し乗算し、二つの異なるスタートが同じ70/30に収束する様子を観察します。

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20

import jax.numpy as jnp

# Chibany's transition matrix. Rows = today (T, H); columns = tomorrow (T, H).
P = jnp.array([[0.65, 0.35],
               [0.82, 0.18]])

def power_iterate(v, P, steps):
    """Multiply the distribution v by P, `steps` times, recording every iterate."""
    rows = [v]
    for _ in range(steps):          # steps is a plain Python int — an ordinary loop
        v = v @ P                   # @ is matrix multiply: step the distribution forward one day
        rows.append(v)
    return jnp.stack(rows)

traj_T = power_iterate(jnp.array([1.0, 0.0]), P, 20)   # always start on tonkatsu
traj_H = power_iterate(jnp.array([0.0, 1.0]), P, 20)   # always start on hamburger

for k in [0, 1, 2, 5, 20]:
    print(f"step {k:2d}: from T ({traj_T[k,0]:.3f}, {traj_T[k,1]:.3f})   "
          f"from H ({traj_H[k,0]:.3f}, {traj_H[k,1]:.3f})")

出力：

step  0: from T (1.000, 0.000)   from H (0.000, 1.000)
step  1: from T (0.650, 0.350)   from H (0.820, 0.180)
step  2: from T (0.709, 0.290)   from H (0.681, 0.319)
step  5: from T (0.701, 0.299)   from H (0.701, 0.299)
step 20: from T (0.701, 0.299)   from H (0.701, 0.299)

ステップ5で二つのスタートは区別できなくなります——連鎖は始まりを忘れ、$\pi \approx (0.70, 0.30)$に落ち着きました。

πは固有値1の固有ベクトル（純粋なNumPy）

べき乗反復は固有ベクトルに忍び寄ります。線形代数ルーチンに直接求めて、二つが一致することを確認することもできます。$P$の左固有ベクトルは転置行列$P^{\mathsf T}$の固有ベクトルです。

 1
 2
 3
 4
 5
 6
 7
 8
 9
10

import numpy as np

vals, vecs = np.linalg.eig(np.array(P).T)   # left eigenvectors of P = eigvecs of P^T
idx = int(np.argmin(np.abs(vals - 1.0)))    # the eigenvalue closest to 1 (eig returns floats)
pi = np.real(vecs[:, idx])
pi = pi / pi.sum()                          # normalize so it sums to 1

print(f"eigenvalue   = {np.real(vals[idx]):.3f}")
print(f"pi (eigen)   = ({pi[0]:.3f}, {pi[1]:.3f})")
print(f"pi (run it)  = ({traj_T[20,0]:.3f}, {traj_T[20,1]:.3f})")

出力：

eigenvalue   = 1.000
pi (eigen)   = (0.701, 0.299)
pi (run it)  = (0.701, 0.299)

固有値は正確に1（行確率的の保証通り）であり、固有ベクトルはべき乗反復が見つけたものと一致します。二つの経路、一つの答え。

行列はサンプラーである（GenJAX）

次は生成的な部分です。推移行列は単なる表ではありません——サンプルを引きます。GenJAXの@genモデルを書いて、$P$の現在の行から次の状態をサンプリングして一歩踏みます。GenJAXのcategorical分布は対数確率からカテゴリを選ぶので、jnp.log(P)[state]——現在の行の対数——を渡します。

系列全体をサンプリングするために、固定ステップ数のモデルを構築します。日数を受け取り、ちょうどその日数の生成モデルを返す小さなファクトリ関数（make_chain）を使います。

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26

import jax
import jax.random as jr
from genjax import gen, categorical

LOGP = jnp.log(P)   # categorical wants log-probabilities; row `state` = current state

def make_chain(n_steps):
    """Build a generative model that samples `n_steps` transitions of Chibany's chain."""
    @gen
    def chain(start):
        state = start
        states = [state]
        for t in range(n_steps):                  # n_steps is a captured Python int
            state = categorical(LOGP[state]) @ f"x_{t}"
            states.append(state)
        return jnp.array(states)
    return chain

chain20 = make_chain(20)

# Five independent runs, each starting on tonkatsu (state 0). vmap runs them in parallel.
labels = ["T", "H"]   # plain Python list — JAX arrays can't hold strings
keys = jr.split(jr.key(0), 5)
runs = jax.vmap(lambda k: chain20.simulate(k, (0,)).get_retval())(keys)
for r in runs:
    print(" ".join(labels[int(s)] for s in r))

出力：

T H H T T H T H H T T T T T T H H T H T T
T H T T T T T H T T T T T H H T T T T H T
T T T H T T T T H T H T H T T H T H T T T
T T T T T H H T T H T T T T T H H T T T H
T T T T T T H T T T H T T T H T H T T H H

これらは「連鎖を走らせる」のセクションの系列そのものです——作り上げたのではなく、生成したものです。それぞれがChibanyの連鎖からの本物のサンプルです。

線形代数ではなくサンプリングによって長期的な70/30を確認するために、一つの長い連鎖を走らせてカウントします。数千ステップの長い実行にはjax.lax.scanを使います。JAXが値を多くのステップを通じて効率的に引き継ぐ方法です。そのstep関数を「現在のstateと新鮮な乱数キーを与えられ、次の状態を二回返す」と読んでください。最初のコピーは次のステップに渡されるキャリー、二番目はすべての訪問した状態の出力配列に収集される値です。

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11

def run_long(key, start, n):
    """Sample n steps by carrying the current state through a scan."""
    def step(state, k):
        nxt = jr.categorical(k, LOGP[state])   # one transition from the current row
        return nxt, nxt
    _, visited = jax.lax.scan(step, start, jr.split(key, n))
    return visited

visited = run_long(jr.key(1), 0, 5000)
frac_tonkatsu = float(jnp.mean((visited == 0).astype(float)))
print(f"long-run fraction tonkatsu (5000 steps) ~ {frac_tonkatsu:.2f}")

出力：

long-run fraction tonkatsu (5000 steps) ~ 0.71

サンプリングと線形代数が一致します。連鎖の長い実行はその日々の約70%をとんかつに費やします——連鎖を走らせることで見つけた定常分布です。

三状態連鎖

最後に、三状態行列$A$に同じpower_iterateを適用し、二つの異なるスタートから$\pi \approx (0.42, 0.13, 0.45)$を再現します。

1
2
3
4
5
6
7
8

A = jnp.array([[0.0, 0.1, 0.9],
               [0.5, 0.0, 0.5],
               [0.8, 0.2, 0.0]])

from_1 = power_iterate(jnp.array([1.0, 0.0, 0.0]), A, 40)
from_2 = power_iterate(jnp.array([0.0, 1.0, 0.0]), A, 40)
print("from state 1:", tuple(round(float(x), 2) for x in from_1[40]))
print("from state 2:", tuple(round(float(x), 2) for x in from_2[40]))

出力：

from state 1: (0.42, 0.13, 0.45)
from state 2: (0.42, 0.13, 0.45)

二つのスタート、一つの答え——状態2が最も孤独で、矢印が予測した通りです。機構は方法の一行も変えることなく二状態から三状態へとスケールしました。

練習問題

付属のノートブックがこれらすべてを対話的に解説しています。

📓 Colabで開く：13_markov_chains.ipynb

このチュートリアルシリーズへの寛大なご支援にJPCCAに特別の感謝を申し上げます。

参考文献

• Assaf, S., Diaconis, P., & Soundararajan, K. (2011). A rule of thumb for riffle shuffling. The Annals of Applied Probability, 21(3), 843–875. https://doi.org/10.1214/10-AAP701
• Bayer, D., & Diaconis, P. (1992). Trailing the dovetail shuffle to its lair. The Annals of Applied Probability, 2(2), 294–313. https://doi.org/10.1214/aoap/1177005705
• Berger, P. D. (1973). On the distribution of hand patterns in bridge: Man-dealt versus computer-dealt. The Canadian Journal of Statistics, 1(2), 261–266.
• Conger, M., & Viswanath, D. (2006). Riffle shuffles of decks with repeated cards. The Annals of Probability, 34(2), 804–819. https://doi.org/10.1214/009117905000000675
• Diaconis, P. (2003). Mathematical developments from the analysis of riffle shuffling. In A. A. Ivanov, M. W. Liebeck, & J. Saxl (Eds.), Groups, combinatorics and geometry (pp. 73–97). World Scientific.
• Thorp, E. O. (1973). Nonrandom shuffling with applications to the game of Faro. Journal of the American Statistical Association, 68(344), 842–847. https://doi.org/10.1080/01621459.1973.10481434