🤖

Pythonでマルコフチェーンを用い、天気の遷移確率を求めよう

2023/03/25に公開

Python

Streamlit

マルコフ連鎖

tech

本記事の目的

ベイズ推定記事のための前記事第二弾。前回はモンテカルロ法を用いて円周率を導出した。

このモンテカルロ法と本記事のマルコフチェーンを組み合わせた手法をMCMC法（マルコフ連鎖モンテカルロ法）といい、ベイズ推定を行う際によく用いられている。なので、ベイズ推定を学ぶ前にマルコフチェーンについても学ぶ必要があると思い、記事にまとめた。

マルコフチェーンとは

表面だけ説明していく。マルコフチェーンとは、事後データが前データの中でも直前のデータのみに依存するモデルを指す。数学的に書くなら、時刻t+1のデータ $X_{t+1}$ は $X_t$ のみに依存し、 $X_{t-1}$ などは関係ない。より数学的に表すと下記の通りになる。

P(X_{t+1} \mid X_t,X_{t-1}, \cdots , X_1) = P(X_{t+1} \mid X_t)

$P(X_{t+1} \mid X_t)$ を遷移確率といい、これを意識しながら、プログラミングしてみる。

理論：天気の遷移確率の導出

パッと思い浮かんだのが、「前日の天気が今日の天気にどのくらい関係しているか」だったので、これを行っていく。こちらのサイトを頼らせていただいた。

天気は簡潔にするため、0:晴、1:曇、2:雨の3通りとし、1日目の天気を配列の0番目、2日目の天気を配列の1番目、･･･、としていって遷移確率を求める。また、遷移確率と難しく書いているが(前日の天気,今日の天気)の各確率、つまり(0,0),(0,1),(0,2),(1,0),(1,1),(1,2),(2,0),(2,1),(2,2)の9通りの確率を求めているに過ぎない。

実践：天気の遷移確率の導出

早速、コードを記す。

trans_prob.py

import numpy as np
import itertools

weather=np.array([1, 1, 2, 1, 0, 1, 2, 0, 1, 1, 0, 2, 0, 0, 2, 0, 0, 0, 0,
                0, 0, 0, 1, 0, 0, 2, 1, 0, 1, 0, 0, 1, 0, 1, 1, 0, 2, 1,
                0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 1, 1, 2, 2, 2, 1, 1, 1, 1, 
                0, 0, 1, 1, 0, 0, 0, 0, 1, 0, 0, 0, 1, 1, 1, 0, 0, 0, 1, 
                1, 2, 1, 0, 0, 1, 0, 1, 2, 2, 1, 0, 0, 1, 1, 1, 0, 0, 1, 
                1, 1, 1, 0, 0, 0, 0, 0, 2, 1, 1, 1, 0, 1, 2, 1, 1, 1, 1, 
                1, 1, 1, 0, 1, 1, 1, 1, 1, 0, 1, 0, 0, 0, 0, 1, 1, 0, 0, 
                0, 0, 1, 0, 0, 1, 0, 1, 1, 1, 1, 2, 0, 1, 2, 1, 0, 2, 1, 
                0, 1, 1, 1, 0, 1, 0, 0, 1, 1, 2, 2, 2, 0, 0, 2, 1, 1, 2, 
                1, 1, 0, 1, 0, 0, 0, 1, 1, 2, 1, 1, 1, 1, 1, 1, 0, 2, 1, 
                0, 0, 0, 0, 1, 1, 1, 2, 2, 1, 1, 1, 1, 1, 0, 0, 0, 0, 1, 
                2, 2, 0, 0, 1, 0, 0, 2, 0, 1, 0, 2, 1, 0, 2, 0, 0, 2, 2, 
                0, 0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 1, 0, 1, 0, 0, 1, 2, 1, 
                0, 2, 1, 0, 2, 2, 2, 2, 0, 1, 1, 1, 1, 1, 0, 0, 2, 1, 1, 
                1, 2, 1, 1, 1, 1, 1, 2, 0, 2, 2, 1, 1, 1, 1, 1, 1, 2, 2, 
                2, 2, 1, 1, 1, 0, 1, 0, 0, 1, 2, 1, 2, 1, 1, 1, 1, 1, 1, 
                1, 1, 1, 2, 1, 1, 1, 1, 1, 1, 0, 0, 1, 1, 1, 2, 1, 0, 2, 
                1, 1, 1, 1, 1, 2, 1, 1, 1, 1, 1, 2, 1, 0, 0, 0, 2, 0, 1, 
                0, 0, 0, 0, 1, 0, 0, 1, 1, 0, 0, 1, 1, 1, 0, 0])

