「動くか」ではなく「戻るか」:マルコフ連鎖で状態を分類する【第22回】

はじめに

第14章はマルコフ連鎖です。今回から確率過程編に入ります。

前回(第21回・標本調査法)まで、線形モデル編では「手元にある固定したデータ」を扱ってきました。今回からは対象が変わります。時間とともに状態が移っていくものが相手です。

実は筆者は、この章の内容を一部すでに使ったことがあります。学生時代の研究で隠れマルコフモデル(HMM)を扱っていたので、マルコフ性や推移確率行列は馴染みがありました。定常分布の手計算も前日に練習して解けるようになっていました。

なので「復習で済むだろう」と思って、確認のために小さなドリルを6問解いてみたのです。結果は4か所で間違えました。しかも間違えたのは計算ではなく、全部用語の定義でした。

いちばん象徴的だったのがこれです。

「そこから動かないから一時的(transient)だろう」と考えたら、完全に逆だった。

動かない状態こそが再帰的で、活発に行き来している状態が一時的でした。字面から受ける印象と定義が食い違っていたわけです。

この記事は、その間違いをそのまま載せる形で書いています。ドリルの問題も全部載せるので、読みながら解いてみてください。同じところで転ぶ人は多いはずだと思っています。

なお筆者は統計の専門家ではありません。理解の誤りが含まれる可能性があるため、試験対策として読む場合は必ず公式テキストで確認してください。

この記事で使う言葉

先に用語をまとめます。略語も先に開いておきます。

用語読み方・スペルアウト意味
マルコフ性Markov property次に何が起きるかが、いまの状態だけで決まる(過去の経路を忘れる)
推移確率行列transition probability matrix状態 ii から状態 jj へ1歩で移る確率を並べた行列 PP
定常分布stationary distributionπP=π\pi P = \pi を満たす分布。方程式の解
極限分布limiting distributionPnP^nnn \to \infty で実際に収束した先
既約irreducibleどの状態からどの状態へも到達できる
周期periodその状態に戻れる歩数の最大公約数
再帰的recurrent(再発する)出ていったら確率1で戻ってくる
一時的transient(一過性の)戻らない可能性がある
吸収状態absorbing state入ったら出られない状態
C-K方程式Chapman-Kolmogorov equationnn 歩の推移を途中で分解する式
HMMHidden Markov Model=隠れマルコフモデル状態が見えず、そこからの出力だけが見える構造
M/M/1Markovian arrival / Markovian service / 1 server到着と処理がともに指数的で窓口が1つの待ち行列

何に使うのか

定義より先に、この章が答える問いを並べます。マルコフ連鎖で聞かれることは、実は3種類しかありません

第一に、放っておいたら最終的にどうなるか。これが定常分布です。Webサイトの回遊、在庫の水準、検索順位、サーバーの混雑度がここに入ります。

第二に、そこに着くまでどれくらいかかるか。平均到達時間や平均回帰時間です。

第三に、そもそも着けるのか、着けるとしたらどの確率で。これが吸収確率です。

そして重要なのは、第一と第三が別の道具だということです。定常分布は「いつまでも続く鎖」の話で、吸収確率は「いつか終わる鎖」の話です。混ぜると解けません(後で実際に解けなくなる例を出します)。

知りたいこと道具鎖の性質
長期的な状態の割合定常分布 πP=π\pi P = \piずっと続く
戻ってくるまでの平均時間1/πi1/\pi_iずっと続く
どちらの結末になるか吸収確率 hi=jPijhjh_i = \sum_j P_{ij}h_jいつか終わる
終わるまでの平均時間ki=1+jPijkjk_i = 1 + \sum_j P_{ij}k_jいつか終わる

マルコフ性は条件付き確率の特殊な形

まず定義です。第2回の条件付き確率をそのまま使います。

P(Xn+1=jXn=i, Xn1=in1, , X0=i0)=P(Xn+1=jXn=i)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)

左辺は「これまでの全履歴を知ったうえでの確率」、右辺は「いまの状態だけを知ったうえでの確率」です。この2つが等しいというのがマルコフ性です。

言い換えると、過去の情報は現在の状態に要約しつくされているということです。どういう経路でいまここに来たかは、次に何が起きるかに影響しない。

第5回の幾何分布の無記憶性が近い発想です。実際、ある状態に留まり続ける回数は幾何分布になります。状態 ii に留まる確率が PiiP_{ii} なら、そこに kk 回連続で留まる確率は Piik1(1Pii)P_{ii}^{k-1}(1-P_{ii}) で、これは幾何分布そのものです。

この右辺を並べたものが推移確率行列 PP です。PijP_{ij} が「状態 ii から jj へ1歩で移る確率」で、各行の和が1になります(どこかへは必ず行くので)。

ここで1つ注意です。行の和が1なのですが、定常分布を求めるときは列を縦に読みます。筆者はここで一度間違えました。πP=π\pi P = \pi を成分で書くと

πj=iπiPij\pi_j = \sum_i \pi_i P_{ij}

で、右辺は jj 列目を縦に足しています。「その状態に流れ込んでくる量の合計」が「その状態の量」に等しい、というつりあいの式だからです。

チャップマン・コルモゴロフの方程式は行列の掛け算なのか

名前が大げさなので、正体を確認しておきます。nn 歩での推移確率を Pij(n)P^{(n)}_{ij} と書くと、

Pij(m+n)=kPik(m)Pkj(n)P^{(m+n)}_{ij} = \sum_k P^{(m)}_{ik} P^{(n)}_{kj}

これがC-K方程式です。中身は「途中の時点で必ずどこかの状態にいる」という全確率の分解にすぎません。mm 歩目にどこにいたかで場合分けして足しているだけです。

そして結論から言えば、有限個の状態を扱う限り、これは行列の積そのものです。 実際に確かめました。

P=(0.70.20.10.30.40.30.20.30.5)P = \begin{pmatrix} 0.7 & 0.2 & 0.1 \\ 0.3 & 0.4 & 0.3 \\ 0.2 & 0.3 & 0.5 \end{pmatrix}

について kPik(2)Pkj(3)\sum_k P^{(2)}_{ik}P^{(3)}_{kj} を素朴な二重ループで計算し、P5P^5 を行列べきで計算して比べると、最大の差は 1.1×10161.1 \times 10^{-16} でした。浮動小数点の誤差だけで、完全に同一です。

では、なぜわざわざ名前が付いているのか。理由は2つあります。

1つは、行列で書けない場合にも成り立つことです。状態が連続的なとき(次回のブラウン運動など)、和は積分になります。

pm+n(x,y)=pm(x,z)pn(z,y)dzp_{m+n}(x,y) = \int p_m(x,z)\,p_n(z,y)\,dz

行列の積という表現はもう使えませんが、C-K方程式としては同じ形です。行列積はC-K方程式の特殊ケースという関係です。

もう1つは、この式がマルコフ性から導かれるという点です。マルコフ性がなければ、途中の時点で場合分けしたときに「そこに至った経路」が効いてしまい、単純な積に分解できません。

試験では「PnP^n を求めよ」「nn 歩後の分布を求めよ」の形で出るので、実質は行列べきの計算だと思って構いません。名前に怯む必要はありませんでした。

状態の分類:4つの用語は判定しているものが違う

ここが今回の本題です。既約・周期的・再帰的・一時的という4つの用語が出てきます。筆者はこれで混乱しました。

混乱の原因がわかったので先に書きます。この4つは判定している対象のレベルが違います。

マルコフ連鎖の状態の分類を4つのパネルで対比した状態遷移図。(a)既約は3つの状態A・B・Cが一周できる形で、A→B→C→Aの矢印と各状態の自己ループを持つ。(b)既約でないは{A,B}と{C,D}という2つの島に分かれ、島の間を結ぶ矢印がない。(c)周期3はA→B→C→Aだけが確率1で繋がり自己ループがない。(d)一時的と再帰的では「試用」「有料」から「解約」へ一方向に進み、解約だけが確率1の自己ループを持つ。

用語判定の対象判定しているもの関係するもの
既約鎖全体どの状態からどの状態へも到達できるか閉じたかたまりが2つ以上になると定常分布が一意でなくなる
周期状態(既約なら全体で共通)戻れる歩数のリズム極限分布が存在しなくなる
再帰的 / 一時的個々の状態出ていったら戻ってくるか

