吸収マルコフ連鎖はAbsorbing Markov chain (Wikipedia)に説明がありますが、かなり分かりづらいので「一次元のランダムウォーク」の例を使って説明していきます
【例題 1】（ランダムウォーク）酔っ払いが家に帰るかバーに行くか？
ある男性が街の0~4の地点を歩いています（下図参照）。もし彼が1、2、3の時には1/2の確率で左か右へ歩きます。彼は4（バー）か0（家）に到達するまで歩き続けます。自宅かバーのどちらかに到達した場合もう動きません。
地点1にいた男性がゴール（バーまたは家）に到達するまで歩く回数の期待値を求めよ。
このようにもう動かないゴールがあるものを吸収マルコフ連鎖と呼び、その定義は以下のようになります
- 少なくとも1つの吸収状態がある
- 任意の状態から少なくとも1つの吸収状態へ有限数のステップで移行することが可能
なんとかイメージは出来たと思いますので。これを確率行列 (Wikipedia)を使って見ていきます。
例題 1の確率行列Pを求める
(今回の計算では行が移動元、列が移動先となる右確率行列を使います)
numpyを使って確率行列Pを求めます。
このPを見ると
- 0,4行を見ると自分への確率が１なので吸収状態である
- 0,4列を見ると他の状態からの遷移がある
このように確率行列を見れば吸収マルコフ連鎖であることの判別が出来ます。
# Transition matrix of Drunkard’s Walk
import numpy as np
n = 5 # number of state
P = np.zeros((n+1, n+1))
for i in range(n+1):
if i in [0,n]: # absorbing states
P[i, i] = 1
continue
for d, p in [(-1,1/2),(1,1/2)]: # transient states
P[i, i+d] = 1/2
print(f"Transition matrix P = \n {P}")
# Transition matrix P =
# [[1. 0. 0. 0. 0. ]
# [0.5 0. 0.5 0. 0. ]
# [0. 0.5 0. 0.5 0. ]
# [0. 0. 0.5 0. 0.5]
# [0. 0. 0. 0. 1. ]]
ここからWikipediaの以下の式を使ってゴールまでの期待値を求めていきます。
吸収マルコフ連鎖の標準形(Canonical form)にする
要は$P$の行と列の順番を変えて遷移状態を左上$Q$に、吸収状態を右下$I$に持ってきます。吸収マルコフ連鎖の定義から分かるように$I$は単位行列です
\begin{align}
\large Pc
= \begin{pmatrix} Q & R\\\ 0 & I \end{pmatrix}
\end{align}
多分例を見たほうがわかりやすいと思うのでプログラムで例題１のPからQとRを求めます。
Qを求める
まず遷移状態は$1,2,3$なのでこの行と列だけ（真ん中の３行３列）を抜き出します。
Q = P[1:-1,:][:,1:-1] # remove 0,n row/coloumn
print(f"Q = {Q}")
# Q = [[0. 0.5 0. ]
# [0.5 0. 0.5]
# [0. 0.5 0. ]]
Rを求める
Rは遷移状態$\rightarrow$吸収状態の確率行列なので行は$1,2,3$、列は$0,5$を抜き出します
R = P[1:-1,:][:,[0,-1]]
print(f"R = {R}")
# R = [[0.5 0. ]
# [0. 0. ]
# [0. 0.5]]
遷移状態にいる回数の期待値N
$Q$を用いて遷移状態にいる回数の期待値の行列$N$は
\begin{align}
\large N := \sum_{k=0}^{\infty}Q^k=(I-Q)^{-1}
\end{align}
で求まるということなので、
I = np.identity(n-1)
N = np.linalg.inv(I-Q)
print(f"N = {N}")
# N = [[1.5 1. 0.5]
# [1. 2. 1. ]
# [0.5 1. 1.5]]
吸収状態に到達するまでの回数の期待値t
この$N$の各行の和を求めればよいのですが、行列的に書くと、
\begin{align}
&\large t := N1 \\
&1はすべて1のベクトル
\end{align}
t = N@np.ones(n-1)
print(f"t = N1 = {t}")
# t = N1 = [3. 4. 3.]
従って例題１の答えは、1から吸収状態に到達するまでの回数の期待値なので3回となります
どの吸収状態に行くかの確率B
\begin{align}
&\large B := NR \\
\end{align}
B = N@R
print(f"B = NR = {B}")
# B = NR = [[0.75 0.25]
# [0.5 0.5 ]
# [0.25 0.75]]
これは$1$からスタートしたら$0$に到達する確率が$0.75$、$4$に到達する確率が$0.25$であることを表しています。
今回はマルコフ連鎖の遷移行列が与えられたときに、吸収状態に到達するまでの回数の期待値を求める方法を紹介しました。その手順をまとめると、
- 確率行列（マルコフ連鎖の遷移行列）を求める
- 吸収マルコフ連鎖であることを確認して、遷移状態と吸収状態を識別する
- 吸収マルコフ連鎖の標準形の$Q,R$を求める
- これから公式を用いて$N, t$を求める
- tが各状態からのゴールへのステップ数の期待値
【例題 ２】チェイス・ゲーム
偶数のn人が輪になって座り、まず向かい合っている２人がサイコロを各々１個ずつ持って振ります。もし１なら左の人に、６なら右の人にサイコロを渡し、それ以外ならそのまま持ちます。ゲームは誰かがサイコロを２個持ったら負けで終わります。これを6人で行ったときの終わるまでの回数の期待値を求めよ
まず２つのサイコロの持っている人の位置の差（0-5）を状態として考えます。この差が０になるとゲームは終わるので、0が吸収状態、1-5が遷移状態です
２つのサイコロを振った結果位置の差がの変化分$(-2,-1,0,1,2)$の各々の確率は
dice2 =[(-2,1/36),(-1,8/36),(0,18/36),(1,8/36),(2,1/36)]
これを使って確率行列を求めます
n = 6
P = np.zeros((n,n)) # transition matrix
for i in range(n):
for d, p in dice2:
P[i][(i+d)%n] += p
0が吸収状態なので0列0行を除いてQを求めます
Q = P[1:,:][:,1:] # delete row 0, column 0 (0 is absorbing state)
print(Q)
# [[0.5 0.22222222 0.02777778 0. 0.02777778]
# [0.22222222 0.5 0.22222222 0.02777778 0. ]
# [0.02777778 0.22222222 0.5 0.22222222 0.02777778]
# [0. 0.02777778 0.22222222 0.5 0.22222222]
# [0.02777778 0. 0.02777778 0.22222222 0.5 ]]
公式に従って$N, t$を求めます
N = np.linalg.inv(np.identity(n-1)-Q) # N = (I-Q)^(-1)
t = N@np.ones(n-1) # t = N1
print(t)
# [ 9.52272727 13.81818182 15.34090909 13.81818182 9.52272727]
サイコロの位置の差の初期値は2なので求める回数の期待値は15.34です
（開発環境：Google Colab）
この考え方はProject Euler, Problem 227 The Chaseを解くのに役に立ちます