F@N Ad-Tech Blog

株式会社ファンコミュニケーションズ nend・nex8のエンジニア・技術ブログ

機械学習を用いた天気予報 その2

こんにちは、気象予報士のy_kawasaki(取得してもう早、15年経つらしい)です。

前回は、精度が向上するどころか、悪化するという大惨事で終わりましたが、精度向上を目指します。

ここで、日本付近の気象学的知識を投入して考えましょう。そう、天気は西から変わるんです!というわけで、安直に、静岡くらいの天気情報を入れたいと思います。

東京と静岡を読み込んで、風向をDropします。

df_tokyo = pd.read_csv('Tokyo.csv')
df_shizuoka = pd.read_csv('Shizuoka.csv')

df_tokyo = df_tokyo.drop(['WindDirection', 'WindDirection_most', 'Date'], axis=1)
df_shizuoka = df_shizuoka.drop(['WindDirection', 'WindDirection_most', 'Date'], axis=1)

このままでは、カラム名が、東京と静岡で被ってしまうので、変更します。

df_tokyo.columns = ['T_Temp_ave', 'T_Temp_max', 'T_Temp_min', 'T_Prec', 'T_Prec_is', 'T_SunDuration', 'T_SunDuration_is', 'T_WindSpeed', 'T_WindSpeed_max', 'T_RH', 'T_RH_min', 'T_Cloud']
df_shizuoka.columns = ['S_Temp_ave', 'S_Temp_max', 'S_Temp_min', 'S_Prec', 'S_Prec_is', 'S_SunDuration', 'S_SunDuration_is', 'S_WindSpeed', 'S_WindSpeed_max', 'S_RH', 'S_RH_min', 'S_Cloud']

ラベルを作ります。

y_label = []
for i in range(len(df_tokyo) -1):
    y_label.append(0 if df_tokyo.T_Prec[i + 1] < 1. else 1)
df_label = pd.DataFrame({'label': y_label})

ラベルを結合して、

df = pd.concat([df_tokyo, df_shizuoka, df_label], axis=1).dropna()

シャッフルして、

df = df.reset_index(np.random.permutation(df.index), drop=True)

切片をつけて

df['intercept'] = 1

cross validationをします。

cross_val_score(LogisticRegression(), df.drop("label", axis=1), df.label, cv=5).mean()
0.71747329443213614

だめです。まだだめですね。ここでめげてはいけません。名古屋を追加して、風向をDropして、東京と静岡と名古屋とラベルを結合します。そしてシャッフル。切片項を付与!

df_nagoya = pd.read_csv('Nagoya.csv')

df_nagoya = df_nagoya.drop(['WindDirection', 'WindDirection_most', 'Date'], axis=1)

df = pd.concat([df_tokyo, df_shizuoka, df_nagoya, df_label], axis=1).dropna()
df = df.reset_index(np.random.permutation(df.index), drop=True)
df['Intercept'] = 1

で、恒例のcross validation、

cross_val_score(LogisticRegression(), df.drop("label", axis=1), df.label, cv=5).mean()
0.74850792065819993

精度が改善しました!念のため大阪も入れてみましょう。

df_osaka = pd.read_csv('Nagoya.csv')

df_osaka = df_osaka.drop(['WindDirection', 'WindDirection_most', 'Date'], axis=1)

df = pd.concat([df_tokyo, df_shizuoka, df_nagoya, df_osaka, df_label], axis=1).dropna()
df = df.reset_index(np.random.permutation(df.index), drop=True)
df['Intercept'] = 1
cross_val_score(LogisticRegression(), df.drop("label", axis=1), df.label, cv=5).mean()
0.74257704829365201

特に大きな変化は見られません。
流石に変数が多すぎて、良いモデルとはいえないので、変数選択をして、変数を減らしてみました。

df_t1 = df_tokyo['T_Prec_is']
df_s1 = pd.concat([df_shizuoka['S_SunDuration'], df_shizuoka['S_WindSpeed'], df_shizuoka['S_WindSpeed_max'], df_shizuoka['S_Cloud']], axis=1)
df_n1 = pd.concat([df_nagoya['N_Prec'], df_nagoya['N_Prec_is'], df_nagoya['N_WindSpeed_max'], df_nagoya['N_RH']], axis=1)
df_o1 = pd.concat([df_osaka['O_RH'], df_osaka['O_Cloud']], axis=1)

