マルコフ連鎖の定常分布を解く

前回に引き続き、マルコフ連鎖に関する記事です。

マルコフ連鎖において、かなり大きな回数状態を更新したときに、ある時点でどんなに状態を更新しても変わらなくなることがあります。

この変わらなくなった状態のことをエルゴート状態と言います（厳密に言うと、もっと難しい）

もちろん、どんなに更新してもずっと変わり続けることもあります。例えば以下の推移行列のときは、常に更新され続きます。これはエルゴート状態とは言えません。

$$ P=\begin{bmatrix} 0 & 1 \\ 1 & 0 \end{bmatrix} $$

エルゴート状態となる例のうちの1つとして、以下のものが挙げられます。

$$ P=\begin{bmatrix} 0.4 & 0.6 & 0.0 \\ 0.2 & 0.5 & 0.3 \\ 0.1 & 0.7 & 0.2 \end{bmatrix} $$

このとき、ある状態\(i\)に到達する確率を\(\pi_i\)と表記し、\(\boldsymbol{\pi} = (\pi_0, \pi_1, …)\)と行ベクトルで表すと、以下の式が成立します。

$$ \begin{eqnarray} \left\{ \begin{array}{lll} \boldsymbol{\pi} P&=&\boldsymbol{\pi} \\ \displaystyle \sum_i \pi_i&=&1 \end{array} \right. \end{eqnarray}$$

でもこれをわざわざ手計算するのは面倒！　なので、Pythonで組みましょう。

手始め

まず、\(\boldsymbol{\pi} P =\boldsymbol{ \pi}\)について、式変形すると、\(\boldsymbol{\pi} (P – I) = \boldsymbol{0}\)となります。

ただし、ここで注意しなければならないのは、\(\det(P – I) = 0\)となっていることです。

このため、一般的な線形代数を用いた連立方程式としてまだ解けません。

ここで、\(\displaystyle \sum_i \pi_i=1 \)を使うことにします。

\( P – I \)の\(j\)列の成分を丸ごとベクトル\(\boldsymbol{\pi}\)に変更した\(Q\)を、また\(\boldsymbol{0}\)の\(j\)番目の要素を\(1\)に変更した\(\boldsymbol{v}\)を用意します。

すると、以下のように線形代数を用いた連立方程式を用意できます。

$$ \boldsymbol{\pi} Q = \boldsymbol{v} $$

エルゴート状態のとき、\(\det(Q)\neq 0\)となっているため（本当はこれも証明しないといけない）、\(\boldsymbol{\pi}\)を求めることができます。

$$ \boldsymbol{\pi} = \boldsymbol{v} Q^{-1} $$

これをPythonで表記したものが以下のようになっています。

import numpy as np








def solve(P: np.ndarray) -> np.ndarray:




    rows, cols = P.shape
    eps = 1e-8
    if rows != cols:
        raise Exception(f'P must be square. <{(rows, cols)=}')
    elif np.any(np.abs(P.sum(axis=1)-1) >= eps):
        raise Exception(f'P.sum(axis=1) must equal 1. <{P.sum(axis=1)=}>')
    elif (det := np.abs(np.linalg.det(np.identity(rows) - P))) >= eps:
        raise Exception(f'|det(I-P)| != 0. <{det=}>')




    Q = P - np.identity(rows)
    Q[:, -1] = 1
    v = np.zeros(rows)
    v[-1] = 1




    return v @ np.linalg.inv(Q)








if __name__ == '__main__':




    P = np.array([[.4,.6,.0],
                  [.2,.5,.3],
                  [.1,.7,.2]])




    π = solve(P)
    print(f'{π=}')

「SymPy使えばもっと楽にできるよ」知ってますが、重いので使いたくない。

結果は以下の通り。

π=array([0.4047619 , 0.42857143, 0.16666667])

ちゃんと算出されてる。

この結果からわかることとしては、無限回遷移させると、状態\(0, 1, 2\)へ向かう確率はそれぞれ\(17/42, 18/42, 7/42\)となっています。

さいごに

プログラムで組むと、どのように計算すれば良いかわかりやすくなり、シミュレーション論の講義の復習にもなりました。

次回は吸収マルコフ連鎖について書ければなあ、と考えています。