最近Rに入門したので、練習課題として、SVM(kernlab)を使って過去の気象情報から明日の天気を予測するということをやってみた。
なお、私は気象学については全く素養が無いので、おそらく専門の人が見たら直ちにグーで殴ってくるようなことをやっているかもしれないが、気象予測というのは今回の本質でないので、大目に見てください。
(とはいえ、あまりにも頓珍漢なことを言っていたらご指摘ください。)
気象データの取得
気象庁の気象データダウンロードサービスから、滋賀県彦根市の2000年1月1日から2013年12月31日までの日別の気象情報を取得。
取得項目は私の勘により、天気に関係のありそうな下記の7項目とした。
- 日最高気温
- 日最低気温
- 日降水量
- 日平均風速
- 日平均気圧
- 日平均相対湿度
- 天気概況(昼)
注意点として、上記のサービスはダウンロード対象のデータサイズがある程度より大きくなるとInternal Server Errorを返すことがあるので、様子を見て期間を区切り、少量ずつダウンロードしなければならない。
(私は14年分を3回に分けてダウンロードした。)
データの整形
ダウンロードした時刻：2014/08/31 13:32:34
,彦根,彦根,彦根,彦根,彦根,彦根,彦根,彦根,彦根,彦根,彦根,彦根,彦根,彦根,彦根,彦根,彦根,彦根,彦根,彦根,彦根,彦根
年月日,最高気温(℃),最高気温(℃),最高気温(℃),最低気温(℃),最低気温(℃),最低気温(℃),降水量の合計(mm),降水量の合計(mm),降水量の合計(mm),降水量の合計(mm),平均風速(m/s),平均風速(m/s),平均風速(m/s),平均現地気圧(hPa),平均現地気圧(hPa),平均現地気圧(hPa),平均湿度(％),平均湿度(％),平均湿度(％),天気概況(昼),天気概況(昼),天気概況(昼)
,,品質情報,均質番号,,品質情報,均質番号,,現象なし情報,品質情報,均質番号,,品質情報,均質番号,,品質情報,均質番号,,品質情報,均質番号,,品質情報,均質番号
2000/1/1,8.5,8,1,3.3,8,1,0.0,0,8,1,4.5,8,1,1013.2,8,1,67,8,1,曇,8,1
2000/1/2,9.1,8,1,-0.1,8,1,0.0,0,8,1,1.4,8,1,1007.8,8,1,73,8,1,晴後曇,8,1
2000/1/3,12.0,8,1,5.0,8,1,0.0,0,8,1,3.4,8,1,1004.8,8,1,75,8,1,晴時々曇一時雨,8,1
...
気象庁から取得したCSVファイルは、5行のヘッダを持ち、指定した項目の他に、それらについての品質情報などを含むので、必要な項目のみ抽出する。
また、天気概況が日本語である上に情報が細かすぎて扱いづらいので、「のち曇り」のような情報を除去し、適当にコード値っぽいものに置き換える。(真面目にやる場合は、こんな乱暴な扱いをしてはいけない。)
#date t_max t_min rain wind atom. hum. weather
2000/1/1 8.5 3.3 0.0 4.5 1013.2 67 Kumori
2000/1/2 9.1 -0.1 0.0 1.4 1007.8 73 Hare
2000/1/3 12.0 5.0 0.0 3.4 1004.8 75 Hare
2000/1/4 9.3 3.9 0.0 4.1 1011.4 64 Kumori
2000/1/5 9.7 1.3 0 1.5 1013.0 71 UsuKumori
...
素性ベクトルの構成
上で整形したデータを使って、SVMに喰わせる素性ベクトルを構成する。
- 予測したいのは明日の天気なので、「今日の気象情報」と「明日の天気」が同じ行に入るようにする。
- 日付については、「月」と「日」のみ素性に含める。
- 14年分の気象データうち、古い10年分を学習用に、新しい4年分を評価用に使う。
- 離散値である天気概況については学習用データか評価用データのどちらか片方にしか含まれない値が存在すると評価ができなくなるので、クラスの構成を調整するなどしてすべての値が両方のデータ群に含まれるようにする。
結果、素性ベクトルの要素は
- 月 (1..12)
- 月内の日 (1..31)
- 今日の最高気温
- 今日の最低気温
- 今日の降水量
- 今日の平均風速
- 今日の平均気圧
- 今日の平均相対湿度
- 今日の天気概況(昼)
となり、これに予測対象の
- 明日の天気概況(昼)
を加えたものが教師データとなる。
天気概況は、「薄曇り」や「快晴」など細かいものを統合して、「雨」「晴れ」「曇り」「雪」の4種類とした。
1 1 8.5 3.3 0.0 4.5 1013.2 67 Kumori Hare
1 2 9.1 -0.1 0.0 1.4 1007.8 73 Hare Hare
1 3 12.0 5.0 0.0 3.4 1004.8 75 Hare Kumori
1 4 9.3 3.9 0.0 4.1 1011.4 64 Kumori Kumori
1 5 9.7 1.3 0 1.5 1013.0 71 Kumori Kumori
...
Rの準備
今回はkernlabパッケージを使うのでインストールする。
$ sudo R
> install.packages("kernlab")
> library(kernlab)
学習
通常、SVMに学習データを与える際には、要素ごとの絶対値の大小によって偏りが発生しないようにベクトルを正規化したり、Factor型のような離散値に適当に番号を振ったりする必要があるが、ksvmはその辺を勝手にやってくれるのでお気楽に使える。
カーネル関数は、データの傾向がよく分からないときに使われる(らしい)RBFカーネルを使う。
wet_train <- read.delim("train.tsv", sep="\t", stringsAsFactors=T, header=F)
# 学習を行う。V10というのが、予測対象とするカラム名。
wet_model <- ksvm(V10 ~ ., data=wet_learn, type = "C-bsvc", kernel = "rbfdot")
評価
wet_test <- read.delim("test.tsv", sep="\t", stringsAsFactors=T, header=F)
# predictを使って予測値を求める。
result <- predict(wet_model, wet_test)
# 予測結果を度数分布表で表示
table(result, wet_test$V10)
result Ame Hare Kumori Yuki
Ame 125 3 35 6
Hare 2 425 198 0
Kumori 67 106 449 14
Yuki 1 0 5 24
| result | Ame | Hare | Kumori | Yuki | 正答率(適合率) |
|---|---|---|---|---|---|
| Ame | 125 | 3 | 35 | 6 | 0.74 |
| Hare | 2 | 425 | 198 | 0 | 0.68 |
| Kumori | 67 | 106 | 449 | 14 | 0.71 |
| Yuki | 1 | 0 | 5 | 24 | 0.80 |
| 正答率(再現率) | 0.64 | 0.80 | 0.65 | 0.55 | 0.70 |
縦が予測結果で、横が実際の天気。
たとえばHare行Kumori列の198というのは、晴れと予報して実際には曇りだった回数が198回であるということを意味する。
考察
結果の度数分布表から、予測結果が正解である確率はおよそ0.7で、一見そこそこの結果であるように見えるが、再現率の方を見ると、雪が降る事象を2回に1回しか予言できていないなど、天気予報としては全く実用的とは言えない。
予測期間を長くすれば性能が向上するかもしれないと考えて、当日の情報に加えて過去3日分の情報を載せたバージョンでも試してみたが、結果はほとんど変わらなかった。そもそもの素性の出来が悪いのだと思う。
今後の課題としては、今回あまりにも適当に作ってしまった素性をちゃんと考えて作るとか、交差検定を使ってパラメータをチューニングするとか、別のカーネル関数で試してみるとかがあると思う。
他のモデル(ナイーブベイズとかパーセプトロンとか)でも試してみたい。
参考リンク
ソースコード
今回データ整形に使ったコード。
$ python jmaformat.py data.csv > formatted_data.tsv
$ grep ^200 formatted_data.tsv | python format_vector.py >train.tsv
$ grep -v ^200 formatted_data.tsv | python format_vector.py >test.tsv
# coding: utf-8
import codecs
class JMACSVReader(object):
def __init__(self, filename, format, encoding='Shift_JIS', header_ofs=5):
self.format = format
self._istream = codecs.open(filename, 'r', encoding)
for i in range(header_ofs):
self._istream.readline() # ヘッダを読み飛ばす
def _parse_line(self,line):
record = {}
for k,v in zip(self.format, line.split(',')):
if k is not None:
record[k] = v
return record
def __iter__(self):
for line in self._istream:
yield(self._parse_line(line))
weather_codes = [
(u'薄', 'Usu'), (u'大', 'Oo'),
(u'暴風雨', 'Bofu'), (u'快晴', 'Kaisei'),
(u'雨', 'Ame'), (u'晴', 'Hare'),
(u'雪', 'Yuki'), (u'曇', 'Kumori'),
(u'みぞれ', 'Mizore'), (u'ひょう', 'Hyo'),
(u'霧', 'Kiri'),
]
def encode_weather(weather):
# 日本語があると何かと面倒なのでコード値にしてしまう。彦根市は上記のルールだけで大丈夫だったが、他の地域や期間に対応するにはパターンを増やす必要があるかも。
# 「晴れ一時雨後曇り」とか後ろについてくるやつは面倒くさいので除去する。
wcode = re.sub( u'(一時|時々|後|、).*','',weather.strip())
for k,v in weather_codes:
wcode = wcode.replace(k,v)
return wcode
def format_file_main(filename, in_format, out_format):
jmacsv = JMACSVReader(filename, in_format)
for record in jmacsv:
line = []
for of in out_format:
if isinstance(of,tuple):
line.append(of[1](record[of[0]]))
else:
line.append(record[of])
print(u'\t'.join(line).encode('utf-8'))
if __name__=='__main__':
# usage: jmaformat.py data.csv
import sys
import re
in_format = [
'date',
'temp_max', None,None,
'temp_min', None,None,
'rainfall', None,None,None,
'windspeed', None,None,
'atom', None,None,
'humidity', None,None,
'weather', None,None,
]
out_format = ['date', 'temp_max', 'temp_min', 'rainfall', 'windspeed', 'atom', 'humidity', ('weather', encode_weather ), ]
for file in sys.argv[1:]:
format_file_main(file, in_format, out_format)
# coding:utf-8
import datetime
# 訓練データ中にしか出現しない珍しい天気などがあると評価できないので、似た天気を適当に束ねて4クラスにまとめる。
weather_code_map = {
'Kumori': 'Kumori',
'UsuKumori': 'Kumori',
'Kiri': 'Kumori',
'Hare': 'Hare',
'Kaisei': 'Hare',
'Ame': 'Ame',
'OoAme': 'Ame',
'Bofu': 'Ame',
'Yuki': 'Yuki',
'OoYuki': 'Yuki',
'Mizore': 'Yuki',
'Hyo': 'Yuki',
}
def normalize_weather(weather):
return weather_code_map[weather]
previous_lines = []
def create_train_vector(line, degree):
""" 気象情報から、学習用の特徴ベクトルを作る。
特徴量として、当日を含めて過去degree日間の気象情報および天気を盛り込む。
過去のデータがdegree日分溜まっておらずベクトルを作れない場合はNoneを返す。
"""
newline = None
line[-1] = normalize_weather(line[-1])
if len(previous_lines) >= degree:
newline = []
lastdate = datetime.datetime.strptime(previous_lines[-1][0], '%Y/%m/%d') # 当日の日付
newline += [str(lastdate.month),str(lastdate.day)]
for i in range(degree):
newline += previous_lines[-1-i][1:] # 過去n日分の気象情報
newline.append(line[-1]) # 最後に明日の天気を追加
previous_lines.append(line)
return newline
if __name__=='__main__':
# usage: python format_vector.py < formatted_data.tsv
import sys
for line in sys.stdin:
newline = create_train_vector(line.rstrip('\n').split('\t'), 1)
if newline:
print '\t'.join(newline)