壊れるものが違うというのが要点です。定常分布の一意性が壊れるか、極限分布の存在が壊れるか。この2つは別の話なので、混ぜないようにします。

1点だけ先に注意しておきます。既約性は一意性の十分条件であって必要条件ではありません。 既約でなくても、閉じたかたまりが1つなら定常分布は一意です(後の節で確かめます)。筆者はここを誤解していました。

覚え方としては、既約は「地図が繋がっているか」、周期は「時計の刻みが揃っているか」、再帰/一時は「行ったら帰れるか」です。

ドリル1:この鎖を分類してください

いきなり問題です。会員の状態を4つに分けた鎖を考えます。数字は1か月後の推移確率です。

現在\翌月新規継続休眠退会
新規00.70.20.1
継続00.80.150.05
休眠00.30.60.1
退会0001

(a) この鎖は既約でしょうか。 (b) 4つの状態それぞれについて、再帰的か一時的かを判定してください。 (c) 定常分布を計算せずに、それが何になるか予想できますか。

筆者の答えはこうでした。

(a) 既約です。一つの塊なので。 (b) 新規と退会はそこから動かないので一時的、残りの二つは再帰的かな (c) 最終的には退会で安定するね

(c) だけが正解で、(a)(b) は間違いです。 順に見ていきます。

既約は「往復できるか」を聞いている

(a) この鎖は既約ではありません。

「一つの塊」という感覚はわかります。バラバラに分断されているわけではないからです。しかし既約の定義はどの状態からどの状態へも到達できることで、往復が要求されます

退会から継続へ戻る矢印はありません。だから「退会 → 継続」が到達不可能で、既約ではない。表の最後の行が 0,0,0,10,0,0,1 になっている時点で決まります。

図の (a) と (b) を見比べてください。(a) は ABCAA \to B \to C \to A と一周できるので既約です。(b) は {A,B}\{A,B\}{C,D}\{C,D\} の間に矢印がないので既約ではありません。一周できる経路があるかで見るのが実用的です。

再帰的と一時的は「動くか」ではなく「戻るか」

(b) が完全に逆でした。 ここが今回いちばんの学びです。

筆者は「そこから動かないので一時的」と考えました。しかし一時的(transient)は「出ていったら戻ってこないかもしれない」という意味で、動かないことではありません。判定基準は帰還確率だけです。

正しく判定するとこうなります。

状態出た後、自分に戻ってくるか判定
新規新規に入る矢印が誰からも無い(列が全部0) → 帰還確率 0一時的
継続退会に抜ける道があるので、いつか必ず抜ける → 帰還確率 <1< 1一時的
休眠同じく退会に抜ける → 帰還確率 <1< 1一時的
退会出られない → 帰還確率 1再帰的(吸収状態)

退会だけが再帰的で、残り3つは全部一時的です。「そこから動かない」退会こそが再帰的で、活発に行き来している継続・休眠が一時的、という逆転が起きます。

なぜ字面に騙されるのか。英語を見ると納得できます。

  • recurrent = 再発する(病気が再発する、と同じ語)
  • transient = 一過性の(一時的な症状、と同じ語)

つまり動くかどうかではなく、繰り返し訪れるか、通り過ぎるだけかを見ています。継続と休眠は「何度も行き来するけれど、最後には必ず抜けていく」ので、通り過ぎるだけの状態です。

そして決定的なことに気づきました。(c) の答えがそのまま (b) の答えになっています。

(c) 「最終的には退会で安定する」は正解で、定常分布は π=(0, 0, 0, 1)\pi = (0,\ 0,\ 0,\ 1) です。これは計算せずに構造から読めます。

ということは、新規・継続・休眠の確率がゼロになるということです。それが「一時的」の意味です。(0,0,0,1)(0,0,0,1) の 0 が並んでいる位置が一時的な状態、1 になっている位置が再帰的な状態。

自分で正しく予想した結論の中に、間違えた問いの答えが入っていたわけです。用語だけが接続されていなかった、ということでした。

訪問回数の期待値が有限か無限か

もう少し厳密な判定基準を書いておきます。試験で使える形です。

一時的な状態は、その状態を訪問する回数の合計が有限です。再帰的なら無限です。「必ず戻ってくるなら無限回訪問する」と考えれば当たり前ですね。

これは計算できます。図 (d) の「試用 → 有料 → 解約」で確かめます。

P=(0.70.20.100.950.05001)P = \begin{pmatrix} 0.7 & 0.2 & 0.1 \\ 0 & 0.95 & 0.05 \\ 0 & 0 & 1 \end{pmatrix}

一時的な状態(試用・有料)だけを取り出した 2×22\times2 の行列を QQ として、

N=(IQ)1=(3.33313.333020)N = (I - Q)^{-1} = \begin{pmatrix} 3.333 & 13.333 \\ 0 & 20 \end{pmatrix}

この NN基本行列と呼びます。NijN_{ij} が「ii から始めたとき jj を訪問する回数の期待値」です。試用から始めると、試用には 3.3333.333 回、有料には 13.33313.333 回滞在します。どちらも有限なので一時的です。

行の合計 3.333+13.333=16.6673.333 + 13.333 = 16.667吸収されるまでの平均ステップ数です。10万回のシミュレーションで 16.7516.75、40万回で 16.698±0.03016.698 \pm 0.030(理論値からのずれは標準誤差の1.04倍)でした。一致しています。

ちなみに「試用から一度でも有料に到達する確率」も同じ考えで出せます。f=0.2+0.7ff = 0.2 + 0.7f を解いて f=0.2/0.3=0.6667f = 0.2/0.3 = 0.6667。シミュレーションで 0.66650.6665 でした。1歩先で場合分けして、自分自身に戻る項を右辺に残すというのがコツで、これが後で吸収確率の型として再登場します。

無限の鎖では直感が効かない:ランダムウォーク

有限の鎖なら「抜け道があるか」を見るだけなので簡単です。厄介なのは状態が無限にある場合です。

ランダムウォークの経路を6本ずつ描いた2枚のグラフ。左はp=0.5の対称な場合で、6本の経路がいずれも0の横線を何度も上下に横切っている。右はp=0.6で上に偏った場合で、経路が右肩上がりに離れていき0の線に戻らないものがある。縦軸は原点からの位置、横軸は歩数で2000歩まで。

直線上のランダムウォークを考えます。確率 pp+1+1、確率 q=1pq = 1-p1-1 動きます。原点は再帰的でしょうか。

p=0.5p = 0.5(左)のとき、原点は再帰的です。確率1で必ず戻ってきます。 ただし、ここに罠があります。

戻るまでの平均時間は無限なのです。確率1で戻るのに、平均待ち時間が無限。この状態を零再帰的と呼びます。

数値で見ると納得できます。2n2n 歩までに一度も原点に戻らない確率は (2nn)/4n\binom{2n}{n}/4^n で、これは 1/πn1/\sqrt{\pi n} に近づきます。つまり戻る確率は

歩数の上限戻る確率
10歩0.7539
100歩0.9204
1,000歩0.9748
10,000歩0.9920
100,000歩0.9975
1億歩0.99992

1に収束しますが、速さが 1/n1/\sqrt{n} という遅さです。第8回n\sqrt{n} の壁がここにも顔を出しています。この遅さのせいで平均が無限に発散します。

p=0.6p = 0.6(右)のとき、原点は一時的です。 上に流されて二度と戻らない経路が存在します。帰還確率は 0.80.8 で、2割の経路は永久に帰ってきません。10万歩以内に戻った割合は 0.80330.8033 でした。

なお、この帰還確率の計算で筆者は間違えました。詳細は最後の「間違えていたこと」に書きます。

再帰的の中の区別をまとめておきます。試験で用語として問われる部分です。

分類帰還確率平均帰還時間定常分布
正再帰的1有限作れる有限の既約な鎖はすべてこれ
零再帰的1無限作れないp=0.5p=0.5 のランダムウォーク
一時的<1< 1無限(正の確率で戻らないので発散)p0.5p \ne 0.5 のランダムウォーク

有限の状態数しかない既約な鎖は、必ず正再帰的になります。だから零再帰という概念が必要になるのは、状態が無限にある場合だけです。実務では出てこないので、用語として知っていれば十分でしょう。

周期は最小公倍数ではなく最大公約数

次の用語です。ここも間違えました。

周期の定義:その状態に戻れる歩数の最大公約数。 1なら非周期的と呼びます。

ドリル2:周期を求めてください

