マルコフ連鎖の定常分布・パラメータ推定|問題演習で理解する統計学【14】

下記などで取り扱った、マルコフ連鎖の定常分布・パラメータ推定に関する問題演習を通した理解ができるように問題・解答・解説をそれぞれ作成しました。
https://www.hello-statisticians.com/explain-terms-cat/markov_chain1.html

基本問題

マルコフ連鎖の定義

・問題
$$
\begin{align}
P(X_{n+1}|X_{n}, X_{n-1}, … , X_{1}, X_{0}) = P(X_{n+1}|X_{n})
\end{align}
$$
マルコフ連鎖の定義にあたっては基本的に上記の式で定義される。この式は一度理解すればシンプルで見やすい一方で、確率の表記に慣れていないとイメージがつきにくい。
よって、マルコフ連鎖の定義式のイメージがつくように演習問題を作成することとした。

下記の問題に答えよ。
i) 全事象が$U={1,2,3,4,5,6,7,8,9,10}$、事象$A$が$A={1,5,7}$のように表されるとき、確率$P(A)$を求めよ。ただし、それぞれの値が観測される確率は同様に確からしいとする。
ⅱ) 0より大きい整数で表される確率変数$X$に関し、$P(X=1)=0.2, P(X=2)=0.3$のとき、$P(X \geq 3)$を求めよ。
ⅲ) 事象Bが起こった前提での事象Aの起こる確率を条件付き確率といい、$P(A|B)$のように表す。条件付き確率$P(A|B)$をBが起こる確率の$P(B)$とAとBが同時に起こる確率の$P(A \cap B)$を用いて表せ。
iv) 季節をS、その日の最高気温をT、場所をPのように表す際に、下記の確率がどのようになり得るかについて考えよ。
$$
\begin{align}
& P(T \geq 15|S=夏、P=東京) \\
& P(T \geq 15|S=冬、P=東京) \\
& P(T \geq 15|S=冬、P=沖縄)
\end{align}
$$
v) マルコフ連鎖はⅲ)、iv)で取り扱った条件付き確率を用いて定義するが、冒頭の式を条件付き確率の理解に基づいて解釈せよ。

・解答
i)
求める確率$P(A)$は下記のようになる。
$$
\large
\begin{align}
P(A) &= \frac{3}{10} \\
&= 0.3
\end{align}
$$

ⅱ)
求める確率$P(X \geq 3)$は下記のようになる。
$$
\large
\begin{align}
P(X \geq 3) &= 1 – (P(X=1) + P(X=2)) \\
&= 1 – (0.2 + 0.3) \\
&= 0.5
\end{align}
$$

ⅲ)
条件付き確率$P(A|B)$は$P(B)$と$P(A \cap B)$を用いて下記のように表すことができる。
$$
\large
\begin{align}
P(A|B) = \frac{P(A \cap B)}{P(B)}
\end{align}
$$

iv)
東京の夏の最高気温は基本的に15度を超えることが多いので、$P(T \geq 15|S=夏、P=東京)$は1に近い値になる。また、冬の最高気温が15度を超えることは少ないため$P(T \geq 15|S=冬、P=東京)$は0に近い値となる。
一方で東京より赤道に近い沖縄では冬もそれなりに暖かいため、$P(T \geq 15|S=冬、P=沖縄)$は1に近い値となる。

v)
マルコフ連鎖は条件付き確率の形式で表されるが、「直前の値に基づいてのみ次の値が決まる」と解釈することができる。

・解説
マルコフ連鎖の理解にあたって数式が抽象的でわかりにくい場合もあると思われたため、数式の表記になれることができるように作成を行いました。特にⅲ)〜ⅴ)で取り扱った条件付き確率については様々なトピックで出てくるので、式の定義がわからないと全くわからないになりがちです。
また、ここで取り扱ったマルコフ連鎖の定義を拡張して強化学習の数式なども表されたりするので、抑えておくことで様々な概念の理解に役立つと思います。

状態確率ベクトル・確率行列・定常分布

・問題
マルコフ連鎖を考えるにあたって主要トピックとなるのは状態確率ベクトル(state probability vector)と確率行列(probability matrix)である。
ここでマルコフ連鎖は「直前の値に基づいてのみ次の値が決まる」ことを定義しているため、状態を確率的に表したベクトルに確率行列を作用させることで具現化することができるというのは抑えておくと良い。
状態確率ベクトルが$N$個の状態を持つと考えた際に、時点$n$における状態確率ベクトル$\pi_n$と確率行列$Q$を下記のように定義する。
$$
\begin{align}
\pi_n &= \left(\begin{array}{ccc} p_n(1) & … & p_n(N) \end{array} \right) \\
Q &= \left(\begin{array}{ccc} p(1,1) & … & p(1,N) \\ … & … & … \\ p(N,1) & … & p(N,N) \end{array} \right)
\end{align}
$$
このとき下記の問題に答えよ。

