Go to Qiita Advent Calendar 2024 Top
LoginSignup
9
7

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?

More than 3 years have passed since last update.

@TeRa_YUKI(寺。)

Pythonで書くマルコフ連鎖の遷移確率

Last updated at Posted at 2019-12-23

1. はじめに

マルコフ連鎖の遷移確率についての記事を書きます。

具体的には、時系列データから遷移確率を求め、行列にする関数を紹介します。
マルコフ連鎖の遷移確率行列をどうこうする関数の紹介ではないことに注意してください。

「2. マルコフ連鎖とは」から細かい情報を書きますが、
コードだけ知りたいという方は、「3. 遷移確率を求める関数」から参照していただければ十分かなと思います。

2. マルコフ連鎖とは

2.1. 遷移確率の定義

下記は、Wikipediaに書いてあるマルコフ連鎖の定義です。

一連の確率変数 X1, X2, X3, ... で、現在の状態が決まっていれば、過去および未来の状態は独立であるものである。

引用元: https://ja.wikipedia.org/wiki/%E3%83%9E%E3%83%AB%E3%82%B3%E3%83%95%E9%80%A3%E9%8E%96

2.2. マルコフ連鎖の遷移確率の式

下記は、Wikipediaに書いてあるマルコフ連鎖の式です。

Pr(X_{n+1}=x | X_n=x_n,..., X_1=x_1,X_0=x_0) = Pr(X_{n+1}=x_{n+1} | X_n=x_n)

Xi のとりうる値は、連鎖の状態空間と呼ばれ、可算集合S をなす。マルコフ連鎖は有向グラフで表現され、エッジにはある状態から他の状態へ遷移する確率を表示する。
マルコフ連鎖の一例に有限状態機械がある。これは、時刻n において状態 y にあるとすると、それが時刻n + 1 において状態x に動く確率は、現在の状態にだけ依存し、時刻n には依存しない。

引用元: https://ja.wikipedia.org/wiki/%E3%83%9E%E3%83%AB%E3%82%B3%E3%83%95%E9%80%A3%E9%8E%96

下記の式で表される確率を遷移確率と呼びます。

Pr(X_{n+1}| X_n)

引用元:
https://mathtrain.jp/markovchain

2.3. 具体例

実際に遷移確率を求めるわけですが、定義や式だけではイメージしづらいかと思います。
ここでは、天気を例に考えてみることにします。

S = {晴れ=0, 曇=1, 雨=2}という1か月間の天気の状態があったとします。

tenki = np.array([0, 0, 2, 0, 0, 0, 0, 0, 1, 0, 0, 2, 0, 1, 1, 0, 1, 1, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1])
print('晴れた日数:{0}日, 曇りの日数:{1}日, 雨の日数:{2}日'.format(np.count_nonzero(data==0), np.count_nonzero(data==1), np.count_nonzero(data==2))))
出力
晴れた日数:20日, 曇の日数:8日, 雨の日数:2日

上記をみるとどうやら晴れが一番多く、次いで曇りそして雨が多いということがぱっと見で分かるかと思います。しかし、**“天気が遷移した確率”**まではわかりません。そこで、遷移確率がわかった状態(図1)の行列が欲しいと考えるわけですが、検索しても出てきません。ということで、関数を作っていきたいと思います。

3. 遷移確率を求める関数

3.1. pythonのバージョン・環境

Python 3.7.1を使用しています。
使用環境はJupyter notebookです。

3.2. 関数の作成

3.2.1 必要なライブラリのimport

import numpy as np
import copy
import itertools
import seaborn as sns

seabornは関数的には入れなくてもよいですが、後で可視化するので入れてます。

3.2.2 データ

今回は先ほどのtenkiを使用します。

3.2.3 関数とその解説