矢印だけ書きます(確率は等分と考えてください)。状態0の周期はいくつでしょうか。

  • A: 01200 \to 1 \to 2 \to 0 のみ
  • B: 010 \to 1020 \to 2101 \to 0202 \to 0
  • C: 012300 \to 1 \to 2 \to 3 \to 0 に加えて 020 \to 2 という近道がある

筆者の答えは A: 3、B: 2、C: 3.5 でした。A と B は正解、C が間違いです。

周期は整数にしかなりません。 3.5 という答えは平均を取っています。定義は最大公約数です。

C で状態0に戻れる歩数を全部並べると、3,4,6,7,8,9,10,11,3, 4, 6, 7, 8, 9, 10, 11, \dots です。02300 \to 2 \to 3 \to 0 の3歩と、012300 \to 1 \to 2 \to 3 \to 0 の4歩があります。gcd(3,4)=1\gcd(3, 4) = 1 なので周期1、つまり非周期的です。

直感的にはこうです。3歩と4歩の両方が使えるなら、組み合わせて 6=3+36 = 3+37=3+47 = 3+48=4+48 = 4+4991010、… とほぼ全部作れてしまいます(5だけ作れませんが、それ以降は全部)。リズムが揃わなくなるので振動が消えるわけです。

実際 P60P^{60} の1行目は (0.2859,0.1428,0.2855,0.2857)(0.2859, 0.1428, 0.2855, 0.2857) で、定常分布 (2/7, 1/7, 2/7, 2/7)=(0.2857,0.1429,0.2857,0.2857)(2/7,\ 1/7,\ 2/7,\ 2/7) = (0.2857, 0.1429, 0.2857, 0.2857) に収束していました。

近道が1本あるだけで周期が壊れるというのが実務的な要点です。現実の鎖はほぼ非周期になります。

戻れる歩数が16と24なら、周期は48ではなく8

ここで疑問が出ました。16と24なら、最小公倍数の48になるのでは?

答えは8です。gcd(16,24)=8\gcd(16,24)=8 で、最小公倍数の48ではありません。 最大公約数なので、大きくならず小さくなります。

そして、ここで直感を一段精密にする必要がありました。「最大公約数が1なら何でも作れる」という理解は正しいのですが、逆向きは成り立ちません。

長さ16の輪と長さ24の輪が状態0を共有する鎖(39状態)を実際に作って、状態0に戻れる歩数を列挙してみました。

16, 24, 32, 40, 48, 56, 64, 72, 80, 16,\ 24,\ 32,\ 40,\ 48,\ 56,\ 64,\ 72,\ 80,\ \dots

すべて8の倍数ですが、8そのものには戻れません。 16と24をどう足しても8は作れないからです。P8[0,0]=0P^8[0,0] = 0 で、P16[0,0]=0.5P^{16}[0,0] = 0.5 から初めて値が立ちました。

つまり「周期が8」の意味はこうです。

「8歩ごとに戻れる」ではなく、「戻れるのは8の倍数の歩数だけ」。

定義が最大公約数になっている理由は、振動の周期がそこで決まるからです。戻れる歩数が全部8の倍数なら、確率は8ステップおきのタイミングでしか立たず、PnP^n は8つの位相を巡回します。

いくつか並べます。最小公倍数と混同しないよう、対比のために併記します。

戻れる歩数周期(gcd\gcd参考:最小公倍数
3 と 41(非周期)12
2 と 31(非周期)6
4 と 6212
6 と 9318
10 と 15530
16 と 24848
16 と 481648

互いに素な長さの経路が2本あれば必ず非周期というのが実用的な判定法です。特に自己ループ(1歩で戻る)が1本でもあれば、gcd(1,何でも)=1\gcd(1, \text{何でも}) = 1 なので即座に非周期になります。

試験で使える判定:既約な鎖なら、推移確率行列 PP の対角成分に正の値が1つでもあれば非周期。

「既約な鎖なら」という条件が必要です。既約な鎖では全状態が同じ周期を共有するので、1箇所で自己ループがあれば全体が非周期になります。既約でないと成り立ちません。反例を挙げます。

P=(100001010)P = \begin{pmatrix} 1 & 0 & 0 \\ 0 & 0 & 1 \\ 0 & 1 & 0 \end{pmatrix}

対角に P11=1P_{11}=1 がありますが、状態2と状態3は互いに往復するだけなので周期2で、P10P11P^{10} \ne P^{11} のまま振動し続けます。自己ループがある状態と、周期を持つ状態が別のかたまりにいるとこうなります。

図 (c) の周期3の鎖に自己ループを1本足すと、固有値の絶対値が (1, 0.707, 0.707)(1,\ 0.707,\ 0.707) になり、P50P^{50} はきれいに (0.5, 0.25, 0.25)(0.5,\ 0.25,\ 0.25) へ収束しました。自己ループ1本で振動が消えるわけです。

おまけですが、周期8の鎖でも8ステップまとめて1歩と見るP8P^8 を推移行列とする)と振動が消えます。ただし P8P^8 は既約でなくなって8つの閉じたかたまりに分かれるので、正確には「各かたまりの中では非周期になって収束する(行き先は出発点で決まる)」という形です。周期性は「見る間隔が合っていない」だけの現象だと見ることもできます。

定常分布と極限分布は同じものではない

筆者はこの2つを区別せずに理解していました。区別が必要な例を作ります。

左右2枚の折れ線グラフ。縦軸は状態Aにいる確率、横軸はステップ数nで0から20。左の周期的な鎖では青い線(初期分布(1,0,0))が1と0を、赤い線(初期分布(0,1,0))が0と1を交互に往復し、定常分布0.500を示す黒い破線には収束しない。右の非周期な鎖では青と赤の線がどちらも数ステップで黒い破線0.427に重なる。

次の鎖を見てください。

P=(00.70.3100100)P = \begin{pmatrix} 0 & 0.7 & 0.3 \\ 1 & 0 & 0 \\ 1 & 0 & 0 \end{pmatrix}

状態Aからは必ずBかCへ行き、BとCからは必ずAへ戻ります。定常分布を計算すると π=(0.5, 0.35, 0.15)\pi = (0.5,\ 0.35,\ 0.15) で、一意に存在します。検算も通ります(πP=π\pi P = \pi)。

ところが状態Aから出発すると、分布はこうなります。

nn分布
0(1, 0, 0)(1,\ 0,\ 0)
1(0, 0.7, 0.3)(0,\ 0.7,\ 0.3)
2(1, 0, 0)(1,\ 0,\ 0)
3(0, 0.7, 0.3)(0,\ 0.7,\ 0.3)
延々と往復

PnP^n はどこにも収束しません。 定常分布はあるのに極限分布がない状態です。

区別を言葉にするとこうです。

定常分布は「πP=π\pi P = \pi という方程式の解」。極限分布は「PnP^n が実際に収束した先」。

前者は連立方程式が解ければ存在します。後者は収束しなければ存在しません。そして関係は片側だけです。

極限分布が存在するなら、それは必ず定常分布。逆は成り立たない。

区別せずに理解していても大半の場面で困らないのは、実務で出てくる鎖がほぼ非周期だからです。ただし試験では「極限分布が存在するか」を問われるので、区別が必要になります。

周期的でも「長期の滞在割合」は定常分布に一致する

救済策があります。極限分布が無くても、定常分布は意味を持ちます。

対数横軸の折れ線グラフ。青・赤・緑の3本の線が状態A・B・Cの滞在割合を表し、1ステップから2万ステップまでの累積割合を描いている。3本ともそれぞれの定常分布を示す同色の破線(0.5、0.35、0.15)に収束していく。

同じ周期的な鎖で、実際に2万ステップ歩かせて「各状態にいた時間の割合」を数えました。結果は (0.500, 0.350, 0.150)(0.500,\ 0.350,\ 0.150) で、定常分布に一致しました。

つまりこうです。

極限分布が無くても、時間平均は定常分布に収束する。

これが第8回の大数の法則とつながる部分です。マルコフ連鎖版の大数の法則をエルゴード定理と呼びます。主張は「1本の経路を長く追ったときの時間平均が、定常分布による平均に一致する」ことです。

ここは正確に言う必要があります。「多数の経路を用意して、固定した時刻 nn で平均する」ものとは一致しません。 それは π0Pn\pi_0 P^n そのもので、いま見たように周期的な鎖では収束しないからです。エルゴード定理が保証するのは時間方向の平均だけです。

