こんにちは、気象予報士のy_kawasaki(取得してもう早、15年経つらしい)です。
一般的な天気予報は、数値予報モデルといわれる、物理方程式で作られた仮想の地球で数値微分を行い、気温、気圧などを予測したあと、ガイダンスといわれる、機械学習を用いた、いわゆる天気予報への翻訳を行っています。
そこで、行われていることは、http://aitc.jp/events/20160916-Seika/20160916_特別講演_気象庁における機械学習の利用.pdfで説明されているように、カルマンフィルター、ロジスティック回帰、ニューラルネットワークなどが用いられています。特に、専門的に知りたい方は、気象庁数値予報課が毎年刊行している、研修テキストに詳しく載っていますので、そちらを参照してください。
さて、今回は気象庁とは別のアプローチで機械学習で、天気を予報してみようと思います。
予測するもの
翌日の東京の降水の有無
東京の前日の気象情報から翌日の降水の有無の予想
用意するもの
東京の過去の天気データ。http://www.data.jma.go.jp/gmd/risk/obsdl/index.phpこのページで適当にデータをダウンロードします。決して使いやすい形式ではないですが、CSVでダウンロードできます。
今回は、東京の2011/1/1 - 2016/12/31の日別値を手に入れます。
データ加工
ダウンロードしてきたデータを左から、「日付,平均気温,最高気温,最低気温,降水量,降水の有無,日照時間,日照時間の有無,風速,最大風速,風向,最多風向,湿度,最小湿度,雲量」と並べます。
#Date,Temp_ave,Temp_max,Temp_min,Prec,Prec_is,SunDuration,SunDuration_is,WindSpeed,WindSpeed_max,WindDirection,WindDirection_most,RH,RH_min,Cloud 2011/1/1,6.6,10.8,2.2,0,1,7.9,0,2.3,4.6,北西,北西,33,21,5 2011/1/2,7,11,4.1,0,1,8.4,0,2.2,4.6,北北西,北北西,41,27,3.3 2011/1/3,5.9,9.2,3.1,0,1,5.2,0,2.1,4.4,北北西,北北西,48,37,7 2011/1/4,6.3,9.8,3.3,0,1,8.4,0,3,5.1,北西,北北西,40,25,3.5 2011/1/5,7.3,11,2.8,0,1,7.4,0,3.3,7.1,南西,北西,39,27,4.5
こんな感じで。
予測してみる
モジュールのimport
import pandas as pd from sklearn.linear_model import LogisticRegression from sklearn.model_selection import cross_val_score, train_test_split import numpy as np np.random.seed(1234)
データの読み込み
df_tokyo = pd.read_csv('Tokyo.csv')
df_tokyo.head()
正解ラベルを作る
降水量が1mm/day以上あった日を1、なかった日を0にします。そして、結合します。
y_label = [] for i in range(len(df_tokyo) -1): y_label.append(0 if df_tokyo.Prec[i + 1] < 1. else 1) df_label = pd.DataFrame({'label': y_label}) df = pd.concat([df_tokyo, df_label], axis=1).dropna()
データを加工する
そして、今回使わない風向関係のカラムと日付のカラムを削除します。そして切片を追加します。
df = df.drop('WindDirection', axis=1) df = df.drop('WindDirection_most', axis=1) df = df.drop('Date', axis=1) df['intercept'] = 1
データをシャッフルします。
df = df.reset_index(np.random.permutation(df.index), drop=True)
学習用データとテスト用データに分けます。
X_train, X_test, y_train, y_test = train_test_split(df.drop("label", axis=1), df.label, test_size=0.2)
学習する
Pythonでロジスティック回帰をする方法は色々あるのですが、今回は、scikit-learnに含まれる物を使ってやってみたいと思います。
lr = LogisticRegression().fit(X_train, y_train)
評価する
これで、モデルが完成しました。
モデルの評価をしましょう。
lr.score(X_test, y_test)
0.71526195899772205
なんと、72%の正解率!
偶然の可能性も否定出来ないのでcross validationしておきましょう。
cross_val_score(lr, df.drop("label", axis=1), df.label, cv=5).mean()
0.71200112335007959
71%ですね。
これが、いいのか悪いのか。2値分類だから50%超えてるのでいいのでしょうか?
正解率を評価する
なんか、70%って見覚えが。。。そう、たしか、東京の晴天率って70%くらいだったような。ということは、常に晴れといえば、この正答率は達成できるのではないか?ということで、晴れの割合を求めてみましょう。
df[df.label == 0].label.count() / df.label.count()
0.7120036513007759
これは、とりあえず「降らない」といえば当たりますね。
その他の指標も確認しましょう。
指標名 | スコア | 備考 |
見逃し率 | 0.24145785876993167 | 降らないと言ったのに降っちゃった率 0に近い方がいい |
空振り率 | 0.04328018223234624 | 降るといったのに降らなかった率 0に近い方がいい |
スレットスコア | 0.14383561643835616 | まれな現象の評価指標 1に近いほどよい |
バイアススコア | 0.31496062992125984 | 予報の偏り具合 1に近いほうが偏りがない |
スキルスコア | 0.13107849191645696 | 予報難易度を除いた精度 0でランダム、1で完全予報 |
精度を上げる
変数選択
モデルの精度を見て見るために詳細を見てみましょう。ただ、scikit-learnだと余り詳細な情報は取れないため、statsmodelsに含まれる一般化線形モデルのロジスティック回帰を使ってみようと思います。
import statsmodels.api as sm import statsmodels.formula.api as smf
mdl = smf.glm('label ~ Temp_ave + Temp_max + Temp_min + Prec + Prec_is + SunDuration + SunDuration_is + WindSpeed + WindSpeed_max + RH + RH_min + Cloud', data=df[:int(len(df)*0.8)], family = sm.families.Binomial(link=sm.families.links.logit)).fit() mdl.summary()
coef | std err | z | P>abs(z) | ||
Intercept | -3.0389 | 0.567 | -5.363 | 0.000 | |
Temp_ave | 0.3854 | 0.130 | 2.968 | 0.003 | |
Temp_max | -0.1779 | 0.076 | -2.338 | 0.019 | |
Temp_min | -0.2118 | 0.072 | -2.956 | 0.003 | |
Prec | 0.0038 | 0.004 | 0.892 | 0.373 | |
Prec_is | -0.3337 | 0.160 | -2.090 | 0.037 | |
SunDuration | 0.0477 | 0.031 | 1.527 | 0.127 | |
SunDuration_is | 4.761e-15 | 1.53e-15 | 3.106 | 0.002 | |
WindSpeed | 0.1708 | 0.119 | 1.437 | 0.151 | |
WindSpeed_max | -0.1252 | 0.066 | -1.903 | 0.057 | |
RH | -0.0102 | 0.011 | -0.923 | 0.356 | |
RH_min | 0.0220 | 0.010 | 2.105 | 0.035 | |
Cloud | 0.2661 | 0.036 | 7.338 | 0.000 |
となりました。P値の大きいものを削除してモデルを再作成してみましょう。
mdl2 = smf.glm('label ~ Temp_ave + Temp_max + Temp_min + Prec_is + SunDuration_is + RH_min + Cloud', data=df[:int(len(df)*0.8)], family = sm.families.Binomial(link=sm.families.links.logit)).fit() mdl2.summary()
coef | std err | z | P>abs(z) | |
Intercept | -3.1253 | 0.436 | -7.168 | 0.000 |
Temp_ave | 0.3922 | 0.124 | 3.163 | 0.002 |
Temp_max | -0.1847 | 0.067 | -2.749 | 0.006 |
Temp_min | -0.2082 | 0.070 | -2.982 | 0.003 |
Prec_is | -0.1905 | 0.147 | -1.295 | 0.195 |
SunDuration_is | -1.282e-15 | 3.88e-16 | -3.307 | 0.001 |
RH_min | 0.0134 | 0.005 | 2.556 | 0.011 |
Cloud | 0.2318 | 0.030 | 7.835 | 0.000 |
まだ、P値が5%水準に入らないので、削除します。これを繰り返します。そのときにAIC(赤池情報量基準)なども下がっていくのを確認します。
最終的には、
mdl3 = smf.glm('label ~ Temp_ave + Temp_max + Temp_min + SunDuration_is + RH_min + Cloud', data=df[:int(len(df)*0.8)], family = sm.families.Binomial(link=sm.families.links.logit)).fit()
で示されるものが良さそうです。というわけで、これでモデルの精度を見てみましょう。
df = df.drop('Prec', axis=1) df = df.drop('SunDuration', axis=1) df = df.drop('WindSpeed', axis=1) df = df.drop('WindSpeed_max', axis=1) df = df.drop('RH', axis=1) df = df.drop('Prec_is', axis=1) lr2 = LogisticRegression() cross_val_score(lr2, df.drop("label", axis=1), df.label, cv=5).mean()
0.70697829229985121
精度落ちてますがね。というわけで別アプローチを取る必要があるようです。
カテゴリカル・データの投入
前回、風向は、カテゴリカル変数だったので変換するのがめんどいという理由で変数に入れませんでしたが、めんどくさがらず、投入しようと思います。というわけで心機一転データの読み込みからやり直します。
df_tokyo = pd.read_csv('Tokyo.csv') y_label = [] for i in range(len(df_tokyo) -1): y_label.append(0 if df_tokyo.Prec[i + 1] < 1. else 1) df_label = pd.DataFrame({'label': y_label}) df = pd.concat([df_tokyo, df_label], axis=1).dropna() df = df.drop('WindDirection_most', axis=1) df = df.drop('Date', axis=1) df['intercept'] = 1
このあたりまでは、一緒なので、詳細は省略します。
WindDirectionに風向が入っていますが、その前に、気象庁のデータには、統計上許容できる欠損値があった場合を示す記号と、許容できない欠損があったことを示す記号が付けられてしまうことがあるので、それを取り除いて、ダミー変数化して、結合します。
wd = [] for f in df['WindDirection']: if f.endswith(']'): f = f.replace(']', '') elif f.endswith(')'): f = f.replace(')', '') wd.append(f) wd = pd.get_dummies(wd) df = pd.concat([df, wd], axis=1) df = df.drop('WindDirection', axis=1)
モデルの作成から評価
で、cross validationをします。
cross_val_score(lr, df.drop("label", axis=1), df.label, cv=5).mean()
0.7101798400266276
変わらない。。。。一応、変数選択もしてみます。
label ~ Temp_ave + Temp_max + Temp_min + RH_min + Cloud + 北北西
まで削られるようです。
では、そこまでDataFrameから削減して、
df = df.drop(['Prec', 'Prec_is', 'SunDuration', 'SunDuration_is', 'WindSpeed', 'WindSpeed_max', 'RH', '北', '北北東', '北東', '東北東', '東', '東南東', '南東', '南南東', '南', '南南西', '南西', '西南西', '西', '西北西','北西'], axis=1)
さて、cross validationを実行
cross_val_score(LogisticRegression(), df.drop("label", axis=1), df.label, cv=5).mean()
0.70469622741598281
駄目です、先生!
というわけで、更に別のアプローチを取りましょう!
次回へ続く(?)
*1:予報あり現象ありをA、予報なし現象ありをB、予報あり現象なしをC、予報なし現象なしをD、N=A+B+C+D、S=A+D、Sc={(A+B)(A+C)+(C+D)(B+D)}/Nとした時、見逃し率B/N、空振り率C/N、スレットスコアA/(A + B + C)、バイアススコア(A + C)/(A + B)、スキルスコア(S-Sc)/(N-Sc) 予報士試験に出るよ!!