def tp(transition_probability):
    
    data = transition_probability
    zero = np.zeros((np.max(data)+1,np.max(data)+1))

    for i in range(len(data)-1):
        j = copy.deepcopy(i)
        j += 1
        for x, y in itertools.product(range(np.max(data)+1), range(np.max(data)+1)):
            if data[i] == x and data[j] == y:
                zero[x][y] += 1
    
    row_sum = np.sum(zero, axis=1).reshape((np.max(data)+1,1))
    prob    = zero / row_sum
    
    return prob

  • zero

    • tenkiで使用されている状態は「0～2」の3つです。
    • そのため、「np.max(data)+1」をとれば最大値*最大値の0行列が作成できます。
    • (3つ以上でも以下でも対応できます。)
      　

  • forの中の処理

    • 「i, j」はそれぞれ Xn と Xn+1 を表しています。
    • 「x, y」はそれぞれ変数「zero」の列と行を表しています。
    • ネストになっているforの中の処理によって、変数「zero」に具体的な組み合わせが入ります。

  • row_sum

    • 変数「zero」の中はカウントデータです。
    • 各行ごとにカウントの合計を出しています。

  • prob

    • 各行の合計と各行の各個別の値を除算しています。

3.2.3 結果

これにて遷移確率の行列が求められました。

print(tp(tenki))
出力
[[0.65       0.25       0.1       ]
 [0.57142857 0.42857143 0.        ]
 [1.         0.         0.        ]]

3.2.4 可視化

これだけでは少々物悲しいので、可視化してみます。

sns.heatmap(tp(data), cmap='Blues', vmin=0, vmax=1, center=.5,
            square=True, cbar_kws={"shrink": .5},
            xticklabels = 1, yticklabels = 1)

ダウンロード.png

このヒートマップは左の軸が Xn を表していて、下の軸が Xn+1 を表しています。

5. 追記

5.1. Graphvizを使って可視化する

Graphvizを使用するとヒートマップよりもわかりやすく描画できるので、Graphvizを使った方法も紹介します。しかし、Graphvizはくせ者で、インストールしただけでは使えず、パスをいじらないといけません。 こちら を参考にすると良いと思います。

5.1.1. 関数とその解説

では、さっそくコーディングしていきます。

from graphviz import Digraph

def Graphviz(data, node_label):
    states = np.max(data)+1

    g = Digraph()

    for i in range(state):
        g.node(str(i), label=node_label[i])

    edges = np.array([np.sort(np.array([np.arange(states)]*states).flatten()),
                      np.array([np.arange(states)]*states).flatten()]).T

    edge_labels = np.round(tp(data), 2).flatten().astype('str')

    for i, e in enumerate(edges):
        if edge_labels[i] != '0.0':
            g.edge(str(e[0]), str(e[1]), label=edge_labels[i])
            
    return g

  • states

    • tenkiで使用されている状態は「0～2」の3つです。
    • そのため、「np.max(data)+1」にすればすべての状態を取得できます。
      　

  • 「for i in range(state)」の処理

    • 「i」分のnodeを作成し、ラベルをつけています。

  • edges

    • 状態は「states**2」で表すことができます。
    • すべてのnodeからstates分のedgeが伸びるような行列を作成しています。

  • edge_labels

    • 遷移確率を1次元配列にして、str型にしています。

  • 「for i, e in enumerate(edges)」の処理

    • edgeを作成し、ラベル(確率)をつけています。
    • 確率が0の場合はマスクしています。

5.1.2. 結果

Graphvizを使って可視化すると下記のようになります。他にも文字や枠に色をつけたりnodeの形を変えたりできます。 こちら に詳しく載っています。

node_label = ['晴れ', '曇り', '雨']
Graphviz(tenki, node_label)

4. 最後に

マルコフ連鎖の遷移確率についての関数を書いてみましたが、いかかだったでしょうか。
分かりづらい点や間違っている点があれば、ご指摘いただけると幸いです。

参考URL

・マルコフ連鎖(Wikipedia)
https://ja.wikipedia.org/wiki/%E3%83%9E%E3%83%AB%E3%82%B3%E3%83%95%E9%80%A3%E9%8E%96

・マルコフ連鎖の基本とコルモゴロフ方程式(高校数学の美しい物語)
https://mathtrain.jp/markovchain

・Graphviz をインストールする
http://ruby.kyoto-wu.ac.jp/info-com/Softwares/Graphviz/

・Graphviz
https://graphviz.readthedocs.io/en/stable/

9
7
4

Register as a new user and use Qiita more conveniently

  1. You get articles that match your needs
  2. You can efficiently read back useful information
  3. You can use dark theme
What you can do with signing up
Sign upLogin
9
7

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?