数学的には、PnP^n 自体は振動しても、その平均 1nk=0n1Pk\frac{1}{n}\sum_{k=0}^{n-1} P^k(チェザロ平均)は収束します。実際に20万項まで足して確かめると、1行目は (0.5, 0.35, 0.15)(0.5,\ 0.35,\ 0.15) でした。

周期的な鎖でも定常分布は「長期的な滞在割合」として意味を持つと覚えておけば十分です。

定常分布の一意性が壊れるのはどんなときか

もう一つ、既約性のほうを確認します。一意性が壊れるのは、閉じたかたまりが2つ以上あるときだけです。

図 (b) の例で見ます。

P=(0.50.5000.50.500000.30.7000.70.3)P = \begin{pmatrix} 0.5 & 0.5 & 0 & 0 \\ 0.5 & 0.5 & 0 & 0 \\ 0 & 0 & 0.3 & 0.7 \\ 0 & 0 & 0.7 & 0.3 \end{pmatrix}

{A,B}\{A,B\}{C,D}\{C,D\} が完全に分離しています。すると定常分布の候補が無限に出てきます。

  • (0.5, 0.5, 0, 0)(0.5,\ 0.5,\ 0,\ 0) ← 検算OK
  • (0, 0, 0.5, 0.5)(0,\ 0,\ 0.5,\ 0.5) ← 検算OK
  • (0.3, 0.3, 0.2, 0.2)(0.3,\ 0.3,\ 0.2,\ 0.2) ← これも検算OK

混合したものも全部定常分布になってしまいます。固有値を見ると λ=1\lambda = 1 が2重になっていました。

一方、一時的な状態が混ざっているだけなら一意です。 その状態の確率がゼロになるだけで、解は1つに決まります。ドリル1の鎖(一時的な状態が3つある)も、定常分布は (0,0,0,1)(0,0,0,1) で一意でした。

まとめると、判定は「固有値1の重複度」=「閉じたかたまりの数」です。

鎖の構造定常分布極限分布
既約かつ非周期一意存在し、定常分布に一致
既約だが周期的一意存在しない(時間平均なら一致)
閉じたかたまりが2つ以上無限個初期分布によって行き先が変わる
一時的な状態を含む(閉じたかたまりは1つで非周期)一意(一時的な状態は0)存在する

収束するのは「状態」ではなく「分布」

もう1問、ドリルを出します。ここでも間違えました。

ドリル3:極限分布が存在するのはどれか

A=(0110)B=(0.50.510)C=(100.50.5)D=(0100010.500.5)A = \begin{pmatrix} 0 & 1 \\ 1 & 0 \end{pmatrix} \quad B = \begin{pmatrix} 0.5 & 0.5 \\ 1 & 0 \end{pmatrix} \quad C = \begin{pmatrix} 1 & 0 \\ 0.5 & 0.5 \end{pmatrix} \quad D = \begin{pmatrix} 0 & 1 & 0 \\ 0 & 0 & 1 \\ 0.5 & 0 & 0.5 \end{pmatrix}

PnP^n が収束するのはどれでしょうか。

筆者の答えは C だけ でした。理由は「収束するとはどこかの状態に落ち着くことだと思っていた」からです。

正解は A 以外の3つ(B, C, D)です。

ここに根本的な誤解がありました。

収束するのは「分布」であって、「状態」ではない。鎖そのものは永遠に動き回り続ける。

動き回りながら、「いまどこにいるかの確率」だけが一定値に落ち着く。これが定常分布です。

B で実測しました。40万ステップ歩かせると、状態が変わった回数は26.7万回。最後まで激しく動き続けています。にもかかわらず滞在割合は (0.6668, 0.3332)(0.6668,\ 0.3332) で、定常分布 (2/3, 1/3)(2/3,\ 1/3) に一致しました。

C を選んだのは、C だけが「本当にどこかに落ち着く」鎖だからです。C は1番目の状態が吸収状態(自分に確率1で留まる)なので、(1, 0)(1,\ 0) に到達したら文字通り動かなくなります。これは確かに収束しますが、収束の必要条件ではありません

A だけが収束しないのは、周期2だからです。(1,0)(1,0)(0,1)(0,1) を往復します。

D は面白い例です。3つの状態を順に巡って戻ってくる形ですが、3番目の状態には確率 0.50.5 の自己ループがあります(P33=0.5P_{33}=0.5)。この自己ループ1本のおかげで非周期になり、(0.25, 0.25, 0.5)(0.25,\ 0.25,\ 0.5) に収束します。自己ループが無ければ周期3で収束しませんでした。ドリル2のCと同じ理屈(互いに素な帰還歩数が2つできる)です。

この誤解を解いておくと、マルコフ連鎖の応用が理解しやすくなります。PageRankが「順位に落ち着く」のも、サーファーが1ページに留まるのではなく、動き回った結果の滞在割合が安定するという話です。

吸収確率:立式の型と、マイナスが出る理由

ここからは「いつか終わる鎖」の話です。試験でよく問われます。

ドリル4:購入に至る確率を求めてください

無料ユーザーの状態遷移です。

現在\翌週無料検討中購入離脱
無料0.50.30.150.05
検討中0.10.40.40.1
購入0010
離脱0001

「無料」から出発したとき、離脱ではなく購入に至る確率はいくらでしょうか。

筆者はこれを解こうとして、答えがマイナスになりました。原因が2つあったので、両方書きます。

型:1歩先で場合分けして、自分に戻る項を右辺に残す

まず正しい立式です。3ステップの型として覚えます。

① 1歩先で場合分けする。 求めたい確率を hih_i として、

hi=jPijhjh_i = \sum_j P_{ij}\, h_j

「次にどこへ行くかで分解する」だけです。

② 吸収状態の値を決め打ちする。 h購入=1h_{\text{購入}} = 1h離脱=0h_{\text{離脱}} = 0

③ 内部状態について解く。

この問題ならこうなります。

hA=0.5hA+0.3hB+0.15×1+0.05×0h_A = 0.5\,h_A + 0.3\,h_B + 0.15 \times 1 + 0.05 \times 0 hB=0.1hA+0.4hB+0.40×1+0.10×0h_B = 0.1\,h_A + 0.4\,h_B + 0.40 \times 1 + 0.10 \times 0

最大の注意点は、自分自身に戻る 0.5hA0.5\,h_A を書き忘れないことです。1歩進んだ先が自分自身でも、そこからの見通しは同じ hAh_A なので、この項が必要です。

移項すると、

0.5hA0.3hB=0.15,0.1hA+0.6hB=0.400.5\,h_A - 0.3\,h_B = 0.15, \qquad -0.1\,h_A + 0.6\,h_B = 0.40

hBh_B 側の符号がマイナスになる点に注意してください。解くと hA=7/9=0.7778h_A = 7/9 = 0.7778hB=0.7963h_B = 0.7963 です。200万回のシミュレーションで 0.7779±0.00030.7779 \pm 0.0003 でした。

原因1:移項の符号

もし移項で符号を ++ にしてしまうと、

(0.50.30.10.6)(hAhB)=(0.150.40)\begin{pmatrix} 0.5 & 0.3 \\ 0.1 & 0.6 \end{pmatrix}\begin{pmatrix} h_A \\ h_B\end{pmatrix} = \begin{pmatrix} 0.15 \\ 0.40 \end{pmatrix}

を解くことになり、答えは hA=0.111h_A = -0.111 です。マイナスが再現できました。

検算の習慣を2つ書いておきます。吸収確率は必ず0以上1以下なので、マイナスや1超えが出たら即座に立式ミスと判断できます。もう1つは、購入確率と離脱確率の両方を計算して足して1になるかを見ることです(0.7778+0.2222=10.7778 + 0.2222 = 1)。

原因2:吸収状態について立式してしまった

もう1つの原因がこれで、こちらのほうが本質的です。購入と離脱についても式を立てようとすると、行き詰まります。

4状態すべてについて hi=jPijhjh_i = \sum_j P_{ij}h_j を書き出してみます。

