マルコフ連鎖の例に挙げられる明日の天気の実際の確率を調べてみた

マルコフ連鎖でどうこうする記事ではありません

1. はじめに

マルコフ連鎖についての例でよく今日の天気によってのみ明日の天気が決まるというものがあります。
はっきり言って全然現実世界に即していない例で、もはやマルコフ連鎖が使える場所など存在しないと言っているようなものではないかとすら思うのですが、ここで一つ疑問がわきます。
別に確率だけなら求められないか？
ということです。

現実世界において、今日の天気によってのみで明日の天気が決まるということはあり得ませんが、過去のデータを分析することによって、現在に至るまでの確率の値は求めることができるのではないでしょうか。

という訳で、気象庁から過去のデータを所得し、状態遷移図を作成するといったところまでやっていきたいと思います。

ちなみにgithubにて今回のコードのipynbファイルを公開しています。

2. 過去データのスクレイピングと変更

下記リンクの記事が大変参考になりました。
GoogleColaboratoryで気象庁の過去気象データをスクレイピングしてみた。

まずBeautifulSoup等、必要なライブラリをインポートします。

import requests
from bs4 import BeautifulSoup
import csv

place_codeA = [44] #都道府県コード
place_codeB = [47662] #地域コード
place_name = ["東京"]

気象庁のurlを書きます

# URL
base_url = "http://www.data.jma.go.jp/obd/stats/etrn/view/daily_s1.php?prec_no=%s&block_no=%s&year=%s&month=%s&day=1&view=p1"

さて、このまま気象庁のデータを使用しますと画像のようになります

晴後一時曇などよく分かりません。
晴→晴→晴→晴→曇→晴→曇
みたいな感じでしょうか。
このまま状態遷移図を作成しますと滅茶苦茶になりますので晴、曇、雨、雪の4つのデータにします。
具体的には、晴後曇の場合は後以降を無視します。

そのための準備が下記のコードです。

unused_weather=["一時","時々","後","、"]
weather_list=["晴","曇","雨","雪"]
def change_weather(weather_conditions):
  txt=str(weather_conditions)
  
  for i in range(4):
    if(txt.find(unused_weather[i])>=0):
      txt=txt[:txt.find(unused_weather[i])]
  for i in range(4):
    if(txt.find(weather_list[i])>=0):
      txt=weather_list[i]
  if(txt.find("みぞれ")>=0):
      txt=weather_list[3]
  
  return txt

次に、天気概要と日付（必須ではない）を気象庁のサイトからスクレイピングしcsvファイルにします。
今回は2020年1月1日から、昨日(2022年6月2日)までの日付をスクレイピングしました。

for place in place_name:
  All_list = [["年月日","天気（昼）","天気（夜）"]]
  index = place_name.index(place)

  weather_set={"晴"}
  # for文で該当期間抽出
  for year in range(2020,2023):
    for month in range(1,13):
      r = requests.get(base_url%(place_codeA[index], place_codeB[index], year, month))
      r.encoding = r.apparent_encoding

      soup = BeautifulSoup(r.text)
      rows = soup.findAll('tr',class_='mtx')
      rows = rows[4:]

      for row in rows:
        # trのなかのtdをすべて抜き出す
        data = row.findAll('td')


        rowData = [] #初期化

        #日付を所得
        rowData.append(str(year) + "/" + str(month) + "/" + str(data[0].string))

        #天気概況を所得
        weather_set.add(change_weather(data[19].string))
        weather_set.add(change_weather(data[20].string))
        rowData.append(change_weather(data[19].string)) 
        rowData.append(change_weather(data[20].string)) 
        All_list.append(rowData)

  #東京.csvを作成
  with open('東京.csv', 'w',encoding="utf_8_sig") as file:
    writer = csv.writer(file, lineterminator='\n')
    writer.writerows(All_list)
  print(weather_set)

以上のコードを実行することにより

このようなファイルを作成することができました。いえい

3. マルコフ連鎖から状態遷移図を作成する

今回は天気（昼）のみを用います。
それなら天気（昼）を一次元配列に突っ込んでいけばよかっただろうという人がいますでしょうがおっしゃる通りです。

昨日の天気をi、今日の天気をjとし、二次元配列をインクリメントします。
その後、各天気ごとの要素数で割ることにより、行を足すと1になる二次元配列を作成することができます。

import numpy as np
with open('東京.csv') as f:
  reader = csv.reader(f)
  data_count=0
  #晴=0,曇=1,雨=2,雪=3とした二次元配列を作成する
  weather_markov=np.zeros((4, 4))
  tmp=-1
  for row in reader:
    data_count+=1
    if(row[1]=='None'):
      break
    weather_now=-1
    #キモイ変換
    if(row[1]=='晴'):
      weather_now=0
    if(row[1]=='曇'):
      weather_now=1
    if(row[1]=='雨'):
      weather_now=2
    if(row[1]=='雪'):
      weather_now=3
    if(tmp!=-1):
      weather_markov[tmp][weather_now]+=1
    if(weather_now!=-1):
      tmp=weather_now
  for i in range(4):
    sum=0
    for j in range(4):
      sum+=weather_markov[i][j]
    for j in range(4):
      weather_markov[i][j]/=sum
    
  print(data_count)
  print(weather_markov)

データ数と二次元配列の値は以下の画像の通りとなりました。
なんとここ3年で東京の昼、雪がメインだった日は1日だけだそうです。
ちなみに勝手にみぞれも雪に含んでいるので純粋な雪のみだった日はは0日かもしれませんね。
というわけで雪は除きます。

最後に状態遷移図を作成します。
下記リンクの記事が大変参考になりました。
Pythonの状態遷移パッケージ(transitions)で状態遷移表を自動生成する

!pip install transitions
from transitions import Machine
from transitions.extensions import GraphMachine

states = ['Sunny', 'Cloudy', 'Rain']#,'Snow']

transitions = []
for i in range(3):
    for j in range(3):
      d={}
      d['trigger']=str(round(weather_markov[i][j], 2))
      d['source']=states[i]
      d['dest']=states[j]
      transitions.append(d)
print(transitions)

class Model(object):
    pass
model = Model()
machine = GraphMachine(model = model, states=states, transitions=transitions, initial=states[0], auto_transitions=False)
model.get_graph().draw('transitions.png', prog='dot', format='png') #ファイルで出力

これを実行しますと

はい、いい感じですね。
もし今日晴れなら、マルコフ連鎖的には明日は雨ではないと思ってもいいでしょうね。