マルコフ過程 - 機械学習に詳しくなりたいブログ

マルコフ過程とマルコフ連鎖

時刻 $n$ における状態を $x_n$ と表し、 $x_{n+1}$ の確率が $x_1,\cdots,x_n$ から決まるのではなく、１つ前の状態 $x_n$ のみによって定まるときにこれをマルコフ性といい、このような性質を持つ確率過程をマルコフ過程といいます。状態が離散的なとき、特にマルコフ連鎖と呼びます。

状態が３つあるとしたとき、例えば以下の図のようにマルコフ連鎖を表すことができます。

各状態をつなぐ矢印は、時刻が進んだときに遷移できる方向を表し、数字は遷移する確率を表します。

時刻 $n$ に状態 $i$ にいる確率を $p_{i}^{(n)}$ とすれば、行列を用いて

$\displaystyle \left( \begin{array}{c} p_{1}^{(n+1)} \\ p_{2}^{(n+1)} \\ p_{3}^{(n+1)} \end{array} \right) = \left( \begin{array}{ccc} 0 & \frac{1}{4} & \frac{1}{2} \\ \frac{1}{4} & \frac{3}{4} & \frac{1}{2} \\ \frac{3}{4} & 0 & 0 \\ \end{array} \right) \left( \begin{array}{c} p_{1}^{(n)} \\ p_{2}^{(n)} \\ p_{3}^{(n)} \end{array} \right) \tag{1}$

と表すことができます。この式における行列のように、マルコフ連鎖を規定する確率は推移核と呼ばれます。

エルゴード性と定常分布

マルコフ連鎖において、

どの状態からスタートしても有限時間内に任意の状態に遷移可能（既約性）
周期的でない

であることをエルゴード性といい、これを満たすとき、時刻 $n \to \infty$ で

$\mathbf p = \mathbf P \mathbf p \tag{2}$

を満たす分布 $\mathbf p$ に収束します。これを不変分布または定常分布と呼びます。冒頭の例はエルゴード性を満たしていますから、ある確率に収束するはずです。後に掲載しているコードで試した結果、 0.23529412 0.58823529 0.17647059 に収束しました。つまりこれを式(2)のように表現すれば、

$\displaystyle \left( \begin{array}{c} 0.23529412\\ 0.58823529\\ 0.17647059 \end{array} \right) = \left( \begin{array}{ccc} 0 & \frac{1}{4} & \frac{1}{2} \\ \frac{1}{4} & \frac{3}{4} & \frac{1}{2} \\ \frac{3}{4} & 0 & 0 \\ \end{array} \right) \left( \begin{array}{c} 0.23529412\\ 0.58823529\\ 0.17647059 \end{array} \right) \tag{3}$

であるということです。このマルコフ連鎖は、各状態が 0.23529412 0.58823529 0.17647059 の確率で発生する乱数生成器であるとみなすことができるようです。（ただしサンプルの前後間には相関がある）

なお次の図のようなマルコフ連鎖は既約性を満たしません。２の状態になったら他の状態へ遷移することができずに、永遠に同じ状態を繰り返してしまいます。

確率論の参考書は難しくてなかなか理解が進まず、、、。Amazonでわかりやすいと評価が高かったものを購入してみましたが、僕の頭ではとても独学で理解できる代物ではなく、今まで学んできた確率統計とは別の学問としか思えないような難解さ。とりあえずマルコフ過程については、このくらいの理解で先に進んでいこうと思います。

今回実験したコードです。状態遷移をシミュレーションしたものと、式(1)を繰り返し演算した結果がそれぞれ一致することが確認できました。コード内では初期状態を1としていますが、何を初期状態としても収束しました。

# マルコフ連鎖のテスト
import numpy as np


def markov(x):
    p = np.random.uniform()
    if x == 1:
        if p < 0.25:
            return 2
        else:
            return 3
    if x == 2:
        if p < 0.25:
            return 1
        else:
            return 2
    if x == 3:
        if p < 0.5:
            return 1
        else:
            return 2


# 初期値
x = 1

# サンプル数
N = 1000000

history = np.empty(N)

for i in range(N):
    x = markov(x)
    history[i] = x

print(np.count_nonzero(history == 1)/N)
print(np.count_nonzero(history == 2)/N)
print(np.count_nonzero(history == 3)/N)

s = np.array([[0,0.25,0.5],[0.25,0.75,0.5],[0.75,0,0]])
x = np.array([1,0,0]).reshape(3,1)

for i in range(N):
    x = np.dot(s, x)

print(x)