h無料=0.5h無料+0.3h検討+0.15h購入+0.05h離脱h_{\text{無料}} = 0.5\,h_{\text{無料}} + 0.3\,h_{\text{検討}} + 0.15\,h_{\text{購入}} + 0.05\,h_{\text{離脱}} h検討=0.1h無料+0.4h検討+0.4h購入+0.1h離脱h_{\text{検討}} = 0.1\,h_{\text{無料}} + 0.4\,h_{\text{検討}} + 0.4\,h_{\text{購入}} + 0.1\,h_{\text{離脱}} h購入=1h購入h_{\text{購入}} = 1 \cdot h_{\text{購入}} h離脱=1h離脱h_{\text{離脱}} = 1 \cdot h_{\text{離脱}}

下の2本は 0=00 = 0 です。 移項すると何も残りません。方程式の形をしているのに情報量がゼロです。

確認すると、IPI - Pランクは2でした。未知数が4個あるのに使える式が2本しかない。だから解が一意に決まりません

ここで、IPI-P の行列式が 00 であること自体は証拠になりません。行の和が1である行列なら、行列式は必ず 00 です(IP)1=0(I-P)\mathbf{1} = \mathbf{0} なので、全成分1のベクトルが必ず核に入る)。吸収状態のない既約な鎖でも det(IP)=0\det(I-P)=0 で、ランクは n1n-1 です。

本質的なのはどこまでランクが落ちたかです。42=24 - 2 = 2 次元の自由度が残っていて、この数は「閉じたかたまりの数」(この鎖では吸収状態が2つなのでちょうど2)です。だから境界条件が2本必要という関係になります。

つまり、この2つの値は方程式から求めることができず、外から与えるしかないのです。

境界条件 h購入=1h_{\text{購入}}=1h離脱=0h_{\text{離脱}}=0 は天から降ってくる約束事ではなく、方程式が縮退した穴を埋めるために必要な情報。

これが「境界条件」の正体でした。教科書では最初から h=1,0h=1, 0 と与えられているので、なぜそれが必要なのかが見えにくい部分です。

1点だけ精密にしておきます。自己ループそのものは原因ではありません。 「無料」も 0.50.5 の自己ループを持っていますが、何の問題も起こしません。

それどころか、自己ループは吸収確率をまったく変えません。 「無料」の 0.50.5 を残りの行き先に比例配分して消した鎖 (0, 0.6, 0.3, 0.1)(0,\ 0.6,\ 0.3,\ 0.1) で計算すると、hA=7/9h_A = 7/9元とぴったり同じでした。「同じ場所に留まる」という足踏みは、最終的にどちらへ落ちるかに影響しないので当然です。

ただし時間には影響します。吸収までの平均ステップ数は 3.3333.333 から 2.2222.222 へ減りました。足踏みしている分だけ長引く、というのは直感どおりです。行き先は変わらないが、かかる時間は変わる。

問題なのは自己ループの確率が 1.01.0 であること、つまり出られないことです。0.50.5 なら hA0.5hA=0.5hAh_A - 0.5h_A = 0.5h_A が残るので式として機能します。1.01.0 だと hh=0h - h = 0 で完全に消えてしまう。係数が1のときだけ左辺と右辺が打ち消し合う、というのがからくりです。

だから実務的な作法は、吸収状態を最初から未知数から外すことです。基本行列 N=(IQ)1N = (I-Q)^{-1}QQ が「一時的な状態だけを取り出した小さい行列」になっているのは、この理由です。4×44\times4 ではなく 2×22\times2 を扱います。

そして基本行列を使うと、吸収確率が一発で出ます。

B=NR=(IQ)1RB = N R = (I-Q)^{-1} R

RR は「一時的な状態から吸収状態への確率」の部分行列です。計算すると

B=(0.77780.22220.79630.2037)B = \begin{pmatrix} 0.7778 & 0.2222 \\ 0.7963 & 0.2037 \end{pmatrix}

1列目が購入確率、2列目が離脱確率で、各行の和が1になっています。連立方程式を解いた結果と一致します。

定常分布と吸収確率で立式の向きが逆になる

型が似ているので、違いを明示します。どちらも連立方程式ですが、行と列の使い方が逆です。

定常分布吸収確率
πP=π\pi P = \pihi=jPijhjh_i = \sum_j P_{ij} h_j
未知数の並び行ベクトル π\piPP の左から掛ける)列ベクトル hhPP の右から掛ける)
行列の読み方列を縦に読む行を横に読む
意味その状態に流れ込む量の合計1歩先での場合分け
足りない条件の補い方1本捨てて πi=1\sum \pi_i = 1 を足す吸収状態の値を決め打ち

定常分布では「1本捨てて和が1を足す」のに対し、吸収確率では「境界条件が最初から与えられている」ので、後者のほうが手続きは楽です。

賭博者の破産:勝率1%の差が結果を分ける

この型の代表例が賭博者の破産です。試験にも出ます。

左右2枚のグラフ。左は手持ちの金額(横軸0から100)に対する100円到達確率で、p=0.45の赤・p=0.49の橙・p=0.5の青の直線・p=0.51の緑の4本の曲線が描かれ、a=50, p=0.49のとき0.119という注釈が矢印で示されている。右は勝負が決まるまでの平均回数を示す紫の放物線で、a=50で最大値2500回に達する山型。

手持ち aa 円から始めて、1回の勝負で確率 pp で1円増え、確率 q=1pq=1-p で1円減る。NN 円に到達すれば勝ち、0円になれば負け。勝つ確率はいくらか。

型に当てはめます。

ha=pha+1+qha1,h0=0, hN=1h_a = p\,h_{a+1} + q\,h_{a-1}, \qquad h_0 = 0,\ h_N = 1

これを解くと閉じた式になります。

