このエントリーは以前書いた機械学習で積雪の有無を予測してみたの続編です。この時は積雪の有無(1か0か)だけを予測したのですが、もうちょっと頑張って積雪量の変化を予測してみました。
先に結果を記しとくと、こんな感じになりました。横軸が日数、縦軸が積雪量(cm)です。
結果その1（青が実際の積雪量、赤線が予測した積雪量）
結果その２（青が実際の積雪量、赤線が予測した積雪量）
「結果その１」と「結果その２」がそれぞれ何なのかは以下を読んでみてください。
やりたかったこと
以前、機械学習で積雪の有無を予測してみた で scikit-learn を使って積雪の有無を予測してみたのですが、ちょっと欲が出てきて、有無じゃなくてまとまった期間の実際の積雪量（cm）を予測してみたい、と思ってやってみました。
具体的には、気象庁から提供される 積雪量 風速 温度 などの気象データを取得して、そのうち最初の約7500日分のデータを学習に使用し、残り2年間（365x2=730日）の積雪量の変化を予測し、実際の積雪量の変化と比較します。
学習用データを集める
学習用データは気象庁が公開しているものを取得します。具体的には以前書いた機械学習で積雪の有無を予測してみたを参照してみてください。
得られるCSVのデータはこんな感じ。ターゲットにしたのは雪の多い富山の砺波市というところです。
data_2013_2015.csv
ダウンロードした時刻：2016/03/20 20:31:19
,砺波,砺波,砺波,砺波,砺波,砺波,砺波,砺波,砺波,砺波,砺波,砺波,砺波,砺波
年月日時,気温(℃),気温(℃),気温(℃),積雪(cm),積雪(cm),積雪(cm),風速(m/s),風速(m/s),風速(m/s),風速(m/s),風速(m/s),降水量(mm),降水量(mm),降水量(mm)
,,,,,,,,,風向,風向,,,,
,,品質情報,均質番号,,品質情報,均質番号,,品質情報,,品質情報,均質番号,,品質情報,均質番号
2013/2/1 1:00:00,-3.3,8,1,3,8,1,0.4,8,西,8,1,0.0,8,1
2013/2/1 2:00:00,-3.7,8,1,3,8,1,0.3,8,北,8,1,0.0,8,1
2013/2/1 3:00:00,-4.0,8,1,3,8,1,0.2,8,静穏,8,1,0.0,8,1
2013/2/1 4:00:00,-4.8,8,1,3,8,1,0.9,8,南南東,8,1,0.0,8,1
...
基本的な考え方
考え方としては、多分こういう予測では標準的なものではないかと思いますが、何種類かの周辺データと結果としての積雪量をセットにしてモデルに学習させ、結果得られたモデルに周辺データのみを与え、積雪量の予測値を得る、というものです。いわゆる「教師あり学習」ってやつです。
今回の場合は、以下のデータを周辺データとして使用しました。
- 温度
- 風速
- 昨日の積雪量
- 1日前の温度,2日前の温度,3日前の温度
- 1日前の風速, 2日前の風速, 3日前の風速
イメージで表すとこんな感じです。
[温度,風速,昨日の積雪量,1日前の温度,2日前の温度,3日前の温度, 1日前の風速, 2日前の風速, 3日前の風速] → 当日の積雪量
[温度,風速,昨日の積雪量,1日前の温度,2日前の温度,3日前の温度, 1日前の風速, 2日前の風速, 3日前の風速] → 当日の積雪量
[温度,風速,昨日の積雪量,1日前の温度,2日前の温度,3日前の温度, 1日前の風速, 2日前の風速, 3日前の風速] → 当日の積雪量
....
[温度,風速,昨日の積雪量,1日前の温度,2日前の温度,3日前の温度, 1日前の風速, 2日前の風速, 3日前の風速] → 当日の積雪量
で、これをもとに、周辺データのみを与えて予測値を得る
[温度,風速,昨日の積雪量,1日前の温度,2日前の温度,3日前の温度, 1日前の風速, 2日前の風速, 3日前の風速] → （当日の積雪量予測値）
こんな流れで行いました。なお、基本的には与えるのは予測対象日のデータですが、一つだけ 昨日の積雪量 だけは予測対象日の１日前のデータになります。そしてこれが与えるデータの中で一番影響が大きいようでした。ま、考えてみれば当たり前ですけども。
最初に書いた通り、気象庁から得たデータのうち約7500日分のデータを学習に使用し、残り2年間の積雪量の変化を予測し、実際の積雪量の変化と比較します。
予測してみる
実際のコードは以下のようになります。
snow_forecaster.py
import csv
import numpy as np
from matplotlib import pyplot
from sklearn import linear_model
from sklearn import cross_validation
class SnowForecast:
def __init__(self):
u"""各インスタンス変数を初期化"""
self.model = None # 生成した学習モデル
self.data = [] # トレーニング用データの配列
self.target = [] # 実際の積雪量の配列
self.predicts = [] # 積雪量の予測値の配列
self.reals = [] # 実際の積雪量の配列
self.day_counts = [] # 開始日からの経過日の配列
self.date_list = []
self.record_count = 0
def load_csv(self):
u"""学習用CSVファイルを読み込む"""
with open("sample_data/data.csv", "r") as f:
reader = csv.reader(f)
accumulation_yesterday0 = 0
date_yesterday = ""
temp_3days = []
wind_speed_3days = []
for row in reader:
if row[4] == "":
continue
daytime = row[0] # "yyyy/mmdd HH:MM:SS"
date = daytime.split(" ")[0] # "yyyy/mm/dd"
temp = int(float(row[1])) # 温度。微妙に影響あり
wind_speed = float(row[7]) # 風速。微妙に影響あり
precipitation = float(row[12]) # 降水量。影響なし
accumulation = int(row[4]) # 積雪量。昨日の降雪量は大いに影響あり
if len(wind_speed_3days) == 3:
# トレーニング用データ
# [温度,風速,昨日の積雪量,1日前の温度,2日前の温度,3日前の温度, 1日前の風速, 2日前の風速, 3日前の風速]
sample = [temp, wind_speed, accumulation_yesterday0]
sample.extend(temp_3days)
sample.extend(wind_speed_3days)
self.data.append(sample)
self.target.append(accumulation)
if date_yesterday != date:
accumulation_yesterday0 = accumulation
self.date_list.append(date)
wind_speed_3days.insert(0, wind_speed)
if len(wind_speed_3days) > 3:
wind_speed_3days.pop()
temp_3days.insert(0, temp)
if len(temp_3days) > 3:
temp_3days.pop()
date_yesterday = date
self.record_count = len(self.data)
return self.data
def train(self):
u"""学習モデルを生成する。学習用データは元データの約7500日目までを使う"""
x = self.data
y = self.target
print(len(x))
# ElasticNetCV,LassoCV,RidgeCVからもっとも誤差の小さかったElasticNetCVを選択
model = linear_model.ElasticNetCV(fit_intercept=True)
model.fit(x[0:self.training_data_count()], y[0:self.training_data_count()])
self.model = model
def predict(self):
u"""学習モデルを使って予測する。予測は最後の2年分に対して行う"""
x = self.data
y = self.target
model = self.model
for i, xi in enumerate(x):
real_val = y[i]
if i < self.training_data_count() + 1:
self.predicts.append(0)
self.reals.append(real_val)
self.day_counts.append(i)
continue
predict_val = int(model.predict([xi])[0])
# 積雪量予測が0以下の場合は0とする。
if predict_val < 0:
predict_val = 0
self.predicts.append(predict_val)
self.reals.append(real_val)
self.day_counts.append(i)
def show_graph(self):
u"""グラフで予測値と実測値を比較する"""
pyplot.plot(self.day_counts[self.predict_start_num():], self.reals[self.predict_start_num():], "b")
pyplot.plot(self.day_counts[self.predict_start_num():], self.predicts[self.predict_start_num():], "r")
pyplot.show()
def check(self):
u"""トレーニング用データと予測データとの誤差を測定"""
x = np.array(self.data[self.predict_start_num():])
y = np.array(self.target[self.predict_start_num():])
model = self.model
p = np.array(self.predicts[self.predict_start_num():])
e = p - np.array(self.reals[self.predict_start_num():])
error = np.sum(e * e)
rmse_10cv = np.sqrt(error / len(self.data[self.predict_start_num():]))
print("RMSE(10-fold CV: {})".format(rmse_10cv))
def training_data_count(self):
u"""最後の2年分は残して、それ以前を学習用データとして使用する。学習用データ数を返す"""
return self.record_count - 365 * 2
def predict_start_num(self):
u"""最後の2年分は予測させ、実測値との誤差を測定するために使用する。予測開始位置を返す"""
return self.training_data_count() + 1
if __name__ == "__main__":
forecaster = SnowForecast()
forecaster.load_csv()
forecaster.train()
forecaster.predict()
forecaster.check()
forecaster.show_graph()
一番めんどくさいのは、生データから前章のような学習データを作る部分でした。それでもpythonだから楽なもんですけど。
で、実行結果はこの通り（青が実際の積雪量、赤線が予測した積雪量）。最初に示した「結果その１」がこれ。
なんかそれっぽい感じで予測しています。
ここで、ふとこのやり方に疑問が浮かびました。
「でも１日前の積雪量を与えて予測してるんだから、これだと実際に未来予測に使おうとした時、予測できるのは明日の積雪量のみじゃね...?」
いや、わかってますよ？それを言うなら温度や風速だってそうだろうって。ただほら、それらは天気予報とかで...ゲフンゲフン
自身が予測した昨日分積雪量を使って翌日分積雪量を予測するように変更
で、早速そのようにコードを修正してみました。
モデルの学習部分については、特に変更点はありません。積雪量予測を行う時に与えるデータのうち 昨日分積雪量 を実測値ではなく 自身が予測した１日前の予測値 で入れ替えて予測させてみます。
コードは以下の通り。変更したのは predict 関数のみです。
snow_forecaster.py
import csv
import numpy as np
from matplotlib import pyplot
from sklearn import linear_model
from sklearn import cross_validation
class SnowForecast:
def __init__(self):
u"""各インスタンス変数を初期化"""
self.model = None # 生成した学習モデル
self.data = [] # トレーニング用データの配列
self.target = [] # 実際の積雪量の配列
self.predicts = [] # 積雪量の予測値の配列
self.reals = [] # 実際の積雪量の配列
self.day_counts = [] # 開始日からの経過日の配列
self.date_list = []
self.record_count = 0
def load_csv(self):
u"""学習用CSVファイルを読み込む"""
with open("sample_data/data.csv", "r") as f:
reader = csv.reader(f)
accumulation_yesterday0 = 0
date_yesterday = ""
temp_3days = []
wind_speed_3days = []
for row in reader:
if row[4] == "":
continue
daytime = row[0] # "yyyy/mmdd HH:MM:SS"
date = daytime.split(" ")[0] # "yyyy/mm/dd"
temp = int(float(row[1])) # 温度。微妙に影響あり
wind_speed = float(row[7]) # 風速。微妙に影響あり
precipitation = float(row[12]) # 降水量。影響なし
accumulation = int(row[4]) # 積雪量。昨日の降雪量は大いに影響あり
if len(wind_speed_3days) == 3:
# トレーニング用データ
# [温度,風速,昨日の積雪量,1日前の温度,2日前の温度,3日前の温度, 1日前の風速, 2日前の風速, 3日前の風速]
sample = [temp, wind_speed, accumulation_yesterday0]
sample.extend(temp_3days)
sample.extend(wind_speed_3days)
self.data.append(sample)
self.target.append(accumulation)
if date_yesterday != date:
accumulation_yesterday0 = accumulation
self.date_list.append(date)
wind_speed_3days.insert(0, wind_speed)
if len(wind_speed_3days) > 3:
wind_speed_3days.pop()
temp_3days.insert(0, temp)
if len(temp_3days) > 3:
temp_3days.pop()
date_yesterday = date
self.record_count = len(self.data)
return self.data
def train(self):
u"""学習モデルを生成する。学習用データは元データの約7500日目までを使う"""
x = self.data
y = self.target
print(len(x))
# ElasticNetCV,LassoCV,RidgeCVからもっとも誤差の小さかったElasticNetCVを選択
model = linear_model.ElasticNetCV(fit_intercept=True)
model.fit(x[0:self.training_data_count()], y[0:self.training_data_count()])
self.model = model
def predict(self):
u"""学習モデルを使って積雪量を予測する。予測は最後の２年間分に対して行う"""
x = self.data
y = self.target
model = self.model
yesterday_predict_val = None # 昨日分の予測値を入れておく変数
for i, xi in enumerate(x):
real_val = y[i]
if i < self.training_data_count() + 1:
self.predicts.append(0)
self.reals.append(real_val)
self.day_counts.append(i)
continue
# 昨日分積雪量を自身が予測した昨日分の予測値と入れ替える
if yesterday_predict_val != None:
xi[2] = yesterday_predict_val
predict_val = int(model.predict([xi])[0])
# 積雪量予測が0以下の場合は0とする。
if predict_val < 0:
predict_val = 0
self.predicts.append(predict_val)
self.reals.append(real_val)
self.day_counts.append(i)
yesterday_predict_val = predict_val
def show_graph(self):
u"""グラフで予測値と実測値を比較する"""
pyplot.plot(self.day_counts[self.predict_start_num():], self.reals[self.predict_start_num():], "b")
pyplot.plot(self.day_counts[self.predict_start_num():], self.predicts[self.predict_start_num():], "r")
pyplot.show()
def check(self):
u"""トレーニング用データと予測データとの誤差を測定"""
x = np.array(self.data[self.predict_start_num():])
y = np.array(self.target[self.predict_start_num():])
model = self.model
p = np.array(self.predicts[self.predict_start_num():])
e = p - np.array(self.reals[self.predict_start_num():])
error = np.sum(e * e)
rmse_10cv = np.sqrt(error / len(self.data[self.predict_start_num():]))
print("RMSE(10-fold CV: {})".format(rmse_10cv))
def training_data_count(self):
u"""最後の2年分は残して、それ以前を学習用データとして使用する。学習用データ数を返す"""
return self.record_count - 365 * 2
def predict_start_num(self):
u"""最後の2年分は予測させ、実測値との誤差を測定するために使用する。予測開始位置を返す"""
return self.training_data_count() + 1
if __name__ == "__main__":
forecaster = SnowForecast()
forecaster.load_csv()
forecaster.train()
forecaster.predict()
forecaster.check()
forecaster.show_graph()
結果はこうなりました（青が実際の積雪量、赤線が予測した積雪量）。最初に示した「結果その２」。
うーん。さすがに実際の昨日分積雪量を与えた時より不正確になりました。ただ、波形としてはそこまでめちゃくちゃでもない感じかな？
感想など
もっとめちゃくちゃな予測値になるかなーと思っていたのですが、意外にそれっぽく予測できたなって思いました。
ただ、途中ゲフンゲフンでうまく誤魔化したけど、予測するときに与える気温や風速は当日の実測値を使ってるんですよね。でも将来の一定期間の予測をしようと思ったら、当然それらも別途予測された値を使うか、そもそもそういう値を使うのをやめるかしかないわけで、予測された値を使うとしたらもっと精度は下がるでしょうね。しかも未来であればあるほど。なので、こういうことをやりたかったら、予測値を使って予測して、さらにそれを使って予測して、みたいになって、後ろの方ほど、前行程でのちょっとした誤差が大きく拡大するなーって思いました。というわけで気象庁さん頑張って（