df = pd.concat([df_t1, df_s1, df_n1, df_o1, df_label], axis=1).dropna()
df = df.reset_index(np.random.permutation(df.index), drop=True)
df['Intercept'] = 1
cross_val_score(LogisticRegression(), df.drop("label", axis=1), df.label, cv=5).mean()
0.74668767747371056

他の指標も確認しておきます。

指標名 スコア
見逃し率 0.15034168564920272
空振り率 0.07061503416856492
スレットスコア 0.3618421052631579
バイアススコア 0.7107438016528925
スキルスコア 0.3922010819143318

全体的に改善はしています。

このアプローチでは、この程度が限界のようです。なかなか難しいですね。

機械学習を用いた天気予報 その1

こんにちは、気象予報士の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で完全予報

*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) 予報士試験に出るよ!!

ベイズ統計モデリングに触れてみた話

こんにちは、情報科学技術研究所のy_kawasakiです。春になり天気予報が微妙にずれる(時間的に)のが気になる今日このごろです。

  • ベイズ統計モデリングとは
  • 準備
    • ベイズ統計モデリングソフト
    • PyStanのインストール
  • 実例
    • 基礎統計的データ
    • 結果
    • 結果2
    • 今後の課題
  • 参考図書
  • 付録

ベイズ統計モデリングとは

ベイズ統計学による、統計モデリングのことです!それぞれの意味は個々でググってください!

準備

ベイズ統計モデリングソフト

ベイズ統計モデリングをやるには、MCMC(マルコフ連鎖モンテカルロ法)で解く必要があるので様々なパッケージが用意されています。たとえば、BUGS言語を使ったソフトや、Stanというツールなど。ここでは、最近流行りのStanを使うことにします。Pythonから使う必要が(自分には)あるため、PyStanを利用します。

PyStanのインストール


> pip install PyStan

難しいことはありません。

実例

弊社では本社があるビルが手狭になったことから一部の部署が切り出されて、近くのビルに放り出されました(以下、本社とサテライト)。そこで、時々、本社と行き来する必要が発生して、片道どれくらいかかるのかを情報科学技術研究所的には調査する必要性が発生しました。(使命感)

そこで、被験者3名に計測を依頼したところ、A7件、B2件、C4件、計13件の下記の結果を得ることができました。なお、本社、サテライト間には、信号とエレベータが障害として存在しています。

続きを読む

『ユーザの平均継続期間が「1/解約率」』であるための十分条件について

 こんにちは、サービス開発部情報科学技術研究所所属のk_oomoriです。先日、
migi.hatenablog.com
というブログが公開され、アドテク界隈で話題になったようです。私も読ませていただきましたが、平易に書かれていてとても良い記事だと思いました。僭越ながら要約させていただくと、

  1. ある段階でのユーザ数をUとし、その後の流入は考えない
  2. 毎月の解約率(常に一定で変わらないとする!)はC

という仮定をおくとnか月後のユーザ数はU(1-C)^nとなり、それを0から\inftyまで足しあげた延べ継続期間が
\displaystyle
\sum_{n=0}^{\infty}U(1-C)^n=\frac{U}{C}
となるため、平均継続期間=ユーザの延べ継続期間/ユーザ数 として定義した平均継続時間が1/解約率で与えられることが示された、というものです。
 似たような問題は実は他の文脈でも現れます。例えば原子核物理学の分野で不安定な原子核の放射性崩壊という現象があり、ある時刻t=0における不安定原子核の数をN_0崩壊率あるいは崩壊定数(崩壊のしやすさ、詳しくは後ほど説明します)を\lambdaとしたときに、時刻tにおいて崩壊せずに残っている原子核数N(t)N_0\lambdaを用いて表すことができます。この場合、1個の原子核が平均的にどのくらい生き延びるか?という問題を考えることができ、平均寿命と呼ばれています。これを先ほどのユーザの解約問題と比較すると、

  • 解約率 ⇔ 崩壊率 (どちらも減る割合を表す。以後まとめて減衰率と表現することにする)
  • 平均継続期間 ⇔ 平均寿命 (解約/崩壊せずに生き残る時間間隔)

によって対応付けできるため、数学的には同種の問題とみなすことができます。(離散、連続の違いはありますが)

続きを読む