ha={aN(p=0.5)1ra1rN,r=qp(p0.5)h_a = \begin{cases} \dfrac{a}{N} & (p = 0.5) \\[8pt] \dfrac{1-r^a}{1-r^N},\quad r = \dfrac{q}{p} & (p \ne 0.5) \end{cases}

a=50a=50N=100N=100 で計算すると、勝率がわずか1%違うだけで結果が激変します。

勝率 pp100円に到達する確率
0.450.0000(ほぼゼロ)
0.490.1192
0.500.5000
0.510.8808

p=0.49p=0.49p=0.51p=0.510.1190.1190.8810.881 です。1回あたり2%の差が、最終的に7倍以上の差になる。連立方程式を線形代数で解いた結果と閉じた式が一致することも確認しました(どちらも 0.1191750.119175)。

これが「続けるほどカジノ側が有利になる」ことの数学的な内容です。ハウスエッジがわずかでも、繰り返せば確実に効きます。

なお NN が有限なら、プレイヤーが勝つ確率がゼロになるわけではありません(p=0.45p=0.45 でも 4.4×1054.4 \times 10^{-5} あります)。プレイヤーが必ず破産すると言い切れるのは、相手が無限の資金を持っている場合(NN \to \infty)です。このとき p<0.5p < 0.5 なら破産確率がちょうど1になります。

平均到達時間は「+1+1」を足すだけ

平均到達時間も同じ型で解けます。右辺に「1歩使った」を表す +1+1 を足すだけです。

ki=1+jPijkjk_i = 1 + \sum_j P_{ij}\, k_j

賭博者の破産で p=0.5p=0.5 のとき、答えは驚くほどきれいです。

ka=a(Na)k_a = a(N-a)

a=50a=50N=100N=100 なら 50×50=250050 \times 50 = 2500 回。連立方程式を解いても 2500.00002500.0000 でした。図の右側の山型がこれで、真ん中から始めると最も長引くことがわかります。

a=3a=3N=10N=10 なら 3×7=213 \times 7 = 21 回。手計算で確かめられる大きさなので、試験対策として覚えておくと便利です。

平均回帰時間は 1/π1/\pi:逆行列は要らない

吸収確率で逆行列が出てきたので、混同しやすい点を整理します。

ドリル5:平均回帰時間

P=(0.50.50001100)P = \begin{pmatrix} 0.5 & 0.5 & 0 \\ 0 & 0 & 1 \\ 1 & 0 & 0 \end{pmatrix}

状態0を出てから、次に状態0に戻るまでの平均ステップ数はいくらでしょうか。

答えは2です。 そしてこれは逆行列を使わずに解けます

定常分布を求めると π=(0.5, 0.25, 0.25)\pi = (0.5,\ 0.25,\ 0.25)。そして

平均回帰時間=1πi\text{平均回帰時間} = \frac{1}{\pi_i}

なので 1/0.5=21/0.5 = 2。40万ステップのシミュレーションで 2.00012.0001 でした。

なぜ逆数になるのか。 長期的に見て、全ステップのうち πi\pi_i の割合を状態 ii で過ごすということは、1/πi1/\pi_i ステップに1回のペースで状態 ii を訪れるということです。πi=0.5\pi_i = 0.5 なら2ステップに1回訪れる。だから平均間隔は2。

この公式が使えるのはその状態が正再帰的なときです(有限で既約な鎖ならすべての状態が自動的に正再帰的になります)。一時的な状態や零再帰的な状態では πi=0\pi_i = 0 なので 1/πi1/\pi_i が無限大になり、これは「戻ってくる平均時間が無限」という正しい答えではあるものの、公式として使う場面ではありません。

使い分けを明示します。基本行列(逆行列)が必要なのは、吸収がある鎖だけです。

問題必要な計算
平均回帰時間(自分に戻る)定常分布を求めて逆数を取る。逆行列は不要
吸収までの平均時間基本行列 N=(IQ)1N=(I-Q)^{-1} の行和
ある特定の状態への平均到達時間連立方程式 ki=1+jPijkjk_i = 1 + \sum_j P_{ij}k_j(目標を吸収扱いにする)

待ち行列と出生死滅過程

ここからは応用です。試験の比重は小さいので、用語と使いどころを中心にします。

出生死滅過程は「隣にしか動けない鎖」

出生死滅過程(birth-death process)は、状態が整数で、1つ上か1つ下にしか動けない鎖です。人口、在庫、待っている人数などが対象になります。

この構造だと計算が非常に楽になります。詳細つりあいという関係が使えるからです。

πnλn=πn+1μn+1\pi_n \lambda_n = \pi_{n+1} \mu_{n+1}

nn から n+1n+1 へ流れる量」と「n+1n+1 から nn へ流れる量」が釣り合う、という式です。隣同士しか行き来しないので、境界を通る流量が両方向で等しくなります。

これを使うと、πn+1=πnλn/μn+1\pi_{n+1} = \pi_n \cdot \lambda_n/\mu_{n+1}比を掛けていくだけで定常分布が出ます。

席が5席までの店(λ=0.8\lambda = 0.8μ=1.0\mu = 1.0)で計算すると、

π=(0.2711, 0.2168, 0.1735, 0.1388, 0.1110, 0.0888)\pi = (0.2711,\ 0.2168,\ 0.1735,\ 0.1388,\ 0.1110,\ 0.0888)

満席で断られる確率が 0.08880.0888 です。生成行列から πQ=0\pi Q = 0 を解いた結果と 2.2×10162.2 \times 10^{-16} まで一致しました。

M/M/1:稼働率を上げると待ちが爆発する

記号の読み方から。M/M/1 は「到着がマルコフ的(ポアソン到着)/サービス時間がマルコフ的(指数分布)/窓口が1つ」の略で、M は Markovian です。第5回の無記憶性が2箇所に効いています。

知りたいのは、到着率 λ\lambda と処理率 μ\mu から、平均何人待っているか平均何分待つかです。

左右2枚のグラフ。左は稼働率ρ(横軸0から0.97)に対する系内の平均人数Lを示す赤い曲線で、ρ=0.5で1人、0.8で4人、0.9で9人、0.95で19人の点に注釈が付き、右端で急激に立ち上がっている。右はρ=0.8のときの系内人数の分布を示す棒グラフで、青が理論値の幾何分布、橙がシミュレーション値で、n=0から12までほぼ一致している。

公式は3つだけです。

ρ=λμ,W=1μλ,Wq=ρμλ\rho = \frac{\lambda}{\mu}, \qquad W = \frac{1}{\mu - \lambda}, \qquad W_q = \frac{\rho}{\mu-\lambda}

WW が平均滞在時間(待ち時間+サービス時間)、WqW_q が平均待ち時間です。人数は

L=ρ1ρ,Lq=ρ21ρL = \frac{\rho}{1-\rho}, \qquad L_q = \frac{\rho^2}{1-\rho}

分母が μλ\mu - \lambda(さばく速さ引く来る速さ)という形が全部の中心です。余力が小さいほど待ちが長い、というだけの式です。

前提として ρ<1\rho < 1λ<μ\lambda < \mu)が必要です。来る速さがさばく速さを超えると分母が0以下になり、行列は無限に伸び続けて定常状態が存在しません。

そして L=λWL = \lambda W という関係をリトルの法則と呼びます。実測で λWsim=3.925\lambda W_{\text{sim}} = 3.925Lsim=3.909L_{\text{sim}} = 3.909 と一致しました。この法則はM/M/1に限らず非常に広い条件で成り立つので、覚えておくと便利です。

図の左が伝えたいことは1点です。

稼働率を上げると、待ちは線形ではなく爆発的に増える。

稼働率 ρ\rho系内の平均人数 LL実測(40万人)
0.51.01.007
0.84.03.909
0.99.08.687
0.9519.018.318

稼働率を1.9倍(0.5→0.95)にしただけで、系内の人数が19倍になります。待っている人数 LqL_q で見ればもっと激しく、0.50.5 人から 18.0518.05 人で36倍です。「CPU使用率が50%だともったいないから90%まで詰めよう」という判断がレイテンシを壊すのは、この式のせいです。

図の右は、定常分布が幾何分布 P(n)=(1ρ)ρnP(n) = (1-\rho)\rho^n になることの確認です。理論値とシミュレーションが小数第2位まで一致しました。ここでも第5回の幾何分布が再登場します。

ドリル6:問い合わせ対応の待ち時間

問い合わせが1時間あたり平均12件届き、担当者1人が1時間あたり平均15件さばけます。

(a) 稼働率 ρ\rho はいくらか。 (b) 平均待ち時間 WqW_q は何分か。 (c) 処理能力を15件から16件に上げると WqW_q は何分になるか。

(a) ρ=12/15=0.8\rho = 12/15 = 0.8

(b) Wq=0.8/(1512)=0.2667W_q = 0.8/(15-12) = 0.2667 時間 = 16分

(c) ρ=12/16=0.75\rho = 12/16 = 0.75Wq=0.75/(1612)=0.1875W_q = 0.75/(16-12) = 0.1875 時間 = 11.25分

(c) が要点です。処理能力を6.7%上げただけで、待ち時間が29.7%減りました。 分母の余力が 343 \to 4 と1.33倍になるので、非線形に効きます。

これは実務で使える見方です。ρ\rho が1に近いシステムでは、わずかな増強が劇的に効く一方、わずかな負荷増が致命傷になる。 「あと1割の負荷なら耐えられるだろう」という直感が外れる領域です。

なお、この滞在時間は指数分布 Exp(μλ)\mathrm{Exp}(\mu-\lambda) に従います。50万人のシミュレーションで確認すると、平均 19.6219.62 分に対して標準偏差 19.3219.32 分でした。指数分布は平均と標準偏差が等しいのでこれが証拠です。分位点も合いました。

分位点実測指数分布の理論値
50%13.66分13.86分
90%45.44分46.05分
95%58.83分59.91分
99%88.32分92.10分

平均20分でも95%の人は59分以内、99%の人は92分という裾の長さが読めます。平均だけ見ていると危険だという、第5回第3回と同じ話です。

連続時間マルコフ連鎖

ここまでは「1歩、2歩」と離散的に進む鎖でした。M/M/1のように時間が連続的に流れる場合は、推移確率行列 PP の代わりに生成行列 QQ(推移速度行列)を使います。

対応関係だけ書きます。

離散時間連続時間
基本の量推移確率 PijP_{ij}推移速度 qijq_{ij}
行の制約各行の和が 1各行の和が 0(対角が負)
定常分布の式πP=π\pi P = \piπQ=0\pi Q = \mathbf{0}
状態に留まる時間幾何分布指数分布

離散なら幾何分布、連続なら指数分布という対応が、そのまま第5回第6回の関係になっています。無記憶性を持つ分布が、離散版と連続版で1つずつある、というだけの話です。

定常分布の式が πP=π\pi P = \pi から πQ=0\pi Q = \mathbf{0} に変わりますが、これは同じ式の書き換えです。微小時間 dtdt の推移行列が P(dt)=I+Qdt+o(dt)P(dt) = I + Q\,dt + o(dt) と書けるので、πP(dt)=π\pi P(dt) = \pi に代入すると πQdt=0\pi Q\,dt = \mathbf{0}、つまり πQ=0\pi Q = \mathbf{0} になります。生成行列 QQ の正体は「単位時間あたりの推移確率の変化率」で、Q=limdt0(P(dt)I)/dtQ = \lim_{dt \to 0}(P(dt)-I)/dt です。行和が0になるのは、P(dt)P(dt) の行和が1で II の行和も1だからです。

