こちらの記事は筆者がM1の時にデータサイエンスの授業のレポートを記事化したものです。
今回のエントリーでは実際にkaggleをやるわけではないのですが、我々でも簡単に手に入るデータからkaggleっぽい遊びをやってみようというテーマでお送り致します。
気象庁データから翌日の天気を予測する
このテーマでやってみようと思います。
データは気象庁のホームページから入手しました。
こちらのサイトでは、地点・項目・期間を任意に設定して気象に関するデータを獲得できます。
今回は以下の設定でデータを入手しました。
| 地点 | 項目 | 期間 |
|---|---|---|
| 東京 | 気温(平均・最高・最低), 降水量, 日照時間, 日射量, 風速(平均・最大・最大瞬間) , 蒸気圧, 相対湿度, 気圧(現地・海面), 雲量, 天気概況(昼) |
2007~2017年における 9/31~10/31 |
予測する項目(目的変数)は天気概況(昼)で特徴量(説明変数)がそれ以外の項目にあたります。
以上のデータは、先ほどのホームページからcsv形式でダウンロードできます。
ここからこのcsvファイルを読み取り処理していくのですが、ダウンロードしたての生の状態では、ダウンロード時刻や品質番号、均質番号といった余分な項目が含まれているため、こういったものをあらかじめexcelを使って取り除きました。
また翌日の天気を予測するタスクのため、日付と天気概況は一日分進めます。成形後は以下のイメージ。
| 年月日 | 特徴量項目(気温etc...) | 天気概況(昼：06時〜18時) |
|---|---|---|
| 2007/10/1 | hogehoge | 曇時々雨 |
それとcsvファイルは文字コードがshift-jisなのでutf-8に変更することもお忘れなく。(nkfとかツールで)
csvファイルを読み込もう
pandasを使ってcsvファイルを読み込みます。今回はindexとして日付型のデータを用いているのでindex及びparse_datesに'年月日'を指定しましょう。これで、indexがdatetime型で扱えます。
import pandas as pd
df = pd.read_csv('data_tokyo.csv', index_col='年月日', parse_dates=['年月日'])
目的変数の成形
今回のデータでは目的変数である天気概況が詳細に記述されているので、おおざっぱに「晴・曇・雨」に分類します。
データをよく眺めてみると、大雨や晴時々曇のように記述にはパターンがあり、接続詞的なものに「'後','一時','時々',' 、'」が使われていることがわかります。
また修飾語がついている場合は、各トークンの最後の文字に注目するとそれが「晴・曇・雨」のいずれかの文字が使われています。
よって正規表現や文字列のインデックスの指定で「晴・曇・雨」に分類できそうです。
import re
import datetime
#以降、目的変数はcolで指定します。
col = '天気概況(昼：06時〜18時)'
df[col] = df[col].apply(lambda x:re.split('[後|一時|時々|、]',x)[0][-1])
#分類できたか確認
df[col].value_counts()
晴 143
曇 141
雨 57
Name: 天気概況(昼：06時〜18時), dtype: int64
無事、「晴・曇・雨」に分類できました!
ついでに目的変数のラベルを整数に変換します。dataframeの型も忘れずに(これのせいで学習時にエラーが発生し問題解決に1時間とられました...)
#ラベルを整数に変換
df.loc[df[col]=='晴',col] = 0
df.loc[df[col]=='曇',col] = 1
df.loc[df[col]=='雨',col] = 2
#目的変数の型変換
df.loc[:,[col]]=df.loc[:,[col]].astype(int)
データの分割
さて次に訓練データとテストデータを分割します。今回は、天気を予測するということで、2017年のデータと2007~2016年のデータに分割します。
訓練データを2016年以前に、テストデータを2017年のものとします。(indexが日付型なのでこの作業は簡単)
traindf = df.loc[:'2016',:]
testdf = df.loc['2017':]
またそれぞれ説明変数と目的変数でも分割します。
X_train = traindf.loc[:'2016',:].drop(col,axis=1).values
y_train = traindf.loc[:'2016',col].values
X_test = testdf.loc['2017':,:].drop(col,axis=1).values
y_test = testdf.loc['2017':,col].values
特徴量の標準化
ランダムフォレストを使う場合だと、ほぼやる必要は無いですが、境界面を作っていくタイプのモデル(重回帰分析,ロジスティック回帰,SVMとか?)だとデータを標準化する方が予測精度が上がりやすい知見があるのでやります。(標準化についてはググればブログやwikiで解説があります)
scikit-learnではこの標準化が以下のコードで簡単に行えます.
ただしテストデータは未知データであるという前提があるため、標準化における平均と分散の計算には用いてはいけません。(重要)
from sklearn.preprocessing import StandardScaler
scaler = StandardScaler()
#標準化させるために訓練データだけで平均と分散を計算
scaler.fit(X_train)
#訓練データ、テストデータを標準化する。
X_train = scaler.transform(X_train)
X_test = scaler.transform(X_test)
とりあえず学習
では一旦、ここまでで獲得したデータからモデルを作成してみます。
モデルはscikit-learnに標準装備されてるサポートベクターマシン(SVM)で行いました。
from sklearn.svm import SVC
clf = SVC()
clf.fit(X_train, y_train)
y_pred = clf.predict(X_test)
とりあえず評価
まぁ評価してみましょう。各ラベルの予測精度の平均はモデルのメソッドであるscoreで確認できます。
またそれぞれの予測精度はscikit-learnのclassification_reportで確認できます。
平均予測精度 : 0.3225806451612903
precision recall f1-score support
晴 0.24 0.75 0.36 8
曇 0.50 0.30 0.37 10
雨 0.00 0.00 0.00 13
avg / total 0.22 0.29 0.21 31
めちゃくちゃ悪いですね。 特に雨は全く予測できてません、、、
当てずっぽで3分の1なんでそれより低いのは使い物になりません。
作戦変更...とりあえず晴れだけでも予測したい
上の通りです。というわけでラベルを以下のものに書き換えました。
#ラベルを整数に変換
df.loc[df[col]=='晴',col] = 1
df.loc[df[col]=='曇',col] = 0
df.loc[df[col]=='雨',col] = 0
平均予測精度 : 0.7096774193548387
precision recall f1-score support
曇or雨 0.76 0.83 0.79 23
晴 0.33 0.25 0.29 8
avg / total 0.65 0.68 0.66 31
少しまともに見えますね。
パラメータチュニーングする
パラメータについて
モデルを構築する際のキモになるところです。
今回利用する学習モデルのSVMには、境界面の複雑さを制御するパラメータとして$C$と$\gamma$があります。
私のイメージでは、$C$と$\gamma$が高いほどoverfittingする感じです。(こちらのサイトでわかりやすく解説)
またそれぞれのクラスのサンプル数に偏りがあるデータ(不均衡なデータ)の場合、片方の分類精度が極端に悪くなってしまう場合があります。(参考:unbalanced-problems)
そこで各クラスでバランス良く重み付けしてくれるclass_weight='balanced'を設定するか否かという点についてもチューニングして判断します。
時系列データにおける交差検証の注意点
交差検証についての説明は割愛しますが、今回の場合は一般的なランダム抽出の交差検証ではなく、各年度ごとにグループを作成して行います。
このアプローチをとった理由は、時系列データは自己相似性があり、indexの近いデータ同士は類似しやすい特性があるためです。
したがって、ある年度のデータをその他の年度のデータを用いて予測するということを繰り返し、パラメータのチュニーングを行なっていく必要があります。
scikit-learnではこのための交差検証のクラスにGroupKFoldがあります。これによって作成した分割のリストをGridSearchCVの引数cvに投げることで前述した交差検証が行えます。
標準化と学習のパイプライン化
標準化の項目でも説明したようにテストデータは平均と分散の計算に用いていけないという掟があります。
交差検証も同様で検証データ(validation data)は別に平均と分散を計算しなければなりません。
これをscikit-learnのグリッドサーチで実現するために、標準化と学習をパイプライン化します。
この関数もscikit-learnに標準装備されており、そのまんまPipelineクラスで行えます。
scikit-learnは本当に便利ですよね〜。
# group 10-fold cross validation
from sklearn.pipeline import Pipeline
from sklearn.preprocessing import StandardScaler
from sklearn.svm import SVC
from sklearn.grid_search import GridSearchCV
from sklearn.model_selection import GroupKFold
#標準化と学習のパイプライン
pipe = Pipeline([('scaler',StandardScaler()),("svm",SVC(random_state=0))])
pipe.fit(X_train,y_train)
#パラメータは10^nオーダーで変化させる
params = {'svm__C':[0.0001,0.001,0.01,0.1,1,10,100],'svm__gamma':[0.001,0.01,0.1,1,10,100],'svm__class_weight':[None,'balanced']}
#年度は31インデックスごとに切り替わるため
groups = sum([[label]*31 for label in range(10)],[])
gkfold = list(GroupKFold(n_splits=10).split(X_train,y_train,groups))
#グリッドサーチ
grid = GridSearchCV(pipe, param_grid = params, cv = gkfold, scoring='accuracy',n_jobs=-1)
grid.fit(X_train,y_train)
print('Best cross-validation accuracy: {:.2f}'.format(grid.best_score_))
print('Train set score: {:.2f}'.format(grid.score(X_train,y_train)))
print('Best parameters : {}'.format(grid.best_params_))
Best cross-validation accuracy: 0.65
Train set score: 0.67
Best parameters : {'svm__gamma': 0.01, 'svm__class_weight': 'balanced', 'svm__C': 1}
それでは、チューニングしたパラメータでテストデータを確認しましょう。
scaler = StandardScaler()
scaler.fit(X_train)
X_train = scaler.transform(X_train)
X_test = scaler.transform(X_test)
from sklearn.svm import SVC
clf = SVC(C=1.0, gamma=0.01, random_state=0, probability=True, class_weight='balanced')
clf.fit(X_train, y_train)
print('平均予測精度 : {}\n'.format(clf.score(X_test, y_test)))
print(classification_report(y_test, pred,target_names=['曇or雨','晴']))
平均予測精度 : 0.6129032258064516
precision recall f1-score support
曇or雨 0.76 0.83 0.79 23
晴 0.33 0.25 0.29 8
avg / total 0.65 0.68 0.66 31
チューニング前より下がっていました。
まとめ
というわけで、まさかのチューニングしない方がよかったという結果になってしまいました。
色々理由は考えられるのですが、今年のデータはこれまでに無い異常気象だったということでなんでしょうかね。
今回のエントリーはここでおしまいです。ありがとうございました！