HOME > CL LAB > (Japanese text only.) 天気データのみで電力を予測する：線形回帰編

天気データのみで電力を予測する：線形回帰編

2018.03.15.　★ 354

こんにちは、クリエーションラインの朱です。最近はどんな業界でも、どんな会社でもAIという言葉を使い始めましたね。こんな熱いAIの分野で、新人でもありますが、日々精進しています。
今回は「重回帰で時系列データを扱う」というテーマで機械学習の話をしたいと思います。

まず僕のような新人のために重回帰が何かを直感的に簡単に説明します。

／1a+2b=3
＼2a+5b=8
上記のaとbを求めてくださいのような中学ラベルの数学問題が見たことがあると思います。

今回のテーマは天気データのみで電力を予測するということなので、書き換えると以下のテーブルになりますね。

気温	降水量	電力
1	2	3
2	5	8

方程式が以下になります。

a☓気温＋b☓降水量＝電力

統計的手法で上記説明変数「気温」の係数aと説明変数「降水量」の係数bを求めるのが回帰です。ここが多変数なので、重回帰になります。電量は目的変数と呼びます。
統計的手法とは具体的にコスト関数を設定して、コスト関数の値が最小になるようにa,bを調整して行く手法ですが、詳細を説明すると長くなるので割愛します。
上記の方程式をモデルと呼びます。最適なモデルを作ることで、天気予報を方程式に代入することで未来の電力を予測できるわけです。

話を戻します。人間が時間軸で生きていますから、電力データのような時系列データが多いですが、実際に時系列データの分析では、多変数の分析が難しく、単一変数の分析が多いようです。
重回帰を使うと、多変数の時系列データ分析が可能ですが、回帰を正しく行うためには、扱うデータが独立仮説（Assumption of Independence）に従わなければならないです。
時系列データは基本独立仮説に従わないので、従うようにfeature engineeringが必要ですし、データのトレンドを考量した学習アルゴリズムを選択することによって独立仮説の問題を解消することが可能です。

独立仮説とは、データが以下の条件に満たしていると仮定します。
1. 説明同士が独立であります。天気データで例えると、気温と日照量が強く関係しているので、両方を同時に説明変数として使うと、独立仮説に違反します。
2. 一つの説明変数のデータ同士の間が独立であります。例えば明日の来客数が今日の来客数に影響されない等。

但し、データが独立仮説に違反しないかを確認するには、データに関するドメイン知識が必要になりますし、データを可視化し、時系列データの季節調整を行うとか、機械学習の前処理が必要になります。
独立仮説の検証は今回の目的ではありませんので、確認せず、時系列データに対して重回帰を行いました。
では本題に入りたいと思います。

天気データのみで電力を予測する

自分自身は電力系のシステム開発を行った経験がありまして、電力の予測に興味があって、今回のテーマにしました。

早速データを取得します。

電力データと天気データを取得し、教師データ(80%)とテストデータ(20%)にわけます。

import pandas as pd

# tepcoから電力データを取得
energy_data_url="http://www.tepco.co.jp/forecast/html/images/juyo-2017.csv"
energy_data_df = pd.read_csv(energy_data_url, skiprows=[0, 1], encoding="shift-jis")
energy_data_df.columns = ["DATE", "TIME", "energy"]
energy_data_df["time"] = energy_data_df["DATE"] + " " + energy_data_df["TIME"] + ":00"

# 天気データを読みこむ
weather_data_df = pd.read_csv("./sample_data/weather_201709.csv",
                              skiprows=[0, 1, 2, 4],# dont skip 3rd row as it's the header
                              encoding="SHIFT-JIS")

# extract time, rain_amount, temperature data 
weather_data_df = weather_data_df.iloc[:, [0, 1, 7]]
weather_data_df.columns = ["time", "rain_amount", "temperature"]

# merge Energy and Weather data as a data set
data_df = pd.merge(weather_data_df, energy_data_df).reset_index(drop=True)

# split the data set as training(80%) and test(20%)
from sklearn.model_selection import train_test_split
x_train_df, x_test_df, y_train_df, y_test_df = train_test_split(data_df.loc[:,["time", "temperature", "rain_amount"]], 
                                                                data_df.loc[:,["energy"]], 
                                                                test_size=0.2, shuffle=False)
x_test_df = x_test_df.reset_index(drop=True)
y_test_df = y_test_df.reset_index(drop=True)

※ 天気データは気象庁から取得しました。

取得したデータをmatplotlibでプロットします。