def trans_prob(data):
    len_data=len(np.unique(data))
    zero = np.zeros((len_data,len_data))

    for i in range(len(data)-1):
        j=i+1
        for before_we, today_we in itertools.product(range(len_data), range(len_data)):
            if data[i] == before_we and data[j] == today_we:
                zero[before_we][today_we] += 1
    sum = np.sum(zero, axis=1).reshape((len_data,1))
    prob    = zero / sum

    return prob

print(trans_prob(weather))

weather=np.array~で何をやっているか
理論編でも述べたように、適当に0:晴、1:曇、2:雨として配列に入れている。各自、適当に変えて試してみるのもよい。

len_data=len(np.unique(data))で何をやっているか
np.unique()で重複を削除した要素配列を取得している。今回だと[0,1,2]が得られる。これにlenを付けてあげることで、要素数3が得られる。

zero = np.zeros((len_data,len_data))で何をやっているか
基本形はnumpy.zeros(shape, dtype=float, order='C')であるが第二引数と第三引数は大体省略されている。第一引数のshapeは配列の行列数を示しているが、2×3のゼロ行列を作るときはnp.zeros((2,3))のようにshape自体を()で囲むことに注意。

j=i+1で何をやっているか
iは前日の天気を表す要素の配列番号、jは当日の天気を表す要素の配列番号としているので、i+1としている。

理論編で紹介した参考サイトではj = copy.deepcopy(i)を用いていた。自分の浅い考えでは'j=i+1'しか考えつかなかったので、少し調べてみた。コピーには浅いコピーと深いコピーがあるとのこと。複合オブジェクト（リストの中にリストや辞書が存在するようなオブジェクト）をコピーするとき、浅いコピーだと、オリジナルの「オブジェクトの中のオブジェクト」のデータが変わると、その部分だけ同じようにコピー先も変わる（参照）。一方、深いコピーだと最初のコピー状態のままを維持する（コピー）。言葉だけでは分からないと思うので実例を記す。

浅いコピー

a=[1,2,[3,4]]
#[3,4]が上記の「オブジェクトの中のオブジェクト」を指す
b=a.copy()

a[1]=200
a[2][1]=400

print(a)
#[1, 200, [3, 400]]
print(b)
#[1, 2, [3, 400]]

深いコピー

import copy

a=[1,2,[3,4]]
#[3,4]が上記の「オブジェクトの中のオブジェクト」を指す
b=copy.deepcopy(a)

a[1]=200
a[2][1]=400

print(a)
#[1, 200, [3, 400]]
print(b)
#[1, 2, [3, 4]]

知らなかったことを知れてよかった。自分が扱うデータは、基本的にオリジナルが変わってもコピー先は変わってほしくないようなデータなので、これからはcopy.deepcopy()を使うことを心掛ける。

for before_we, today_we in itertools.product(range(len_data), range(len_data))で何をやっているか
ここが本記事の重要ポイント。itertools.productは直積を示している。range(len_data)は[0,1,2]を表しており、数式にすると

[0,1,2]\times[0,1,2]=[(0,0),(0,1),(0,2),(1,0),(1,1),(1,2),(2,0),(2,1),(2,2)]

となる。理論編で述べた遷移確率の各要素ですね。これをfor文で回し、if data[i] == before_we and data[j] == today_weで適当に決めた天気配列weatherを(0,0),(0,1),(0,2),(1,0),(1,1),(1,2),(2,0),(2,1),(2,2)に割り振っている。

sum = np.sum(zero, axis=1).reshape((len_data,1))で何をやっているか
上記はカウント数となっているので、確率にするためには全体数で割る必要がある。前日の天気は必ず0:晴、1:曇、2:雨のどれかになるとし、横方向（同じ行）の要素を合計するnp.sum(zero, axis=1)を用いて遷移確率を求めている。

結果

結果は次のようになると思う。

結果

[[0.51515152 0.34848485 0.13636364]
 [0.31360947 0.55029586 0.13609467]
 [0.21428571 0.51785714 0.26785714]]

分かりにくいが、行番号が前日の天気、列番号が当日の天気を表しているため、(0,1):晴→曇の遷移確率は0.34848485、(2,0):雨→晴の遷移確率は0.21428571などの結果が分かる。もっと見やすくしたい方は、ヒートマップなどを用いてみてください。本記事で伝えたかったのは、直前データから現データの確率を求める手法をマルコフチェーンといい、前記事で述べたモンテカルロ法と組み合わせて、ベイズ推定を行っていること。