Web・ITでの応用:PageRankは定常分布そのもの

ここからは実務での使いどころです。

左は4ページA・B・C・Dのリンク構造を示す状態遷移図。A→B(0.5)、A→C(0.5)、B→C(1)、C→A(1)、D→C(1)の矢印が描かれている。右はPageRankの棒グラフで、A=0.373、B=0.196、C=0.394、D=0.037。Cが最上位でDが最下位。

PageRankは、リンクをランダムに辿り続けたときの定常分布です。 これは比喩ではなく、文字通りそうです。

4ページの小さな例で計算します。A から B と C へリンク、B から C へ、C から A へ、D から C へ。リンク行列を HH とすると、HijH_{ij} は「ページ ii からリンクをランダムに1本選んだとき jj へ行く確率」です。

ただし HH をそのまま使うと壊れます。行き止まりのページがあると行の和が1にならず、既約でないグラフだと一意性が崩れる。そこで実際のPageRankは

G=dH+(1d)11n,d=0.85G = d\,H + (1-d)\frac{\mathbf{1}\mathbf{1}^\top}{n}, \qquad d = 0.85

と混ぜます。これは「15%の確率でランダムなページに飛ぶ」という操作で、全成分を正にして既約かつ非周期にするための細工です。

つまりダンピング係数 dd の正体は、極限分布の存在と一意性を強引に保証する装置でした。この章で学んだ条件がそのまま出てきます。

計算結果は (0.3725, 0.1958, 0.3941, 0.0375)(0.3725,\ 0.1958,\ 0.3941,\ 0.0375)。被リンク3本の C が最上位、被リンク0本の D が最下位です。D の 0.03750.03750.15/40.15/4 にちょうど一致していて、これが「ランダムジャンプで流れ込む分だけ」という下限になっています。

実装はべき乗法です。初期分布に GG を掛け続けるだけ。

反復回数状態Aの確率
10.2500
20.5209
50.3885
100.3751
200.3725
400.3725

40回で固有ベクトルを直接解いた結果と 1.5×10141.5 \times 10^{-14} まで一致しました。PnP^n を掛け続けると定常分布に収束する」という定理が、そのままアルゴリズムになっています。

収束の速さは第2固有値の絶対値で決まります。d=0.85d=0.85 なら第2固有値の絶対値が 0.850.85 以下に抑えられるので、数十回で収束します。実際のWebは数十億ページありますが、この性質のおかげで計算可能になっています。

他のWeb・IT系の応用も挙げておきます。

応用状態何を知りたいか
PageRankページ定常分布=重要度
サイト内回遊の分析ページ定常分布、離脱までの平均ページ数(吸収確率)
コンバージョン分析ユーザーの段階吸収確率(購入 vs 離脱)
解約予測契約状態吸収までの平均時間=顧客生涯
サーバー容量設計待っているリクエスト数M/M/1、ρ\rho と待ち時間
A/Bテストの逐次判定累積の差ランダムウォークの吸収確率
MCMC(次回以降)パラメータ定常分布を「作る」ことで事後分布から標本を得る

最後の MCMC(Markov Chain Monte Carlo=マルコフ連鎖モンテカルロ法)が実は最も重要な応用です。この章では「与えられた PP から定常分布 π\pi を求める」ことをしてきましたが、MCMCは逆向きに、欲しい π\pi を定常分布に持つ PP を設計するという発想です。ベイズ統計の計算がこれで回ります。

HMMとの関係:状態が見えるかどうか

筆者は学生時代の研究でHMM(Hidden Markov Model=隠れマルコフモデル)を使っていました。ただしライブラリに投げて結果を使う使い方で、内部計算は自分でやっていませんでした。

この章との関係を整理しておきます。

上下2段の図。上段には赤い円Z1・Z2・Z3・Z4が横に並び、赤い矢印で左から右へ繋がっている(隠れ状態のマルコフ連鎖・見えない)。下段には青い四角X1・X2・X3・X4が並び、それぞれ上の円から青い縦矢印が下りている(観測・見える)。下段の四角同士を結ぶ横矢印はない。

マルコフ連鎖は「状態が見える」場合、HMMは「状態が見えない」場合という理解で正しいです。図の上段だけがこの章の対象で、HMMは上段が隠れて下段だけが観測できる構造です。

ここで1つ、押さえておくと綺麗な事実があります。

観測列そのものはマルコフ連鎖になりません。 隠れ状態がマルコフ的でも、観測を通すと記憶が漏れます。

数値で確認しました。隠れ状態の推移が A=(0.950.050.100.90)A = \begin{pmatrix} 0.95 & 0.05 \\ 0.10 & 0.90\end{pmatrix}、出力確率が B=(0.90.10.30.7)B = \begin{pmatrix} 0.9 & 0.1 \\ 0.3 & 0.7 \end{pmatrix} のHMMから200万ステップ生成して、観測列の依存を測ります。

条件次に0が出る確率
直前が00.7968
直前が0、その前も00.8369
直前が0、その前は10.6395

直前が同じ0でも、2つ前が違うと確率が 0.8370.8370.6400.640 で大きく変わります。 マルコフ性が成り立つならこの2つは一致するはずなので、観測列はマルコフ連鎖ではありません。

理由は考えれば当然です。観測が0でも、隠れ状態が1なのか2なのか確定しません。過去の観測は「いまどちらの隠れ状態にいそうか」という情報を追加で持っているので、効いてしまいます。

逆に言えば、「観測がマルコフ的でないデータを、マルコフ的な隠れ状態で説明する」のがHMMの狙いです。この章の道具は、その土台にあたります。

なお、隠れ状態の定常分布は (2/3, 1/3)(2/3,\ 1/3) で、観測の定常的な出現確率は πAB=(0.7, 0.3)\pi_A B = (0.7,\ 0.3) と計算できます。隠れ状態の定常分布に出力確率を掛けるだけです。

試験範囲について。 準1級のワークブックでは、HMMは正面から出題されません。前向き・後ろ向きアルゴリズムやバウム・ウェルチ法(EMアルゴリズムの一種)は範囲外です。この章は「状態が見えるマルコフ連鎖」までなので、HMMは深追いしなくて大丈夫です。

試験で問われる計算パターン

まとめとして、典型的な出題の型を整理します。

パターン1:nn 歩後の分布を求める。 PPnn 乗して初期分布に掛けるだけです。C-K方程式の名前で聞かれても同じことです。手計算なら n=2,3n=2, 3 程度なので素直に掛けます。

パターン2:定常分布を求める。 πP=π\pi P = \pi を成分で書き、列を縦に読んで式を立てる。1本捨てて πi=1\sum \pi_i = 1 を足す。最後に πP\pi P に代入して検算する。2状態なら「π\pi の比 = 相手から来る確率の比」で暗算できます。

パターン3:平均回帰時間。 1/πi1/\pi_i。定常分布が出ていれば割り算1回です。

パターン4:吸収確率。 hi=jPijhjh_i = \sum_j P_{ij}h_j に境界条件。自分自身に戻る項を右辺に残すのを忘れない。答えが0〜1に入るか検算する。

パターン5:吸収までの平均時間。 上の式の右辺に +1+1 を足す。または基本行列 N=(IQ)1N=(I-Q)^{-1} の行和。

パターン6:状態の分類。 既約か(一周できるか)、非周期か(既約なら対角に正の成分があるかで判定できる)、各状態が再帰的か一時的か(抜け道があるか)。

パターン7:M/M/1の待ち時間。 ρ=λ/μ\rho = \lambda/\mu を出してから W=1/(μλ)W = 1/(\mu-\lambda)Wq=ρ/(μλ)W_q = \rho/(\mu-\lambda)。単位(時間か分か)に注意。

自分が間違えていたこと

ドリルで4か所、加えて解説の準備中に2か所、自分の誤りが見つかりました。全部書いておきます。

① 再帰的と一時的を「動くか」で判定していた。 「そこから動かないので一時的」と考えましたが逆でした。判定基準は帰還確率だけです。recurrent は「再発する」、transient は「一過性の」という語感で覚えると間違えません。自分が正しく予想した定常分布 (0,0,0,1)(0,0,0,1) の中に、この問いの答えが入っていたというのが決定的でした。0が並ぶ位置が一時的、1になる位置が再帰的です。