i) 時点$0$が初期状態を表し、状態の数が$N=3$かつ状態2の確率が0.8、状態3の確率が0.2のとき、初期分布$\pi_0$を求めよ。「時点0における状態確率ベクトル=初期分布」であることは前提とした。
ⅱ) i)のように初期分布$\pi_0$が与えられ、確率行列$Q$が下記のように与えられるとき、時点1における状態確率ベクトルの$\pi_1=\pi_0 Q$を求めよ。
$$
\begin{align}
Q &= \left(\begin{array}{ccc} 0.5 & 0.3 & 0.2 \\ 0.2 & 0.6 & 0.2 \\ 0.2 & 0.3 & 0.5 \end{array} \right)
\end{align}
$$
ⅲ) $\pi=\pi Q$が成立する$\pi$を定常分布であると考えるとき、ⅱ)で定義した$Q$に対応する定常分布の$\pi$を求めよ。

・解答
i)
初期分布$\pi_0$は下記のように表すことができる。
$$
\large
\begin{align}
\pi_n &= \left(\begin{array}{ccc} 0 & 0.8 & 0.2 \end{array} \right)
\end{align}
$$

ⅱ)
時点1における状態確率ベクトルの$\pi_1=\pi_0 Q$を計算すると下記のようになる。
$$
\large
\begin{align}
\pi_1 &= \pi_0 Q \\
&= \left(\begin{array}{ccc} 0 & 0.8 & 0.2 \end{array} \right) \left(\begin{array}{ccc} 0.5 & 0.3 & 0.2 \\ 0.2 & 0.6 & 0.2 \\ 0.2 & 0.3 & 0.5 \end{array} \right) \\
&= \left(\begin{array}{ccc} 0.16+0.04 & 0.48+0.06 & 0.16+0.1 \end{array} \right) \\
&= \left(\begin{array}{ccc} 0.2 & 0.54 & 0.26 \end{array} \right)
\end{align}
$$

ⅲ)
$\pi = \left(\begin{array}{ccc} a & b & c \end{array} \right)$とおき、$\pi=\pi Q$に代入する。
$$
\large
\begin{align}
\pi &= \pi Q \\
\left(\begin{array}{ccc} a & b & c \end{array} \right) &= \left(\begin{array}{ccc} a & b & c \end{array} \right) \left(\begin{array}{ccc} 0.5 & 0.3 & 0.2 \\ 0.2 & 0.6 & 0.2 \\ 0.2 & 0.3 & 0.5 \end{array} \right) \\
\left(\begin{array}{ccc} a & b & c \end{array} \right) &= \left(\begin{array}{ccc} 0.5a+0.2b+0.2c & 0.3a+0.6b+0.3c & 0.2a+0.2b+0.5c \end{array} \right)
\end{align}
$$
上記より、$b=1.5a=1.5c$が得られる。$a+b+c=1$に代入することで下記を求めることができる。
$$
\large
\begin{align}
a &= \frac{2}{7} \\
b &= \frac{3}{7} \\
c &= \frac{2}{7}
\end{align}
$$
上記より、定常分布の$\pi$は下記のように表せる。
$$
\large
\begin{align}
\pi = \left(\begin{array}{ccc} \frac{2}{7} & \frac{3}{7} & \frac{2}{7} \end{array} \right)
\end{align}
$$

・解説
状態確率ベクトル・確率行列・定常分布について取り扱いました。確率行列については行と列の定義がわかりにくいですが、$p(i,j)$とあれば「$i$から$j$に推移する確率」と解釈するのが良いです。
この行列の設定については転置行列を考えて状態確率ベクトルの左からかけるという定義の仕方も可能なように思われますが、$\pi_1=\pi_0 Q$が成立するように定義されているというのもあるので、基本的に教科書などの定義に合わせて表記を抑えておくと良いと思います。
$(\pi_n Q)^T = Q^{T} \pi_n^{T}$のようにも考えられるとは思いますが、表記は統一しておく方が良さそうです。

発展問題

定常分布の必要十分条件

マルコフ連鎖のパラメータ推定

・問題
マルコフ連鎖のパラメータ推定を考えるにあたっては基本的に「最尤法」と同様に考えれば良い。最尤法は同時確率を尤度と読み替え、尤度を最大にするパラメータを求めるという手法である。確率変数列${X_n}$に関し、観測列${x_n}$が与えられたとき、同時確率は下記のようになる。
$$
\large
\begin{align}
P(X_0=x_0, X_1=x_1, …, X_n=x_n) = p(x_0) \prod_{k=1}^{n} p_{\theta}(x_{k-1}, x_{k}) \qquad (1)
\end{align}
$$
上記において、$p_{\theta}(x_{k-1}, x_{k})$を考えたが、$p_{\theta}(i, j)$を「状態$i$から状態$j$に移る確率をパラメータ$\theta$を用いて表したもの」と考えるものとする。
下記の問題に答えよ。