おまけ：47都道府県の天気遷移確率の導出

上記では、適当にナンバリングしたデータで天気の遷移確率を求めたが、次は気象庁の実データでやってみたいと思う。ここではWEBスクレイピングを用いて、47都道府県の先月の天気遷移確率を導出してみた。ただ、WEBスクレイピングに関しては、やり方次第でサーバーへの負荷を高め、業務妨害になるなど、色々危険なので、プログラミングコードは載せず、出力画面と自分の躓きメモを記すだけにする。自分用ノートみたいなものなので、興味がある方だけ見てもらえたらと。では、出力画面をどうぞ。

出力画面（静岡県）

出力画面（和歌山県）
これらのデータは気象庁「過去の気象データ検索」を加工して作成している。いつも通りの雑なフロントに加えて、計算日の先月データということで、約30個のデータしかなかったため、遷移確率が0となっているデータがあるが、仕方ない。WEBスクレイピングを試せたので良し。
また、今回は、こちらのサイトを参考にして、作成している。
https://expy-style.net/excel/weather/

以下、自分用メモを記す。

HTMLとCSSについて
後々出てくるHTMLとCSSについてまとめておく。Web情報はHTMLとCSSでほとんど構成されている。簡潔にまとめるとHTMLは文章構造、このHTMLで表す文章の色や文字サイズなどの装飾機能をCSSが担う。Google Chromeを用いた実際のコードの確認方法は次の通り。まず、適当なウェブサイト上で右クリックすると「検証」というボタンがあると思う。これをクリックするとウィンドウが現れ、上の方にそのウェブサイトのHTMLコードが確認できると思う。また、そのHTMLコード上で右クリックし、Copy→Copy selectorでHTMLのCSSコードが確認できる。

HTML例

<a>HTMLの例です</a>

CSS例

*{color:red}

WEBスクレイピングのライブラリについて
スクレイピングはデータの取得→抽出というのが一連の流れである。スクレイピングの代表的なライブラリにはRequests、Beautiful Soup、Seleniumがあるが、Requestsは主にデータ取得、Beautiful Soupは主にデータ抽出、Seleniumはデータ取得とデータ抽出共にできるが、遅いというデメリットがある。そのため、RequestsとBeautiful Soupを用いるのが一般的である。こちらのサイトを参考にさせていただいた。

https://ai-inter1.com/beautifulsoup_1/

ライブラリBeautiful Soupの使い方について
Beautiful Soup()の第一引数はライブラリrequestsで取得した解析対象のHTML、第二引数は解析機の種類を指す。解析機には追加ライブラリ不要の"html.parser"や高速に解析できる"lxml"など様々ある。これらの情報に関しても上記のサイトを参考にするとよい。

Beautiful Soupの第一引数を(request.getで取得したHTML).textとした時、下記のエラーコードが生じた。

'<' not supported between instances of 'NoneType' and 'NoneType'

どうやらContent-typeにcharsetがない時、.textだと文字化けが発生するらしい。詳しい説明は下記のサイトを参考にしてほしい。なので、とりあえず(request.getで取得したHTML).contentでBeautiful Soup()に渡し、Beautiful Soup側でデコードしたほうが良いと思われる。

.select_oneの引数について
上記で解析するサイトと解析機の種類を決めたら、次にそのサイトのどこを解析するか細かく設定していく。その設定の仕方の一つが.select_one(CSSセレクタ)である。CSSセレクタの調べ方は上記の「HTMLとCSSについて」を参考にしてほしい。

最初、Copy→Copy selectorで「#tablefix1 > tbody > tr:nth-child(5) > td:nth-child(20)」をコピペしたが上手くいかなかった。よく見ればtbodyには:が無く、tbody要素がない。どうやらソースにtbodyが無くても、デベロッパーツールではtbodyが補完されるらしい。なので、tbodyを消したらうまくいった。

ダメな例

hour = html.select_one('#tablefix1 > tbody > tr:nth-child(5) > td:nth-child(20)').text

上手くいった例

hour = html.select_one('#tablefix1 > tr:nth-child(5) > td:nth-child(20)').text

月の日数の取得方法について
calendar.monthrange(年,月)を用いる。これを行うと、月初日の曜日（月曜が0、日曜が6）と、月の日数の情報がタプルで取得できる。なので、月の日数のみが欲しいときは最後に[1]を付ける。

例

calendar.monthrange(2022, 8)[1]

最後までお読みいただきありがとうございました。お疲れ様でした。