def plot_data(time_df, x_df, y_df, start_datetime, end_datetime, title):
    start_pos = 0
    end_pos = 0
    start_datetime_dt = datetime.strptime(start_datetime, "%Y/%m/%d %H:%M")
    end_datetime_dt = datetime.strptime(end_datetime, "%Y/%m/%d %H:%M")
    for index, row in time_df.iterrows():
        temp_datetime_dt = datetime.strptime(row[0], "%Y/%m/%d %H:%M:%S")
        if temp_datetime_dt == end_datetime_dt:
            end_pos = index
        if temp_datetime_dt == start_datetime_dt:
            start_pos = index

    fig, ax1 = plt.subplots(figsize=(len(time_df)/3, 12))
    time_matrix = dates.date2num([datetime.strptime(t, "%Y/%m/%d %H:%M:%S") for t in time_df.as_matrix().transpose()[0]])
    ax1.plot(time_matrix, x_df.loc[:, ["temperature"]], label="temperature", color="black")
    ax1.plot(time_matrix, x_df.loc[:, ["rain_amount"]], label=rain amount, color="maroon")
    ax1.plot(time_matrix, y_df.loc[:, ["energy"]], label="energy", color="crimson")

    ax1.legend(loc=2)
    ax1.set_xticklabels(time_matrix, rotation=90)
    ax1.tick_params(labelsize=18)
    ax1.xaxis.set_major_locator(dates.HourLocator(interval=1))
    hfmt = dates.DateFormatter("%Y/%m/%d %H:%M")
    ax1.xaxis.set_major_formatter(hfmt)
    plt.title(title)
    plt.show()

time_train_df = x_train_df.loc[:, ["time"]]
x_train_df = x_train_df.loc[:, ["temperature", "rain_amount"]]
plot_data(time_train_df, x_train_df, y_train_df, start_datetime="2017/09/01 01:00", 
          end_datetime="2017/09/15 23:00", title="training data before normalization")

うん、これだと天気データと電力データの関係がわかりませんので、正規化したデータを再プロット。

from sklearn import preprocessing
def range_scale_data(dataframe):
    min_max_scaler = preprocessing.MinMaxScaler(feature_range=(-0.5, 0.5))
    return pd.DataFrame(min_max_scaler.fit_transform(dataframe), index=dataframe.index,
                        columns=dataframe.columns)

plot_data(time_train_df, range_scale_data(x_train_df), range_scale_data(y_train_df), start_datetime="2017/09/01 01:00", 
          end_datetime="2017/09/15 23:00", title="training data before normalization")

気温と電力の間に強い関係がみられるのが分かりますね。これが回帰を上手く行うための必須要素です。但し、プロットするまでもないが、冬の気温が低いのに電力が高いでしょう。（暖房を使うから）
無論こういた季節性要素も考量する必要がありますが、今回は使ったデータが1ヶ月分のみなので、季節性を無視しました。

説明変数の正規化

回帰を行うためには説明変数の正規化が必要です。どういう正規化手法を使うかは正直経験則だと思いますが、正規化の目的としては主に以下2点と思います。
1. スケールの大きい説明変数に釣られず、正しい分析するため
2. 学習アルゴリズムのの計算時間を節約するため

今回は全ての説明変数を-0.5~0.5の範囲に再スケーリングしました。

x_train_df = range_scale_data(x_train_df)
x_test_df = range_scale_data(x_test_df)

モデルの作成及び交差検証

今回は線形回帰ということで、モデルの作成はLassoアルゴリズムを使いました。
また、モデルの性能向上のため、交差検証が定番の手法です。
交差検証の手法は沢山ありますが、ここではscikit-learnのShuffleSplitを利用しました。
教師データから30%のデータを交差検証データにしました。

from sklearn import linear_model
from sklearn.model_selection import cross_val_score, ShuffleSplit

cv = ShuffleSplit(n_splits=10, test_size=0.3, random_state=0)
best_lamda = 0.0
prev_score = 0.0
for lamda in [0.0001, 0.0003, 0.001, 0.003, 0.01, 0.03, 0.1, 0.3, 1.0, 3.0, 10]:
    clf_model = linear_model.Lasso(alpha=lamda)
    clf_model.fit(x_train_df.as_matrix(), x_test_df.as_matrix().transpose()[0])
    # use cross validation to get the best lamda
    score = np.mean(cross_val_score(clf_model, x_train_df.as_matrix(), x_test_df.as_matrix(), cv=cv))
    if prev_score < score:
        best_lamda = lamda
        prev_score = score
    # build model with the best lamda
    clf_model = linear_model.Lasso(alpha=best_lamda)

上記のbest_lambaとはLassoのL1 正則化項の係数です。
ここの交差検証がscikit-learnのcross_val_scoreを利用していて、評価指標はR^2らしいです。
R^2の他にも、MAE、RSME等いろんな指標が存在するので、実際にCase-by-caseで使うのが良いでしょう。
clf_modelが実際に作ったモデルになります。

では機械学習定番の学習曲線を実際にプロットしてみました。ソースコードはscikit-learnの例をそのまま使っていましたので、省略します。
まずは今回使ったLassoの学習曲線です。

次に線形アルゴリズムではないが、機械学習の現場でよく使われるRandom Forestの学習曲線です。

鋭い方はもう理解できていると思いますが、教師データ数(training samples)を増やしても今回のLassoの精度が上がらないことが分かりますね。対して、Random Forestが教師データ数(training samples)を増やせば、精度がどんどん上がることが分かります。
但しアルゴリズムの選択は今回の範囲外なので、割愛します。

学習したモデルで電力の予測を行って見ました

import sklearn.metrics as metrics