② 既約を「バラバラでないこと」と理解していた。 定義は往復できることです。吸収状態があると既約ではなくなります。

③ 周期を平均で答えた。 3歩と4歩で「3.5」としましたが、周期は最大公約数なので1(非周期)です。周期は整数にしかなりません。さらに16と24なら最小公倍数の48ではなく最大公約数の8で、「8歩ごとに戻れる」のではなく「戻れるのは8の倍数だけ」という意味でした(8自身には戻れません)。

④ 収束するのは「状態」だと思っていた。 収束するのは分布です。鎖は永遠に動き続けます。40万ステップ中26.7万回も状態が変わっているのに滞在割合は一定、という実測がこれを示しています。

⑤ ランダムウォークの帰還確率を間違えた。 p=0.6p=0.6 のとき min(1,q/p)=0.667\min(1, q/p) = 0.667 と書きましたが、これは+1+1 の位置から 00 へ下りる確率」であって原点からの帰還確率ではありません。正しくは1歩目で場合分けして

f00=pmin ⁣(1,qp)+qmin ⁣(1,pq)=0.6×0.667+0.4×1=0.8f_{00} = p \cdot \min\!\left(1, \frac{q}{p}\right) + q \cdot \min\!\left(1, \frac{p}{q}\right) = 0.6 \times 0.667 + 0.4 \times 1 = 0.8

シミュレーションの 0.80330.8033 と一致します。公式を丸暗記すると「どこからどこへの確率か」を取り違えるという典型例でした。

⑥ 定常分布の一意性が壊れる条件を誤解した。 「誰からもリンクされないページがあると一意性が壊れる」と考えましたが誤りです。一意性が壊れるのは閉じたかたまりが2つ以上あるときだけで、一時的な状態があるだけなら一意です(その状態の確率が0になる)。固有値1の重複度で確認できます。閉じたかたまりが2つある鎖では重複度2、一時的な状態があるだけの鎖では重複度1でした。

つまり既約性は一意性の十分条件であって必要条件ではない、ということです。記事の前半で「既約性が崩れると一意性が壊れる」と書きたくなったのですが、正しくは「閉じたかたまりが2つ以上になると壊れる」でした。

この記事の要点

  • 状態の分類の4用語は判定対象のレベルが違う。 既約と周期は鎖全体、再帰的/一時的は個々の状態。壊れるものも別で、閉じたかたまりが2つ以上になると定常分布の一意性が壊れ、周期性が残ると極限分布の存在が壊れる
  • 既約性は一意性の十分条件であって必要条件ではない。 既約でなくても閉じたかたまりが1つなら定常分布は一意(一時的な状態の確率が0になるだけ)
  • 再帰的/一時的は「動くか」ではなく「戻るか」。 recurrent=再発する、transient=一過性の。吸収状態こそが再帰的で、活発に行き来する状態が一時的になる逆転が起きる
  • 定常分布が一意で (0,0,0,1)(0,0,0,1) なら、0が並ぶ位置が一時的で1の位置が再帰的。 分類と定常分布は同じことを言っている
  • 一時的な状態は訪問回数の期待値が有限。 基本行列 N=(IQ)1N=(I-Q)^{-1} で計算でき、行和が吸収までの平均時間(実測 16.698±0.03016.698 \pm 0.030 対 理論 16.66716.667
  • 周期は最大公約数。最小公倍数ではない。 16と24なら8。しかも「戻れるのは8の倍数だけ」で8自身には戻れない
  • 既約な鎖なら、対角成分に正の値が1つでもあれば非周期。 自己ループ1本で振動が消える(既約でないと成り立たない。反例あり)
  • 定常分布は方程式の解、極限分布は PnP^n の収束先。 極限分布があれば定常分布だが、逆は成り立たない。周期的な鎖では π=(0.5,0.35,0.15)\pi=(0.5,0.35,0.15) が一意に存在するのに PnP^n は往復し続ける
  • 極限分布が無くても時間平均は定常分布に一致する(エルゴード定理)。実測 (0.500,0.350,0.150)(0.500, 0.350, 0.150)。ただし保証されるのは時間方向の平均だけで、固定時刻での多数経路の平均は収束しない
  • 収束するのは分布であって状態ではない。 40万ステップ中26.7万回状態が変わっても滞在割合は (0.6668,0.3332)(0.6668, 0.3332) で一定
  • C-K方程式は有限状態なら行列積そのもの(差 1.1×10161.1\times10^{-16})。名前が付いているのは連続状態では積分になるから
  • 定常分布は列を縦に読み、吸収確率は行を横に読む。 向きが逆
  • 吸収状態について立式すると 0=00=0 になる。 IPI-P が4状態でランク2まで落ちるので、境界条件 h=1,0h=1,0 は縮退した穴を埋める必需品(det(IP)=0\det(I-P)=0 自体はどんな確率行列でも成り立つので証拠にならない。残った自由度の数=閉じたかたまりの数を見る)
  • 自己ループは吸収確率を変えない(比例配分で消しても 7/97/9 のまま)。係数が 1.01.0 のときだけ左右が打ち消し合って式が消える
  • 吸収確率の検算は「0〜1に入るか」「2つ足して1になるか」。 移項の符号ミスで 0.111-0.111 が出た
  • 賭博者の破産は勝率1%差で結果が激変。 a=50,N=100a=50, N=100p=0.49p=0.49 なら 0.1190.119p=0.51p=0.51 なら 0.8810.881
  • 平均到達時間は同じ式に +1+1 を足すだけ。 p=0.5p=0.5 なら a(Na)a(N-a)50×50=250050\times50=2500
  • 平均回帰時間は 1/πi1/\pi_i で逆行列は不要(既約な鎖に限る)。逆行列が要るのは吸収がある鎖だけ
  • M/M/1は稼働率を上げると待ちが爆発する。 ρ\rho を1.9倍(0.5→0.95)にすると系内人数 LL は19倍、待ち人数 LqL_q は36倍。逆に処理能力を6.7%上げると待ちが29.7%減る
  • M/M/1の滞在時間は指数分布(平均19.62分・標準偏差19.32分)。平均20分でも99%点は92分
  • 出生死滅過程は隣にしか動けない鎖。 詳細つりあいで比を掛けるだけで定常分布が出る
  • PageRankは定常分布そのもの。 ダンピング係数 0.850.85 の正体は「既約かつ非周期にして極限分布の存在と一意性を保証する装置」
  • 観測列はマルコフ連鎖にならない。 HMMで実測すると2つ前の観測が効く(0.8370.8370.6400.640)。準1級ではHMMは範囲外

おわりに

前提知識があるので復習で済むと思っていたのに、6問中4問で用語の定義を間違えました。ただ、間違え方に一貫性があったので学びは大きかったです。

概念は理解できていて、用語のラベルだけが繋がっていなかった、という状態でした。定常分布が (0,0,0,1)(0,0,0,1) になると正しく予想できたのに、それを「一時的」という言葉に翻訳できなかったのが典型です。

この連載で何度も同じことが起きています。第16回では決定係数の用語が分からずに詰まりました。用語を「知っている概念の名前」として受け取り直す作業が、独学ではいちばん抜けやすいのだと思います。

そして今回、実務との接続が最も濃い章でもありました。PageRankがこの章の内容そのままだと確認できたのは収穫です。ブログを運営しているので、検索順位の仕組みの土台が定常分布だというのは面白い発見でした。M/M/1の「稼働率を上げると待ちが爆発する」も、サーバー設計で直接使える形です。

次回は第15章・確率過程の基礎です。今回の「離散的に1歩ずつ進む鎖」が、時間だけを連続にしたポアソン過程と、時間も状態も連続にしたブラウン運動に広がります。

主役は2つあります。1つはポアソン過程で、今回のM/M/1で「客が到着する」と言っていた部分の中身です。λ\lambda が単位時間あたりの率として意味を持つ理由がわかります。もう1つがブラウン運動で、「分散が時間に比例する」という新しい性質が加わります。第8回n\sqrt{n} の壁が、σt\sigma\sqrt{t} という形で顔を出す回です。

今回の再帰的・一時的という用語も、ランダムウォークで再登場します。1次元と2次元は再帰的なのに3次元以上は一時的という、次元をまたぐ断崖をシミュレーションで確かめます。