i) 冒頭で表した式$(1)$が成立するにあたってはマルコフ性が成立していなければならない。このことを下記のマルコフ連鎖の定義式を用いて説明せよ。
$$
\begin{align}
P(X_{n+1}|X_{n}, X_{n-1}, … , X_{1}, X_{0}) = P(X_{n+1}|X_{n})
\end{align}
$$
ⅱ) 状態の数を$N=3$、確率行列$Q_{\theta}$が下記のように表すことができるとする。
$$
\begin{align}
Q_{\theta} &= \left(\begin{array}{ccc} 0.5 & \theta & 0.5-\theta \\ 0.3-\theta & 0.7 & \theta \\ 0.1 & 0.2 & 0.7 \end{array} \right)
\end{align}
$$
このとき$p_{\theta}(1,2), p_{\theta}(1,3), p_{\theta}(2,3)$をそれぞれ$\theta$を用いて表せ。
ⅲ) 観測列$\{x_n\}$が$\{1,1,2,3,2,1\}$のように観測される確率$P(X_0=1,X_1=1,X_2=2,X_3=3,X_4=2,X_5=1)$を$\theta$を用いて表せ。
iv) 尤度を$L(\theta)=P(X_0=1,X_1=1,X_2=2,X_3=3,X_4=2,X_5=1)$のようにおくとき、$L(\theta)$を最大にする$\theta$と$\log{L(\theta)}$を最大にする$\theta$は同じであることを説明せよ。
v) $L(\theta)$を最大にする$\theta$を求めよ。

・解答
i)
同時確率$P(X_0,X_1,X_2)$を考えるとき、条件付き確率の公式に基づいて下記のように表すことができる。
$$
\large
\begin{align}
P(X_0,X_1,X_2) &= P(X_1,X_2|X_0)P(X_0) \\
&= P(X_2|X_1,X_0)P(X_1|X_0)P(X_0)
\end{align}
$$
上記の$P(X_2|X_1,X_0)$はマルコフ連鎖の定義式より、$P(X_2|X_1,X_0)=P(X_2|X_1)$のように変形することができるので、$P(X_0,X_1,X_2)$は下記のように変形できる。
$$
\large
\begin{align}
P(X_0,X_1,X_2) &= P(X_2|X_1,X_0)P(X_1|X_0)P(X_0) \\
&= P(X_2|X_1)P(X_1|X_0)P(X_0) \\
&= P(X_0) \prod_{k=1}^{2} p_{\theta}(X_{k-1}, X_{k})
\end{align}
$$
$P(X_2|X_1,X_0)=P(X_2|X_1)$がマルコフ性に基づいており、同様に$P(X_0,X_1,X_2,…,X_n)$についても考えることができる。

ⅱ)
確率行列$Q_{\theta}$より、$p_{\theta}(1,2), p_{\theta}(1,3), p_{\theta}(2,3)$はそれぞれ下記のように表せる。
$$
\large
\begin{align}
p_{\theta}(1,2) &= \theta \\
p_{\theta}(1,3) &= 0.5-\theta \\
p_{\theta}(2,3) &= \theta
\end{align}
$$

ⅲ)
確率行列$Q_{\theta}$より、$P(X_0=1,X_1=1,X_2=2,X_3=3,X_4=2,X_5=1)$は下記のように計算できる。
$$
\large
\begin{align}
P(X_0=1,X_1=1, &X_2=2,X_3=3,X_4=2,X_5=1) \\
&= 0.5 \times \theta \times \theta \times 0.2 \times (0.3-\theta) \\
&= 0.1 \theta^2 (0.3-\theta)
\end{align}
$$

iv)
対数関数$\log{x}$は定義域$x>0$において単調増加関数であるので、$L(\theta)$を最大にする$\theta$と$\log{L(\theta)}$を最大にする$\theta$は同じとなる。

v)
iv)より、$L(\theta)$を最大にする$\theta$を求めるにあたっては、$\log{L(\theta)}$を考えればよい。
$$
\large
\begin{align}
\log{L(\theta)} &= \log{0.1 \theta^2 (0.3-\theta)} \\
&= 2 \log{\theta} + \log{(0.3-\theta)} + Const
\end{align}
$$
上記を$\theta$に関して微分し、微分した値が0となるような$\theta$を求める。
$$
\large
\begin{align}
\frac{\partial \log{L(\theta)}}{\partial \theta} &= 0 \\
\frac{2}{\theta} – \frac{1}{0.3-\theta} &= 0 \\
\frac{2}{\theta} &= \frac{1}{0.3-\theta} \\
2(0.3-\theta) &= \theta \\
0.6 &= 3 \theta \\
\theta &= 0.2
\end{align}
$$
上記より、$L(\theta)$を最大にする$\theta$は$\theta=0.2$となる。

マルコフ連鎖と強化学習

参考書籍

・統計学実践ワークブック