def evaluate_model(model_name, clf_model, X_test_matrix, Y_test_matrix, columns, title):
    if hasattr(clf_model, "coef_"):
        importances = clf_model.coef_
    else:
        importances = clf_model.feature_importances_
    # plot weights
    if 1 == (importances.shape)[0]:
        importances = importances[0]
    indices = np.argsort(importances)
    plt.figure(figsize=(30, 8))
    plt.subplot(1, 3, 1)
    plt.title("coefficient")
    plt.barh(range(len(indices)), importances[indices], color='b', align='center')
    plt.yticks(range(len(indices)), np.array(list(columns))[indices])
    plt.xticks(fontsize=18)

    # predictions
    predictions = clf_model.predict(X_test_matrix)
    # mean absolute error
    #print(Y_test_matrix)
    #print(predictions)
    mae = metrics.mean_absolute_error(Y_test_matrix, predictions)

    plt.subplot(1, 3, 2)
    plt.title(title + " Accurracy")
    plt.scatter(Y_test_matrix, predictions, label="Predicted")
    if Y_test_matrix.max() > predictions.max():
        plotmax = Y_test_matrix.max()
    else:
        plotmax = predictions.max()
    x = np.linspace(0, plotmax)
    y = x
    plt.plot(x, y, "r-")      # draw a diagonal line
    accuracy = round(metrics.r2_score(Y_test_matrix, predictions) * 100, 2)
    legend = plt.legend(loc=2, title="[{0}:{1:.2f}%]".format(model_name, accuracy))
    plt.setp(legend.get_title(), fontsize='20')
    plt.xlabel("real values")
    plt.ylabel("predicted values")
    plt.xticks(fontsize=18)
    plt.yticks(fontsize=18)
    plt.show()
    return mae

mae_lasso = utilities.evaluate_model("Lasso", clf_model, x_test_df.as_matrix(), y_test_df.as_matrix(), 
                                     ["rain_amount", "temperature", "weekend", "time_zone"], 
                                     "Energy Prediction")
print("MAE of Random Forest:"+ str(float("{0:.5f}".format(mae_lasso))))

MAE of Lasso:340.9775

※　余談ですが、ソフトウェア開発出身の関数/メソットエンジニアなので、習慣的に関数で書いてしまうので、関数が好きではない方には申し訳ない。

気温の係数(coefficient)が確かに高く、電力との関連性が高いことを示していますが、精度のほうがいまいちですね。正規化後のグラフを再び確認してみると、昼と夜はトレンドが違うので、区別したほうが良さそうです。また週末企業が休みなので、電力も変わるでしょう。
このように仮説を立ててみました。当り前の仮説を立ててみましたが、実際結構当り前のことでも思いつかない場合が多いので、データを分析しながら仮説を立てるのが非常に重要です。

では、実際に上記仮説に対して説明変数２つを追加してみました。これがいわゆるfeature engineeringですね。

time_array = weather_data_df["time"].as_matrix()
weekend_array = weather_data_df["time"].as_matrix()

# add weekend as a new feature
for i in range(0, weekend_array.shape[0]):
    weekend_array[i] = float(
        datetime.strptime(weekend_array[i], '%Y/%m/%d %H:%M:%S').weekday())
    if weekend_array[i] == 5 or weekend_array[i] == 6:
        weekend_array[i] = 1
    else:
        weekend_array[i] = 0

for i in range(0, time_array.shape[0]):
    time_array[i] = float(datetime.strptime(time_array[i], '%Y/%m/%d %H:%M:%S').hour)
    if time_array[i] > 6.0:
        time_array[i] = 0.0
    else:
        time_array[i] = 1.0

weather_data_df["weekend"] = weekend_array
weather_data_df["time_zone"] = time_array

追加した特徴も正規化し、再びモデルを作成、予測を行いました。

MAE of Lasso:215.32479

time_zone(day/nightのフラグ)がそれなりに電力と関連しているのが分かりました。予測性能も改善しました。
まだまだ改良が必要ですが、精度が70%を突破しましたので、今回はここまで終了とします。

因みに、気になる方もいらっしゃると思いますが、精度がまあまあ達成していて、MAEの値がありえないぐらい高いのが何故でしょうか。
ここでは説明しませんが、興味がある方は自分で考えてみてください。

終わり

はい、以上で今回は線形回帰で電力を予測してみました。
time_zone(day/nightのフラグ)を特徴に追加するだけでそんなに性能が変わりますね。feature engineeringの重要さは再び感じています。
今回の結果的には予測が大幅にはずれていないが、まだまだfeature engineeringが必要です。

また、前節で紹介した通りデータが独立仮説に従っているかが判断し辛い場合は取り敢えず回帰でやってみるのも手ですね。
これ自体も仮説検証なので、とても重要なステップです。データが独立仮説に従っていないことが分かれば別の手法に切り替えれば良いと思います。

今回のソースを全てgithubに入れているので、興味のある方はご参照ください
https://github.com/george-j-zhu/timeseriesprocessing

次回は時系列データ分析でよく使われる深層学習のLSTMアルゴリズムを試す予定です。