2. アメダス気象データ分析チャレンジ!(Python版) 2日目 実践編

Copyright 2021 気象ビジネス推進コンソーシアム、岐阜大学 吉野純
(C) 2021 WXBC、岐阜大学 吉野純

<利用条件>
本書は、本書に記載した要件・技術・⽅式に関する内容が変更されないこと、および出典 を明示いただくことを前提に、無償でその全部または⼀部を複製、翻案、翻訳、転記、引 用、公衆送信等して利用できます。なお、全体を複製、翻案、翻訳された場合は、本書にあ る著作権表⽰および利用条件を明⽰してください。

<免責事項>
本書の著作権者は、本書の記載内容に関して、その正確性、商品性、利用目的への適合性 等に関して保証するものではなく、特許権、著作権、その他の権利を侵害していないことを 保証するものでもありません。本書の利用により生じた損害について、本書の著作権者は、 法律上のいかなる責任も負いません。

次に,過去の気象データ過去の電力消費量データを分析することで統計的モデル(ニューラルネットワーク)を構築して将来の電力消費量を予測します.このような分析を予測的分析(Predictive Analytics)と言います.予測的分析により,例えば,商品ごとの需要を予測分析することによって,需要が上がるタイミングで商品を発注し,そして需要が下がるタイミングで値下げを行い,企業の潜在的な商機やリスクを特定することが可能となります.前回のアメダス気象データ分析チャレンジ!(Python版)の1日目の理解を前提としていますので,必要があれば前回の資料を復習しながら読み進めてください.

2.1 使用するデータ

2.1.1 気象庁AMeDASの気象データ

気象データビジネスデータにより統計的モデルを作ってみましょう.このプログラムでは,気象データとして,気象庁AMeDASの2017年1月1日1時~2019年1月1日0時の1時間毎2年間分東京気温露点温度を使用します.露点温度は,気体を冷却することで凝結(結露)する温度のことを指しており,空気の湿り具合を表す指標の1つです.前半の1年間つまり2017年1月1日1時~2018年1月1日0時の間のデータは予測の学習データとして,後半の1年間つまり2018年1月1日1時~2019年1月1日0時の間のデータは予測の検証データとして使用されます.ファイル名は,amedas.csvとします.ダウンロードの詳細については,「アメダス気象データ分析チャレンジ!(Python版)の1日目」の講義内容をご覧下さい.

気象庁のホームページのリンク: http://www.data.jma.go.jp/gmd/risk/obsdl/

今回のチャレンジ!では事前に用意したcsv形式のファイル(amedas.csv)を使用しますので,データをダウンロードしていただく必要はありませんが,異なる地点や期間,異なる気象要素を各自でダウンロードして分析してみて下さい.

2.1.2 東京電力の電力消費量データ

また,このプログラムでは,東京電力の電力消費量データを使用します.気象庁AMeDASの気温データと同様に,2017年1月1日0時~2018年12月31日23時の1時間毎2年間分東京電力電力消費量を入手します.前半の1年間つまり2017年1月1日1時~2018年1月1日0時の間のデータは予測の学習データとして,後半の1年間つまり2018年1月1日1時~2019年1月1日0時の間のデータは予測の検証データとして使用されます.ファイル名は,tepco.csvとします.このデータは,統計的モデルにおける目的変数として使われるだけでなく,検証用データ(正解データ)としても使用されます.ダウンロードの詳細については,「アメダス気象データ分析チャレンジ!(Python版)の1日目」の講義内容をご覧下さい.

東京電力でんき予報のリンク: http://www.tepco.co.jp/forecast/index-j.html

今回のチャレンジ!では事前に用意したcsv形式のファイル(tepco.csv)を使用しますので,データをダウンロードしていただく必要はありませんが,異なる地域の電力消費量データやその他のオープンデータを各自でダウンロードして分析してみて下さい.

2.2 Pythonプログラムの構成と実行

このプログラムは,本日配布した資料一式の中のpythonプログラムpredictive_analytics.pyを,Jupyter Notebook用に再収録したものです.どちらも同じ内容のpythonプログラムとなっています.このプログラムは,以下のような構成となっています.

 1. ライブラリのインポート
 2. 気象庁AMeDASの気象データ(csv形式)を読む
 3. 東京電力の電力消費量データ(csv形式)を読む
 4. 2つのデータを結合して整理する
 5. 学習データと検証データの作成
 6. データの標準化
 7. 統計的モデルの設定と最適化
 8. 時系列図と散布図の作成

早速,1つずつ手を動かしてみましょう.

2.2.1 ライブラリのインポート

以下のセルから,プログラムが始まります.上のプルダウンメニューが「Code」となっているセルがpythonプログラムに相当します.セルにカーソルを合わせて上にある「Run」ボタンを押して実行して下さい.以降のプログラムは上のプログラムを前提にして組み立てられてますので,上から順番に実行する必要がありますので注意して下さい.

上の「View」タブの中の「Toggle Line Numbers」をクリックすると,セルに行番号が表示されるので表示させて読み進めて下さい.

まずは,必要なライブラリをインポートします.2行目で数値計算ライブラリnumpyをインポートします.import numpy as npとすることで,プログラム中ではnpと省略してインポートすることができます.作図のための,matplotlibseabornも4行目と6行目でインポートします.また,8行目のdatetimeは日付(日時や時間)の処理ができる便利なライブラリで,気象データ分析には必須となります.また,10行目でデータの前処理や整理に不可欠なpandasもインポートします.機械学習に必要となる標準化などの前処理のためのライブラリpreprocessingを16行目でインポートします.最後に,18行目の線形重回帰モデル(linear_model)と20行目のニューラルネットワーク(neural_network)を作成するのに必要なsklearn (scikit-learn)もインポートします.

Jupyter Notebookの中だけで有効なマジックコマンド%precison 3(小数点3桁まで)と%matplotlib inline(インラインで画像を表示)も忘れずにつけましょう.

In [1]:
# numpyをnpという別名でインポートします.
import numpy as np
# matplotlibをpltという別名でインポートします.
import matplotlib.pyplot as plt
# Seabornをインポートします.
import seaborn as sns
# datetimeは日時データを処理する際に便利なメソッドです.インポートします.
from datetime import datetime
# pandasをpdという別名でインポートします.
import pandas as pd
# matplotlibで時系列図を作成するときには以下をインポートします
from pandas.plotting import register_matplotlib_converters
# これを登録しておきます.
register_matplotlib_converters()
# sklearn(scikit-learn)の前処理のためのライブラリをインポートします
from sklearn import preprocessing
# sklearn(scikit-learn)の線形重回帰モデルをインポートします.
from sklearn import linear_model
# sklearn(scikit-learn)のニューラルネットワークのライブラリをインポートします
from sklearn import neural_network
%precision 3
%matplotlib inline

2.2.2 気象庁AMeDASの気象データ(csv形式)を読む

まずは,前回の診断的分析と同様に,気象庁AMeDASの気象データ(amedas.csv)を読み込んで,pandasのデータフレームとしてメモリすることからはじめます.表計算ソフトのExcelで例えるならば,Excelを起動してcsvファイルを開いた状態です.

このデータを見ると,上から5行目まではデータとは関係のないヘッダー情報となっているようです.6行目から最終列(17525列目)までのデータ部が必要な情報です.また,A列は日時と時刻,B列が気温(℃),C列とD列は気温の品質や均質に関する参考情報となっています.また、E列は露点温度(℃),F列とG列は露点温度の参考情報となっています.今回のデータ分析に必要なデータは3列分(A列,B列,E列)のみです.ここでは,前回の診断的分析に基づき,気温と露点温度が電力消費量の変化に影響するものと仮設を立ててデータを用意しています.そして,分析に不要な情報は削除し,必要な情報だけをメモリする必要があります.

インポートしたpandasには,csvファイルを読み込むメソッドpd.read_csvが備わっています.これを使うとcsvファイルをそのままデータフレームに割り当てることができます.後のデータ整理ではpandasのデータフレームを使うと便利なので好都合です.ここでは,以下のセルのようにreadamedasという関数(3行目)の中でデータの読み込みを行います.この関数は分かりやすく複数行に別けて書かれていますが,メソッドpd.read_csvを使えば,実質1行分のプログラミングでデータを読み込むことができます.このpd.read_csv(13行目)の引数として,14~19行目の6種類を入力する必要があります.

pd.read_csv(

1. ファイル名(ここでは,amedas.csv)
2. encoding:文字エーコーディング(ここでは,Shift_JIS)
3. skiprows:上から何行分を読み飛ばすか(ここでは,5行)
4. names:列のラベル(ここでは,A列~D列を['date_time', 'T', 'dummy1', 'dummy2', 'Td', 'dummy3', 'dummy4']とつけます)
5. parse_dates:datetime型で読み込む列(ここでは,'date_time'の列をdatetime型で読み込む)
6. index_col:インデックスとする列(ここでは,'datetime')

)

この様にして読み込んだデータは,pandasのデータフレーム型で割り当てられて,関数の中のdfという変数から,メインプログラムではamedasという変数に戻されて格納されます(29行目).その直後の31行目では,

amedas=amedas.drop(['dummy1', 'dummy2', 'dummy3', 'dummy4'], axis=1)

のメソッドdropにより,ラベルが'dummy1'~'dummy4'の列が削除されています.axis=1とすると列の方向で削除されるという意味です.axis=0とすれば行の方向で削除されます.

それでは実行してみましょう.

In [2]:
# 気象庁アメダスの気温の時系列データを読み込んで,
# DataFrameに割り当てる関数
def readamedas(filename,skipline): 
    # pandasのread_csvというメソッドでcsvファイルを読み込みます.
    # 引数として,
    # [0]入力ファイル名
    # [1]エンコーディング
    # [2]読み飛ばす行数,
    # [3]column名
    # [4]datetime型で読み込むcolumn名
    # [5]indexとするcolumn名
    # を与える.
    df=pd.read_csv(
       filename,
       encoding='Shift_JIS',
       skiprows=skipline,
       names=['date_time','T','dummy1','dummy2','Td','dummy3','dummy4'],
       parse_dates={'datetime':['date_time']}, 
       index_col='datetime'
       )
    return df

# 気象庁AMeDAS(東京)の気温(℃)と露点温度(℃)のcsvファイルの名前を指定します.
filename1='amedas.csv'
# csvファイルの上から5行目まではデータではないため呼び飛ばします.
skipline1=5
# 気象庁AMeDAS(東京)の気温(℃)と露点温度(℃)のcsvファイルを読み込んで,
# pandasのDataFrame(tepco)に割り当てる関数を呼び出します.
amedas=readamedas(filename1,skipline1)
# DataFrame(amedas)の中のdummy1~dummy7の列を削除する.
amedas=amedas.drop(['dummy1','dummy2','dummy3','dummy4'],axis=1)
amedas
Out[2]:
T Td
datetime
2017-01-01 01:00:00 5.1 0.7
2017-01-01 02:00:00 4.1 -1.3
2017-01-01 03:00:00 4.0 -0.8
2017-01-01 04:00:00 3.0 -1.0
2017-01-01 05:00:00 3.6 -0.8
... ... ...
2018-12-31 20:00:00 5.3 -6.3
2018-12-31 21:00:00 4.8 -6.7
2018-12-31 22:00:00 4.2 -6.1
2018-12-31 23:00:00 4.6 -6.3
2019-01-01 00:00:00 2.7 -5.7

17520 rows × 2 columns

日付('datetime')と気温('T')と露点温度('Td')からなる表が現れれば成功です.dummy1dummy4も消えていることを確認して下さい.これは,amedasというデータフレームの内容を省略して表示しています.2017年1月1日01時00分00秒~2019年1月1日0時0分0秒までの1時間毎の気温データ(合計17520行)を読み込めました.

2.2.3 東京電力の電力消費量データ(csv形式)を読む

次に,同様に,東京電力でんき予報の電力消費量データ(tepco.csv)を読み込んで,pandasのデータフレームとしてメモリしてみましょう.これも,表計算ソフトのExcelで例えるならば,Excelを起動してcsvファイルを開いた状態となります.

このデータを見ると,上から3行目まではデータとは関係のないヘッダー情報となっているようです.4行目から最終列(17523列目)までのデータ部が必要な情報です.また,A列は日時,B列が時刻,C列が電力消費量となっています.気温データのように不必要な列('dummy1'や'dummy2')はありませんが,気温データと違って日時と時刻が2つの列に別れています.こういったデータの特性を考慮しながら情報をメモリする必要があります.

同じく,csvファイルを読み込むメソッドpd.read_csvを使って,csvファイルをそのままデータフレームに割り当てます.ここでは,以下のセルのようにreadtepcoという関数(3行目)の中でデータの読み込みを行います.この関数も分かりやすく複数行で書かれていますが,メソッドpd.read_csvを使えば,実質1行分のプログラミングでデータを読み込むことができます.このpd.read_csv(13行目)の引数として,次の6種類を設定する必要があります.

pd.read_csv(

1. ファイル名(ここでは,tepco.csv)
2. encoding:文字エーコーディング(ここでは,Shift_JIS)
3. skiprows:上から何行分を読み飛ばすか(ここでは,3行)
4. names:列のラベル(ここでは,A列~C列を['date','time','POWER']とする)
5. parse_dates:datetime型で読み込む列(ここでは,'date'と'time'の2つの列をまとめてdatetime型で読み込む)
6. index_col:インデックスとする列(ここでは,'datetime')

)

この様にして読み込んだデータは,pandasのデータフレーム型で割り当てられて,関数の中のdfという変数から,メインプログラムではtepcoという変数に戻されて格納されます(29行目).

それでは実行してみましょう.

In [3]:
# 東京電力の電力消費量の時系列データを読み込んで,
# DataFrameに割り当てる関数
def readtepco(filename,skipline): 
    # pandasのread_csvというメソッドでcsvファイルを読み込みます.
    # 引数として,
    # [0]入力ファイル名
    # [1]エンコーディング
    # [2]読み飛ばす行数,
    # [3]column名
    # [4]datetime型で読み込むcolumn名
    # [5]indexとするcolumn名
    # を与える.
    df=pd.read_csv(
       filename,
       encoding='Shift_JIS',
       skiprows=skipline,
       names=['date','time','POWER'],
       parse_dates={'datetime':['date','time']},
       index_col='datetime'
       )
    return df

# 東京電力の電力消費量(10^6 kW)のcsvファイルの名前を指定します.
filename2='tepco.csv'
# csvファイルの上から3行目まではデータではないため呼び飛ばします.
skipline2=3
# 東京電力の電力消費量(10^6 kW)のcsvファイルを読み込んで,
# pandasのDataFrame(tepco)に割り当てる関数を呼び出します.
tepco=readtepco(filename2,skipline2)
tepco
Out[3]:
POWER
datetime
2017-01-01 00:00:00 2783
2017-01-01 01:00:00 2634
2017-01-01 02:00:00 2520
2017-01-01 03:00:00 2438
2017-01-01 04:00:00 2389
... ...
2018-12-31 19:00:00 3531
2018-12-31 20:00:00 3473
2018-12-31 21:00:00 3376
2018-12-31 22:00:00 3252
2018-12-31 23:00:00 3198

17520 rows × 1 columns

日付('datetime')と電力消費量('POWER')からなる表が現れれば成功です.これは,tepcoというデータフレームの内容を省略して表示しています.2017年4月1日00時00分00秒~2018年12月31日23時0分0秒までの1時間毎の電力消費量データ(合計17520行分)を読み込めました.

2.2.4 2つのデータを結合して整理する

次に,気象データの時系列amedas電力消費量データの時系列tepco結合して,目的変数と説明変数を1つのデータフレームdataに統合します.例えば,表計算ソフトExcelで例えるならば,1つの表からもう1つの表へ時間軸が合うようにコピーした状態であると言えます.また,このデータフレームから,気温('T')と露点温度('Td')から湿数('T-Td') を計算して,また,時刻日にち平日・休日1日前の発電量1時間前の発電量2時間前の発電量3時間前の発電量といったデータを抽出して,説明変数の1つとして結合します.以前の診断的分析の結果により,これらの情報の一部は電力消費量の変化と関係していることが想像されます.気温('T')のみならず空気の湿り具合(湿数'T-Td')も電力消費量に関係しているものと想像されますので湿数も説明変数として準備しましょう.以下のセルの4行目からの関数dataprocessの中でその一連の処理を行っています.

まず,pandasには2つのデータフレームを結合するメソッドpd.merge()があります.次のセルの6行目にあるように,

df=pd.merge(y_data, x_data, how="outer", left_index=True, right_index=True)

と実行することで,y_data (tepco)のインデックスdatetimeとx_data (amedas)のインデックスdatetimeと紐付けながら,2つのデータを横方向に結合することができます.howは,outerとすれば全てのインデックスが残るように結合することになり(完全外部結合)innerとすれば両方のデータに含まれるインデックスだけが残るように結合することになります(内部結合)left_index=Trueと設定すると,左側のデータフレームのインデックスを結合のキーとして用います.right_index=Trueと設定すると,右側のデータフレームのインデックスを結合のキーとして用います.

次に8行目のdropna()メソッドにより,上記で結合されたデータフレームdataから欠測値NaNを取り除きますhow='any'は'T'と'Td'と'POWER'のどれかひとつに欠測がある時刻(行)は削除されることになります.一方,how=allの場合には,'T'と'Td'と'POWER'のすべてが欠測となる時刻(行)が削除されることになります.

次に,10行目では,気温('T')と露点温度('Td')の差を取ることで湿数('T-Td')を計算し,新たな列として追加しています.湿数(単位:℃)とは,空気の湿り具合を表す指標であり,0℃に近いほど湿潤であり(相対湿度100%に近づく),0℃より大きくなることで乾燥する(相対湿度0%へと近づく)ことを意味します.空気の湿り具合を表す1つの指標として湿数を導入しています.露点温度('Td')のデータは使いませんので,11行目でdrop()メソッドで削除しています.

ここからは,データフレームのdatetimeの情報の中から,説明変数として使えそうな情報を抽出します.16行目ではデータフレームdfのインデックスから .hour により 時刻 を取り出し,キーを'HOUR'として,hourというデータフレームに格納しています.

hour=pd.DataFrame({'HOUR': df.index.hour}, index=df.index)

このデータフレームhourでも,データフレームdfのインデックスdf.indexを割り当てています.17~24行目では,得られた各時刻のHOURから3時ごとの時間区分(キー'HOUR1','HOUR2',…,'HOUR8')に分類し,該当すれば1,該当しなければ0としています.

hour['HOUR1']=( ( hour['HOUR'] >= 0 ) & ( hour['HOUR'] < 3 ) ).astype(int)
hour['HOUR2']=( ( hour['HOUR'] >= 3 ) & ( hour['HOUR'] < 6 ) ).astype(int)
hour['HOUR3']=( ( hour['HOUR'] >= 6 ) & ( hour['HOUR'] < 9 ) ).astype(int)
hour['HOUR4']=( ( hour['HOUR'] >= 9 ) & ( hour['HOUR'] < 12 ) ).astype(int)
hour['HOUR5']=( ( hour['HOUR'] >= 12 ) & ( hour['HOUR'] < 15 ) ).astype(int)
hour['HOUR6']=( ( hour['HOUR'] >= 15 ) & ( hour['HOUR'] < 18 ) ).astype(int)
hour['HOUR7']=( ( hour['HOUR'] >= 18 ) & ( hour['HOUR'] < 21 ) ).astype(int)
hour['HOUR8']=( ( hour['HOUR'] >= 21 ) & ( hour['HOUR'] < 24 ) ).astype(int)

その後,25行目では,不要になった'HOUR'の列をデータフレームhourから.drop()メソッドで削除し,27行目では,pd.merge()メソッドによりdfhourを結合しています.

次に,32行目では,データフレームdfのインデックスから .day により 日にち を取り出し,キーとして'DAY'を割り当て,dayというデータフレームに格納しています.

day=pd.DataFrame({'DAY': df.index.day}, index=df.index)

さらに,33~35行目では,得られた各時刻のDAYから上旬(キー'DAY1'),中旬(キー'DAY2'),下旬(キー'DAY3')に分類し,該当すれば1,該当しなければ0としています.

day['DAY1']=( ( day['DAY'] >= 1 ) & ( day['DAY'] < 11 ) ).astype(int)
day['DAY2']=( ( day['DAY'] >= 11 ) & ( day['DAY'] < 21 ) ).astype(int)
day['DAY3']=( ( day['DAY'] >= 21 ) & ( day['DAY'] < 32 ) ).astype(int)

その後,36行目では,不要になった'DAY'の列をデータフレームdayから.drop()メソッドで削除し,38行目では,pd.merge()メソッドによりdfdayを結合しています.

また,41行目では,データフレームdfのインデックスから .weekday により 曜日 を取り出し,キーとしてWEEKを割り当て,weekというデータフレームに格納しています. .weekdayは,月曜日は0,火曜日は1,水曜日は2,木曜日は3,金曜日は4,土曜日は5,日曜日は6と返します.

week=pd.DataFrame({'WEEK': df.index.weekday}, index=df.index)

さらに,42行目では,得られた各時刻のWEEKから平日(月0,火1,水2,木3,金4)と休日(土5,日6)に分類して,平日ならば0,休日ならば1となるような修正を行います.

week['WEEK']=( week['WEEK'] >= 5 ).astype(int)

その後,44行目では,pd.merge()メソッドによりdfweekを結合しています.

最後に,49行目では,データフレームdfのインデックスから .month により を取り出し,キーとしてMONTHを割り当て,monthというデータフレームに格納しています.

month=pd.DataFrame({'MONTH': df.index.month}, index=df.index)

さらに,50~53行目では,得られた各時刻のMONTHから12月・1月・2月(キー'DJF'),3月・4月・5月(キー'MAM'),6月・7月・8月(キー'JJA'),9月・10月・11月(キー'SON')に分類して,それぞれ該当すれば1,該当しなければ0とします.

month['DJF']=( ( month['MONTH'] == 12 ) | ( month['MONTH'] == 1 ) | ( month['MONTH'] == 2 ) ).astype(int)
month['MAM']=( ( month['MONTH'] == 3 ) | ( month['MONTH'] == 4 ) | ( month['MONTH'] == 5 ) ).astype(int)
month['JJA']=( ( month['MONTH'] == 6 ) | ( month['MONTH'] == 7 ) | ( month['MONTH'] == 8 ) ).astype(int)
month['SON']=( ( month['MONTH'] == 9 ) | ( month['MONTH'] == 10 ) | ( month['MONTH'] == 11 ) ).astype(int)

その後,54行目では,不要になった'MONTH'の列をデータフレームmonthから.drop()メソッドで削除し,56行目では,pd.merge()メソッドによりdfmonthを結合しています.

また,57~61行目では,目的変数の電力消費量('POWER')のデータから1日前の発電量(power24, キー:'POWER24') を説明変数として加えています.電力消費量を予測するモデルを作ろうとしているのに,電力消費量を説明変数として加えるのはおかしいのではないかと感じる人もいるかもしれませんが,1日前の過去の電力消費量は予測時点では既知のデータとなっていると仮定して,過去の電力消費量をも目的変数としているわけです.同じ要領で,62~66行目では1時間前の発電量(power1,キー:'POWER1'),67~71行目では2時間前の発電量(power2,キー:'POWER2'),72~74行目では3時間前の発電量(power3,キー:'POWER3') といった説明変数をデータフレームdfに追加しています.ちなみに,「1時間前の発電量」 は,目的変数であるキー'POWER'の列データを1つ下にずらすことで作成できます.1つ下にずらすには,pandasのshift()メソッド を使います(64行目).つまり,

df['POWER'].shift(+1)

とすることで,'POWER'列のデータを1つ下にずらした列データとなります.

最後に,再び,76行目ではdropna()メソッドにより,データフレーム中の欠測値(NaN)のある行を削除します(1つ下にずらすことで上端で欠測が生じるため).

それでは,以下のセルを実行してみましょう.

In [4]:
# 2つのデータフレームを結合し,欠測値を除去し,予測モデルの説明変数を
# 列として追加する関数.ここでは,湿数,時刻,日にち,曜日,月,過去の発電量
# を抽出・計算し,データフレームdfに連結している.
def dataprocess(x_data, y_data):
    # 2つのデータフレーム(x_dataとy_data)を結合してデータフレームdfとします.
    df=pd.merge(y_data, x_data, how="outer", left_index=True, right_index=True)
    # NaNがある行を取り除く
    df=df.dropna(how='any')
    # T-Tdにより湿数を計算して、T-Tdの列を作成し、もう使用しないTdは削除する.
    df['T-Td']=df['T']-df['Td']
    df=df.drop(['Td'],axis=1)
    # indexからhourを取り出しデータフレームhour(キー:'HOUR')する.
    # 3時間ごとに時刻の区分(キー)をHOUR1, HOUR2, …, HOUR8に分ける.
    # 該当すれば1,該当しなければ0とする.後の学習の説明変数で使用する.
    # もう使用しないHOURは削除する.
    hour=pd.DataFrame({'HOUR': df.index.hour}, index=df.index)
    hour['HOUR1']=( ( hour['HOUR'] >= 0 ) & ( hour['HOUR'] < 3 ) ).astype(int)
    hour['HOUR2']=( ( hour['HOUR'] >= 3 ) & ( hour['HOUR'] < 6 ) ).astype(int)
    hour['HOUR3']=( ( hour['HOUR'] >= 6 ) & ( hour['HOUR'] < 9 ) ).astype(int)
    hour['HOUR4']=( ( hour['HOUR'] >= 9 ) & ( hour['HOUR'] < 12 ) ).astype(int)
    hour['HOUR5']=( ( hour['HOUR'] >= 12 ) & ( hour['HOUR'] < 15 ) ).astype(int)
    hour['HOUR6']=( ( hour['HOUR'] >= 15 ) & ( hour['HOUR'] < 18 ) ).astype(int)
    hour['HOUR7']=( ( hour['HOUR'] >= 18 ) & ( hour['HOUR'] < 21 ) ).astype(int)
    hour['HOUR8']=( ( hour['HOUR'] >= 21 ) & ( hour['HOUR'] < 24 ) ).astype(int)
    hour=hour.drop(['HOUR'],axis=1)
    # データフレームdfにデータフレームhour列を追加する
    df=pd.merge(df, hour, how="outer", left_index=True, right_index=True)
    # indexからdayを取り出しデータフレームday(キー:'DAY')とする.
    # 10日間ごとに日付けの区分(キー)をDAY1, DAY2, DAY3に分ける.
    # 該当すれば1,該当しなければ0とする.後の学習の説明変数とする.
    # もう使用しないDAYは削除する.
    day=pd.DataFrame({'DAY': df.index.day}, index=df.index)
    day['DAY1']=( ( day['DAY'] >= 1 ) & ( day['DAY'] < 11 ) ).astype(int)
    day['DAY2']=( ( day['DAY'] >= 11 ) & ( day['DAY'] < 21 ) ).astype(int)
    day['DAY3']=( ( day['DAY'] >= 21 ) & ( day['DAY'] < 32 ) ).astype(int)
    day=day.drop(['DAY'],axis=1)
    # データフレームdfにデータフレームday列を追加する
    df=pd.merge(df, day, how="outer", left_index=True, right_index=True)
    # indexからweekdayを取り出しデータフレームweek(キー:'WEEK')とする.
    # 土日(5,6)ならば1、平日(0,1,2,3,4,)ならば0とする。 
    week=pd.DataFrame({'WEEK': df.index.weekday}, index=df.index)
    week['WEEK']=( week['WEEK'] >= 5 ).astype(int)
    # データフレームdfにデータフレームweek列を追加する
    df=pd.merge(df, week, how="outer", left_index=True, right_index=True)
    # indexからmonthを取り出しデータフレームmonth(キー:'MONTH')とする.
    # 3ヶ月ごとに月の区分(キー)をDJF, MAM, JJA, SONに分ける.
    # 該当すれば1,該当しなければ0とする.後の学習の説明変数とする.
    # もう使用しないMONTHは削除する.
    month=pd.DataFrame({'MONTH': df.index.month}, index=df.index)
    month['DJF']=( ( month['MONTH'] == 12 ) | ( month['MONTH'] ==  1 ) | ( month['MONTH'] ==  2 ) ).astype(int)
    month['MAM']=( ( month['MONTH'] ==  3 ) | ( month['MONTH'] ==  4 ) | ( month['MONTH'] ==  5 ) ).astype(int)
    month['JJA']=( ( month['MONTH'] ==  6 ) | ( month['MONTH'] ==  7 ) | ( month['MONTH'] ==  8 ) ).astype(int)
    month['SON']=( ( month['MONTH'] ==  9 ) | ( month['MONTH'] == 10 ) | ( month['MONTH'] == 11 ) ).astype(int)
    month=month.drop(['MONTH'],axis=1)
    # データフレームdfにデータフレームmonth列を追加する
    df=pd.merge(df, month, how="outer", left_index=True, right_index=True)
    # indexから'POWER'を取り出し24時間後にずらしてデータフレームpower24(キー:'POWER24')とする.
    # 後の学習の説明変数とする. 
    power24=pd.DataFrame({'POWER24': df['POWER'].shift(+24)}, index=df.index)
    # データフレームdfにデータフレームpower3列を追加する
    df=pd.merge(df, power24, how="outer", left_index=True, right_index=True)
    # indexから'POWER'を取り出し1時間後にずらしてデータフレームpower1(キー:'POWER1')とする.
    # 後の学習の説明変数とする. 
    power1=pd.DataFrame({'POWER1': df['POWER'].shift(+1)}, index=df.index)
    # データフレームdfにデータフレームpower1列を追加する
    df=pd.merge(df, power1, how="outer", left_index=True, right_index=True)
    # indexから'POWER'を取り出し2時間後にずらしてデータフレームpower2(キー:'POWER2')とする.
    # 後の学習の説明変数とする. 
    power2=pd.DataFrame({'POWER2': df['POWER'].shift(+2)}, index=df.index)
    # データフレームdfにデータフレームpower2列を追加する
    df=pd.merge(df, power2, how="outer", left_index=True, right_index=True)
    # indexから'POWER'を取り出し3時間後にずらしてデータフレームpower3(キー:'POWER3')とする.
    # 後の学習の説明変数とする. 
    power3=pd.DataFrame({'POWER3': df['POWER'].shift(+3)}, index=df.index)
    # データフレームdfにデータフレームpower3列を追加する
    df=pd.merge(df, power3, how="outer", left_index=True, right_index=True)
    # NaNがある行を取り除く 
    df=df.dropna(how='any')
    return df

# 2つのデータフレーム(amedasとtepco)を結合して,
# 欠測値を除去し,目的変数と説明変数の対応関係を示す表を作成する.
data=dataprocess(amedas,tepco)
data
Out[4]:
POWER T T-Td HOUR1 HOUR2 HOUR3 HOUR4 HOUR5 HOUR6 HOUR7 ... DAY3 WEEK DJF MAM JJA SON POWER24 POWER1 POWER2 POWER3
datetime
2017-01-02 01:00:00 2383.0 4.8 6.2 1 0 0 0 0 0 0 ... 0 0 1 0 0 0 2634.0 2494.0 2669.0 2780.0
2017-01-02 02:00:00 2335.0 4.3 5.0 1 0 0 0 0 0 0 ... 0 0 1 0 0 0 2520.0 2383.0 2494.0 2669.0
2017-01-02 03:00:00 2296.0 4.3 4.8 0 1 0 0 0 0 0 ... 0 0 1 0 0 0 2438.0 2335.0 2383.0 2494.0
2017-01-02 04:00:00 2273.0 5.0 4.2 0 1 0 0 0 0 0 ... 0 0 1 0 0 0 2389.0 2296.0 2335.0 2383.0
2017-01-02 05:00:00 2308.0 4.5 3.8 0 1 0 0 0 0 0 ... 0 0 1 0 0 0 2394.0 2273.0 2296.0 2335.0
... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ...
2018-12-31 19:00:00 3531.0 5.8 12.5 0 0 0 0 0 0 1 ... 1 0 1 0 0 0 3683.0 3574.0 3516.0 3190.0
2018-12-31 20:00:00 3473.0 5.3 11.6 0 0 0 0 0 0 1 ... 1 0 1 0 0 0 3656.0 3531.0 3574.0 3516.0
2018-12-31 21:00:00 3376.0 4.8 11.5 0 0 0 0 0 0 0 ... 1 0 1 0 0 0 3554.0 3473.0 3531.0 3574.0
2018-12-31 22:00:00 3252.0 4.2 10.3 0 0 0 0 0 0 0 ... 1 0 1 0 0 0 3373.0 3376.0 3473.0 3531.0
2018-12-31 23:00:00 3198.0 4.6 10.9 0 0 0 0 0 0 0 ... 1 0 1 0 0 0 3211.0 3252.0 3376.0 3473.0

17464 rows × 23 columns

データフレームdataに,予測モデルの作成に必要となる目的変数('POWER') と多数の 説明変数('T', 'T-Td', 'HOUR1', 'HOUR2', 'HOUR3', 'HOUR4', 'HOUR5', 'HOUR6', 'HOUR7', 'HOUR8', 'DAY1', 'DAY2', 'DAY3', 'WEEK', 'DJF', 'MAM', 'JJA', 'SON', 'POWER24','POWER1','POWER2','POWER3') を1つに集約することができました.次の節では,これらの説明変数から予測モデルの学習用データや検証用データを選択する方法について説明します.

2.2.5 学習データと検証データの作成

次に,作成されたデータフレームdataから予測モデルを作成する際に使用する説明変数目的変数を設定し,また,前半1年間分のデータを学習用データ(x_train,y_train)として,後半1年間分のデータを検証用データ(x_test,y_test)としてnp.ndarrayのfloat型で割り当てます.まず,22~42行目のリストx_colsで予測モデルの説明変数を設定しています.ここで1つを設定すれば単回帰モデルになりますが,複数を設定すれば重回帰モデルとなります.リストx_cols.extend()メソッドにより,必要な説明変数が追加されていきます.デフォルトでは,3つの説明変数

x_cols=['T','WEEK', 'HOUR1', 'HOUR2', 'HOUR3', 'HOUR4', 'HOUR5', 'HOUR6', 'HOUR7', 'HOUR8']

が設定されていることになります.前回の診断的分析による知見に基づいて,気温('T')曜日('WEEK')時刻( 'HOUR1', 'HOUR2', 'HOUR3', 'HOUR4', 'HOUR5', 'HOUR6', 'HOUR7', 'HOUR8') が指定されています.次に,45行目のリスト y_colsで予測モデルの目的変数 を設定しています.

y_cols=['POWER']

このような回帰問題では目的変数は1つだけ指定します.また,47行目では,前半1年間の学習用データと後半1年間の検証用データの境目となる行番号dt1を,以下の

dt1=int(len(data)/2)

により取得します.まず,len(data)で全体の行数を取得 して,それを 2で割った数(int型) を割り当てています.

上記3つのデータを引数として,2行目以降の関数preprocessで,まず,説明変数xと目的変数yを設定しています.5行目では,データフレームdfdataのこと)から先に設定したリストx_colsに該当する列の値(.values)を抽出し,np.ndarrayのfloat(浮動小数点)型xとして割り当てています.また,9行目では,データフレームdfdataのこと)からリストy_colsに該当する列の値(.valuesメソッド)を抽出し,同様にnp.ndarrayのfloat(浮動小数点)型yとして割り当てています.ただし,後の予測モデル作成の際に,変数yは行データである必要があるため,.ravel()メソッドにより1次元(ベクトル) に変換しています.

抽出された説明変数xと目的変数yに対して,numpyスライスの機能 を使って,学習用データ(前半)と検証用データ(後半)の境目の行数splitdt1のこと)を指定することにより,11行目と13行目のように,

x_train= x[:split]
y_train= y[:split]

とすることで,インデックスが先頭行~split-1行までが学習用データとなります.一方で,15行目と17行目のように,

x_test= x[split:]
y_test= y[split:]

とすることで,インデックスがsplit行~最終行までが検証用データとなります.

それでは以下のセルを実行してみましょう.

In [5]:
# データを学習用(train)と検証用(test)に別ける.
def preprocess(df, x_cols, y_cols, split):
    # データフレームdfからx_colsの列のデータの値(.values)を
    # float(浮動小数点型)として取り出す.
    x=df.loc[:,x_cols].values.astype('float')
    # データフレームdfからy_colsの列のデータの値(.values)を
    # float(浮動小数点型)として取り出す.
    # ravel()メソッドで1次元化する.
    y=df.loc[:,y_cols].values.astype('float').ravel()
    # 0からsplit-1までのxを学習用データとする
    x_train= x[:split]
    # 0からsplit-1までのyを学習用データとする.
    y_train= y[:split]
    # splitから終わりまでのxを検証用データとする.
    x_test= x[split:]
    # splitから終わりまでのyを検証用データとする.
    y_test= y[split:]
    return x_train,y_train,x_test,y_test

# データフレームdataの中から予測モデルの説明変数とするものを選ぶ
# さまざまな説明変数で試してみましょう.
x_cols=[]
# 気温
x_cols.extend(['T'])
# 湿数
#x_cols.extend(['T-Td'])
# 曜日
x_cols.extend(['WEEK'])
# 時間
x_cols.extend(['HOUR1','HOUR2','HOUR3','HOUR4','HOUR5','HOUR6','HOUR7','HOUR8'])
# 日にち
#x_cols.extend(['DAY1','DAY2','DAY3'])
# 月
#x_cols.extend(['DJF','MAM','JJA','SON'])
# 1日前の消費電力
#x_cols.extend(['POWER24'])
# 1時間前の消費電力
#x_cols.extend(['POWER1'])
# 2時間前の消費電力
#x_cols.extend(['POWER2'])
# 3時間前の消費電力
#x_cols.extend(['POWER3'])
# データフレームdataの中から予測モデルの目的変数とするものを選ぶ.
# ここでは発電量データ('POWER')を目的変数とする.
y_cols=['POWER']
# 学習用データと検証用データに別けるために,全体の行(len(data))を2で割る
dt1=int(len(data)/2)
# データフレームdataを訓練データ(x_train,y_train)と
# 検証データ(x_test,y_test)に別ける
# dt1より前半(主に2017年)のデータを訓練用データとする.
# dt1より後半(主に2018年)のデータを検証用データとする.
# リストx_colsで設定した列が説明変数として扱われる.
# リストy_colsで設定した列が目的変数として扱われる.
x_train, y_train, x_test, y_test=preprocess(data, x_cols, y_cols, dt1)
print('x_train.shape=',x_train.shape)
print('y_train.shape=',y_train.shape)
print('x_test.shape=',x_test.shape)
print('y_test.shape=',y_test.shape)
x_train.shape= (8732, 10)
y_train.shape= (8732,)
x_test.shape= (8732, 10)
y_test.shape= (8732,)

割り当てられたデータのサイズを.shapeメソッドで確認してみましょう.まず,説明変数xについては,学習用データx_train と検証用データx_testに対しては,デフォルト設定で8732行10列のデータが割り当てられているはずです(ただし,x_colsの設定を変えれば列数が変わります).つまりこの場合,入力層は10ニューロンです.また,目的変数yについては,学習用データy_trainと検証用データy_testに対しては,8732行のデータが割り当てられていれば問題ありません.つまり,出力層は1ニューロンです.

2.2.6 データの標準化

次に,データの標準化を行います.機械学習により予測モデルを作成する場合,入力するデータのスケール変換が必要となります.異なる桁数のデータを扱う場合にはスケール変換はほぼ必須となります.ニューラルネットワークやサポートベクターマシンではこの種の前処理をしないと学習の進みが遅いという問題もあります.そこで,機械学習ライブラリsklearn(scikit-learn)に実装されている前処理ライブラリpreprocessingを利用します.preprocessingには,様々なスケール変換のクラスが用意されています.例えば,

・StandardScaler: 標準化(平均0, 標準偏差1)
 $X$は変換前の元データ,$Y$はスケール変換後のデータ,$\mu$は平均値,$\sigma$は標準偏差(データの散らばり具合の尺度).データをの平均値が0,標準偏差が1になるように加工します.大部分のデータは-1から1の範囲に収まりますが,一部のデータはここからはみ出ます.この処理によりデータが正規分布になるわけではない点に注意を要します.
\begin{align} Y = \frac{X- \mu}{ \sigma } \end{align}

・RobustScaler: 外れ値に頑健な標準化
 $X$は変換前の元データ,$Y$はスケール変換後のデータ,$Q_1$は第1四分位(上位1/4位の値),$Q_2$は第2四分位(中央値),$Q_3$は第3四分位(下位1/4位の値).外れ値の影響を受けにくくなるような形でデータを加工します.
\begin{align} Y = \frac{X- Q_2 }{ Q_3 - Q_1 } \end{align}

・MinMaxScaler: 正規化(最大1, 最小0)
 $X$は変換前の元データ,$Y$はスケール変換後のデータ,$x_{max}$は最大値,$x_{min}$は最小値.データを0から1の範囲になるようにデータを加工します.
\begin{align} Y = \frac{X- x_{min} }{ x_{max} - x_{min} } \end{align}

などがあります.ここでは,平均0,標準偏差1となるようなStandardScaler(標準化)preprocessing.StandardScalerクラスを利用しましょう.まず,以下のセルのように,3行目のように,

scaler = preprocessing.StandardScaler()

scalerというインスタンスを作成します.そして,5行目では,

scaler.fit(x_train)

のように,.fit()メソッドにより,説明変数x_trainに対して平均$\mu$と標準偏差$\sigma$を計算して記憶しています.この時点ではまだスケール変換は行われていません.そして,8行目と12行目で,

x_train=scaler.transform(x_train)
x_test=scaler.transform(x_test)

とすることで,記憶された平均と標準偏差からx_trainとx_testに対してそれぞれ標準化を行います.

それでは以下のセルを実行してみましょう.

In [6]:
# 説明変数に対して,平均0,標準偏差1となるような標準化を行う.
# preprocessing.StandardScaler()からscalerというインスタンスを作成する.
scaler = preprocessing.StandardScaler()
# fitメソッドにより説明変数x_trainの平均と分散を計算して記憶する.
scaler.fit(x_train)
# 説明変数x_trainに対して標準化を行い,変換後の配列を返す.
print("original x_train=", x_train)
x_train=scaler.transform(x_train)
print("scaled x_train=", x_train)
# 説明変数x_testに対して標準化を行い,変換後の配列を返す.
print("original x_test=", x_test)
x_test=scaler.transform(x_test)
print("scaled x_test=", x_test)
original x_train= [[4.8 0.  1.  ... 0.  0.  0. ]
 [4.3 0.  1.  ... 0.  0.  0. ]
 [4.3 0.  0.  ... 0.  0.  0. ]
 ...
 [3.3 1.  0.  ... 0.  0.  1. ]
 [2.7 1.  0.  ... 0.  0.  1. ]
 [2.7 1.  0.  ... 0.  0.  1. ]]
scaled x_train= [[-1.339 -0.633  2.648 ... -0.378 -0.378 -0.378]
 [-1.399 -0.633  2.648 ... -0.378 -0.378 -0.378]
 [-1.399 -0.633 -0.378 ... -0.378 -0.378 -0.378]
 ...
 [-1.519  1.581 -0.378 ... -0.378 -0.378  2.645]
 [-1.592  1.581 -0.378 ... -0.378 -0.378  2.645]
 [-1.592  1.581 -0.378 ... -0.378 -0.378  2.645]]
original x_test= [[2.3 0.  1.  ... 0.  0.  0. ]
 [1.5 0.  1.  ... 0.  0.  0. ]
 [1.  0.  1.  ... 0.  0.  0. ]
 ...
 [4.8 0.  0.  ... 0.  0.  1. ]
 [4.2 0.  0.  ... 0.  0.  1. ]
 [4.6 0.  0.  ... 0.  0.  1. ]]
scaled x_test= [[-1.64  -0.633  2.648 ... -0.378 -0.378 -0.378]
 [-1.737 -0.633  2.648 ... -0.378 -0.378 -0.378]
 [-1.797 -0.633  2.648 ... -0.378 -0.378 -0.378]
 ...
 [-1.339 -0.633 -0.378 ... -0.378 -0.378  2.645]
 [-1.411 -0.633 -0.378 ... -0.378 -0.378  2.645]
 [-1.363 -0.633 -0.378 ... -0.378 -0.378  2.645]]

学習用データの説明変数x_trainと検証用データの説明変数x_testに対して,平均0,標準偏差1の標準化ができていることが確認できます.

2.2.7 統計的モデルの設定と最適化

これで入力データの準備が完了しました.ここでは機械学習ライブラリsklearn(scikit-learn)の中の統計的モデルを用いて予測モデルを構築してみましょう.ここの説明では,ニューラルネットワークモデルのneural_network.MLPRegressorクラス を統計的モデルとして導入します.neural_network.MLPRegressorクラスの詳細については,アメダス気象データ分析チャレンジ!(Python版)2日目の基本編を御覧ください.

11行目のneural_network.MLPRegressor()クラスからmodelインスタンスを作成しています.ユーザーは,多数の引数を設定し,試行錯誤で ハイパーパラメータ の最適値を設定する必要があります.

model = neural_network.MLPRegressor(solver="adam", activation="relu", hidden_layer_sizes=(10,10,), max_iter=10000, tol=1e-6, random_state=1, verbose=True)

以下,ハイパーパラメータの設定について説明します.

solverでは,学習の最適化アルゴリズムを選択します.ここでの最適化アルゴリズムに基づいて予測誤差が最小になるように重み付けを決定します.ここの選択を誤ると学習速度が遅くなったり,最終的な学習結果が最適な場所(最小値)に行き着かない可能性があります. ここでは,solver=adamを設定します. solver=adamは,同じく確率的勾配降下法をベースとする手法(ADAptive Momentum estimation)で,「AdaGard法」(過去の学習によって更新されてこなかったパラメータを優先的に更新する方法)と「Momentum法」(更新量の急変を抑え滑らかに更新する方法)を組み合わせた方法で高い性能で最適化できる手法です.

activateは,活性化関数を選択します. ここでは,activate=reluを設定します. これは,ランプ関数(Rectified Linear Unit: ReLU)と呼ばれます.式で表すと,$h(x)=\max(0,x)$となります.ReLU関数の導関数はは$x>0$で常に1となるので,ニューラルネットワークの学習で長年問題となっていた「勾配消失」という問題が生じないという利点があります.ここ数年で最も人気でシンプルな活性化関数です.

hidden_layer_sizesは,中間層の数ニューロンの数を指定します. ここでは,中間層が2層で10ニューロンずつ配置することとし,hidden_layer_sizes=(10,10,)と設定します. 中間層のニューロンの数を入力層の次元よりも小さくすると次元圧縮のように機能し,主成分分析のようなことが可能になります.逆に入力層の次元よりも大きくするとスパースなデータ表現を得ることができます.中間層の層数やニューロン数を増やすことで,複雑な非線形回帰問題を解くことができます.

max_iterは,反復学習の最大回数を設定します.常に最大回数学習が回るというわけではなく,途中で学習が完了したと判断された場合はこれよりも早く終了します.ここでは,max_iter=10000と設定し,10000回を反復学習の最大回数とします.

tolは,学習の収束値を設定します.10エポック(学習回数のこと)連続でこの値よりもlossやscoreが向上しなかった場合には,学習が収束したと判断して,学習を終了します.大きくしすぎると学習途中で終了してしまいます.ここでは,tol=1e-5と設定し,lossやscoreが1e-6を下回るまで反復学習を継続します.

random_stateは,乱数発生器の番号をしてします.重み付けやバイアスの初期値やバッチサンプリング(ソルバーsgdadamのみ)で乱数を発生させる場合に,指定した番号の乱数を発生させます.これをNone(デフォルト)とすることで,毎回異なる結果となりますが,乱数発生器の番号を指定することで再現性のある同じ結果が得られます.ここでは,random_state=1と乱数を固定しています.

verboseは,進行状況のメッセージ を標準出力として出力するか出力しないかを指定します.ここでは,verbose=Trueとします.

学習用データの入力(説明変数x_train)と出力(目的変数y_train)による予測モデルの最適化計算は,下のセルの16行目の

model.fit(x_train, y_train)

により行います.そして,そのモデルによる当てはまりの良さ(決定係数$R^2$)は,学習用データ(18~19行目)と検証用データ(21~22行目)のそれぞれに対して計算されます(score_trainおよびscore_test).決定係数$R^2$が1に近いほど精度の高いモデルであると言えます.

構築された予測モデルによる予測結果は,24行目や27行目において,model.predictメソッドで計算でき,学習用データ(x_train)と検証用データ(x_test)を引数とすることで,それぞれ,

y_train_predict=model.predict(x_train)
y_test_predict=model.predict(x_test)

と計算されます.最後に,30行目や33行目において,先にpreprocessing.StandardScaler()クラスによって先にスケール変換された説明変数に対して,scaler.inverse_transform()メソッドによって変換前のオリジナルのスケールに戻すことができます.

x_train=scaler.inverse_transform([x_train])
x_test=scaler.inverse_transform([x_test])

それでは以下のセルを実行してみましょう.

In [7]:
# modelインスタンスを作成
# ニューラルネットワークモデルの場合:
#  solver="adam" 最適化手法(lbfgs, sgd, adam)
#  activateion="relu" 活性化関数(identify, logistic, tanh, relu)
#  max_iter=10000 反復の最大回数
#  tol=1e-5 学習の収束値
#  hidden_layer_sizes=(10,10,) 隠れ層のノード数(多層化可)
#  alpha=0.0001 L2正則化のペナルティー係数
#  batch_size='auto' バッチサイズの設定
#  random_state=1 重み係数やバイアスの初期値やバッチサンプリングに用いられる乱数の設定
model = neural_network.MLPRegressor(solver="adam", activation="relu", hidden_layer_sizes=(10,10,), max_iter=10000, tol=1e-5, random_state=1, verbose=True)
# 線形重回帰モデルの場合:
#model = linear_model.LinearRegression()
# fitメソッドにより目的変数y_train,説明変数x_trainによる
# 予測モデルの最適化
model.fit(x_train, y_train)
# 決定係数R2(学習データ)
score_train=model.score(x_train, y_train)
print('train R2 score=', score_train)
# 決定係数R2(検証データ)
score_test=model.score(x_test, y_test)
print('test R2 score=', score_test)
# 予測データの作成(学習データ)
y_train_predict=model.predict(x_train)
print('y_train_predict=',y_train_predict)
# 予測データの作成(検証データ)
y_test_predict=model.predict(x_test)
print('y_test_predict=',y_test_predict)
# 標準化したデータを元に戻す(学習データ)
x_train=scaler.inverse_transform([x_train])
print('original x_train=',x_train)
# 標準化したデータを元に戻す(検証データ)
x_test=scaler.inverse_transform([x_test])
print('original x_test=',x_test)
Iteration 1, loss = 5526009.39155152
Iteration 2, loss = 5521030.96747619
Iteration 3, loss = 5513918.77747958
Iteration 4, loss = 5503360.93906644
Iteration 5, loss = 5488744.91312786
Iteration 6, loss = 5468943.96139021
Iteration 7, loss = 5442011.25500002
Iteration 8, loss = 5405929.41143430
Iteration 9, loss = 5359532.71401770
Iteration 10, loss = 5302170.27276899
Iteration 11, loss = 5233092.36529831
Iteration 12, loss = 5151750.38233558
Iteration 13, loss = 5057722.97596615
Iteration 14, loss = 4950656.46586016
Iteration 15, loss = 4830359.24612216
Iteration 16, loss = 4697141.95840205
Iteration 17, loss = 4551118.73438669
Iteration 18, loss = 4392935.24560285
Iteration 19, loss = 4223366.77915567
Iteration 20, loss = 4043265.82035440
Iteration 21, loss = 3853553.71138942
Iteration 22, loss = 3655956.97768650
Iteration 23, loss = 3451733.74481843
Iteration 24, loss = 3242670.69626944
Iteration 25, loss = 3030745.68314583
Iteration 26, loss = 2817734.61542888
Iteration 27, loss = 2605063.60637487
Iteration 28, loss = 2395281.73540515
Iteration 29, loss = 2189678.48972077
Iteration 30, loss = 1990307.26379359
Iteration 31, loss = 1798498.30315708
Iteration 32, loss = 1615895.43304424
Iteration 33, loss = 1443722.35376947
Iteration 34, loss = 1282865.21203096
Iteration 35, loss = 1133869.04603910
Iteration 36, loss = 997496.44216448
Iteration 37, loss = 873799.12635523
Iteration 38, loss = 762485.47872092
Iteration 39, loss = 663750.35095017
Iteration 40, loss = 576534.29230934
Iteration 41, loss = 500571.25605225
Iteration 42, loss = 435008.91478709
Iteration 43, loss = 378693.28979490
Iteration 44, loss = 330826.35323368
Iteration 45, loss = 290377.48515395
Iteration 46, loss = 256485.77485695
Iteration 47, loss = 228214.30717182
Iteration 48, loss = 204663.30061767
Iteration 49, loss = 185129.71972775
Iteration 50, loss = 168988.04507984
Iteration 51, loss = 155562.19064881
Iteration 52, loss = 144429.01446109
Iteration 53, loss = 135136.27433365
Iteration 54, loss = 127342.24079211
Iteration 55, loss = 120835.43979339
Iteration 56, loss = 115331.11643160
Iteration 57, loss = 110644.54130416
Iteration 58, loss = 106669.99594952
Iteration 59, loss = 103270.43087533
Iteration 60, loss = 100352.31134427
Iteration 61, loss = 97837.35850474
Iteration 62, loss = 95674.43315839
Iteration 63, loss = 93798.49584594
Iteration 64, loss = 92178.53615280
Iteration 65, loss = 90763.14326763
Iteration 66, loss = 89546.02341607
Iteration 67, loss = 88485.35739127
Iteration 68, loss = 87557.26022417
Iteration 69, loss = 86739.67979828
Iteration 70, loss = 86040.14034749
Iteration 71, loss = 85405.53232848
Iteration 72, loss = 84851.87165165
Iteration 73, loss = 84367.70329660
Iteration 74, loss = 83921.10732234
Iteration 75, loss = 83534.65406575
Iteration 76, loss = 83187.31404234
Iteration 77, loss = 82862.90752547
Iteration 78, loss = 82571.77679455
Iteration 79, loss = 82295.76370597
Iteration 80, loss = 82055.05119918
Iteration 81, loss = 81819.90244336
Iteration 82, loss = 81605.90693108
Iteration 83, loss = 81398.66234563
Iteration 84, loss = 81198.73277596
Iteration 85, loss = 81003.30178102
Iteration 86, loss = 80815.78273902
Iteration 87, loss = 80643.93141851
Iteration 88, loss = 80470.55025254
Iteration 89, loss = 80303.92621614
Iteration 90, loss = 80139.62252890
Iteration 91, loss = 79978.02368561
Iteration 92, loss = 79820.68006669
Iteration 93, loss = 79671.09965625
Iteration 94, loss = 79514.92702785
Iteration 95, loss = 79367.85794928
Iteration 96, loss = 79210.99698420
Iteration 97, loss = 79055.79500687
Iteration 98, loss = 78905.85826320
Iteration 99, loss = 78738.66281846
Iteration 100, loss = 78567.49098945
Iteration 101, loss = 78408.10486622
Iteration 102, loss = 78243.10149558
Iteration 103, loss = 78066.35698586
Iteration 104, loss = 77924.96583569
Iteration 105, loss = 77740.37925394
Iteration 106, loss = 77584.29423179
Iteration 107, loss = 77429.43197369
Iteration 108, loss = 77262.78204945
Iteration 109, loss = 77109.23431246
Iteration 110, loss = 76953.93605286
Iteration 111, loss = 76795.33066611
Iteration 112, loss = 76651.73606408
Iteration 113, loss = 76508.46389127
Iteration 114, loss = 76369.26744072
Iteration 115, loss = 76236.55900612
Iteration 116, loss = 76102.08334963
Iteration 117, loss = 75971.36976012
Iteration 118, loss = 75855.57136230
Iteration 119, loss = 75729.82227621
Iteration 120, loss = 75625.84830808
Iteration 121, loss = 75474.02620521
Iteration 122, loss = 75344.76903397
Iteration 123, loss = 75227.99107586
Iteration 124, loss = 75096.20836315
Iteration 125, loss = 74962.92532583
Iteration 126, loss = 74835.01570336
Iteration 127, loss = 74705.34725505
Iteration 128, loss = 74580.48105994
Iteration 129, loss = 74451.01372253
Iteration 130, loss = 74331.47075704
Iteration 131, loss = 74204.75322531
Iteration 132, loss = 74098.77761793
Iteration 133, loss = 73983.97372114
Iteration 134, loss = 73859.77581502
Iteration 135, loss = 73758.51377042
Iteration 136, loss = 73659.32591739
Iteration 137, loss = 73545.15006385
Iteration 138, loss = 73471.63900830
Iteration 139, loss = 73337.44100377
Iteration 140, loss = 73231.82276117
Iteration 141, loss = 73133.99406410
Iteration 142, loss = 73022.61669153
Iteration 143, loss = 72936.92592504
Iteration 144, loss = 72834.07431306
Iteration 145, loss = 72735.76123111
Iteration 146, loss = 72636.48071571
Iteration 147, loss = 72517.56404604
Iteration 148, loss = 72433.67601485
Iteration 149, loss = 72337.99619461
Iteration 150, loss = 72238.84764281
Iteration 151, loss = 72160.70160511
Iteration 152, loss = 72066.64892797
Iteration 153, loss = 71992.44241407
Iteration 154, loss = 71887.18994627
Iteration 155, loss = 71807.36876855
Iteration 156, loss = 71699.43049341
Iteration 157, loss = 71632.91312805
Iteration 158, loss = 71546.25886137
Iteration 159, loss = 71462.25368900
Iteration 160, loss = 71367.39732927
Iteration 161, loss = 71273.97756735
Iteration 162, loss = 71193.50037361
Iteration 163, loss = 71082.04329324
Iteration 164, loss = 70992.15180136
Iteration 165, loss = 70906.88475990
Iteration 166, loss = 70802.88541209
Iteration 167, loss = 70702.16623198
Iteration 168, loss = 70595.39340006
Iteration 169, loss = 70489.68362005
Iteration 170, loss = 70394.54828364
Iteration 171, loss = 70302.42069323
Iteration 172, loss = 70174.02194652
Iteration 173, loss = 70098.82919836
Iteration 174, loss = 69982.38719662
Iteration 175, loss = 69883.17642068
Iteration 176, loss = 69807.11231402
Iteration 177, loss = 69731.35728357
Iteration 178, loss = 69621.08784036
Iteration 179, loss = 69558.14571034
Iteration 180, loss = 69432.76103948
Iteration 181, loss = 69352.99873004
Iteration 182, loss = 69269.98065073
Iteration 183, loss = 69179.21510569
Iteration 184, loss = 69081.16655269
Iteration 185, loss = 68989.19362428
Iteration 186, loss = 68886.31489542
Iteration 187, loss = 68783.99693936
Iteration 188, loss = 68674.03463633
Iteration 189, loss = 68578.40936636
Iteration 190, loss = 68470.96478004
Iteration 191, loss = 68362.62656792
Iteration 192, loss = 68208.87183829
Iteration 193, loss = 68066.49211085
Iteration 194, loss = 67892.56851577
Iteration 195, loss = 67725.96677211
Iteration 196, loss = 67556.43600465
Iteration 197, loss = 67336.02851186
Iteration 198, loss = 67147.09983360
Iteration 199, loss = 66942.97297238
Iteration 200, loss = 66736.98241246
Iteration 201, loss = 66566.14379184
Iteration 202, loss = 66406.04489834
Iteration 203, loss = 66201.34643076
Iteration 204, loss = 66024.38015032
Iteration 205, loss = 65864.38634688
Iteration 206, loss = 65674.31266108
Iteration 207, loss = 65480.31169217
Iteration 208, loss = 65298.26380797
Iteration 209, loss = 65101.08432525
Iteration 210, loss = 64745.57705548
Iteration 211, loss = 64490.94574600
Iteration 212, loss = 64218.63943805
Iteration 213, loss = 63987.11486465
Iteration 214, loss = 63752.45348520
Iteration 215, loss = 63549.22773275
Iteration 216, loss = 63301.06543883
Iteration 217, loss = 63061.75275380
Iteration 218, loss = 62846.88726884
Iteration 219, loss = 62599.90823028
Iteration 220, loss = 62338.09004689
Iteration 221, loss = 62063.34692032
Iteration 222, loss = 61804.86268407
Iteration 223, loss = 61560.08931889
Iteration 224, loss = 61311.51770224
Iteration 225, loss = 61063.89909210
Iteration 226, loss = 60820.89053198
Iteration 227, loss = 60545.20468957
Iteration 228, loss = 60318.54870925
Iteration 229, loss = 60066.48306553
Iteration 230, loss = 59780.91445733
Iteration 231, loss = 59529.60090102
Iteration 232, loss = 59265.00204052
Iteration 233, loss = 58959.58170429
Iteration 234, loss = 58679.19893954
Iteration 235, loss = 58378.81982631
Iteration 236, loss = 58077.57739621
Iteration 237, loss = 57777.65256198
Iteration 238, loss = 57476.65028310
Iteration 239, loss = 57178.02014488
Iteration 240, loss = 56866.76890916
Iteration 241, loss = 56530.54722118
Iteration 242, loss = 56196.04955146
Iteration 243, loss = 55850.31893267
Iteration 244, loss = 55491.70464333
Iteration 245, loss = 55041.90814633
Iteration 246, loss = 54600.15839112
Iteration 247, loss = 54176.72783143
Iteration 248, loss = 53800.75324281
Iteration 249, loss = 53445.13063151
Iteration 250, loss = 53104.40777543
Iteration 251, loss = 52717.61987868
Iteration 252, loss = 52372.07594424
Iteration 253, loss = 51997.55868817
Iteration 254, loss = 51615.81684575
Iteration 255, loss = 51238.98108781
Iteration 256, loss = 50820.24023888
Iteration 257, loss = 50388.94802731
Iteration 258, loss = 49978.70967263
Iteration 259, loss = 49527.14485004
Iteration 260, loss = 49119.00228441
Iteration 261, loss = 48743.86260742
Iteration 262, loss = 48352.41391749
Iteration 263, loss = 47970.41742218
Iteration 264, loss = 47607.15891848
Iteration 265, loss = 47211.77044119
Iteration 266, loss = 46838.00589840
Iteration 267, loss = 46449.57011788
Iteration 268, loss = 46029.97342753
Iteration 269, loss = 45633.66129242
Iteration 270, loss = 45184.98149400
Iteration 271, loss = 44734.34048843
Iteration 272, loss = 44269.49975708
Iteration 273, loss = 43770.51067911
Iteration 274, loss = 43290.66470321
Iteration 275, loss = 42827.30525583
Iteration 276, loss = 42368.04750025
Iteration 277, loss = 41929.42640964
Iteration 278, loss = 41477.24494143
Iteration 279, loss = 41033.33329225
Iteration 280, loss = 40604.62902712
Iteration 281, loss = 40174.58453113
Iteration 282, loss = 39732.87110891
Iteration 283, loss = 39308.44395813
Iteration 284, loss = 38896.54456412
Iteration 285, loss = 38462.44511707
Iteration 286, loss = 38066.03745752
Iteration 287, loss = 37655.54971989
Iteration 288, loss = 37249.23818971
Iteration 289, loss = 36882.79728508
Iteration 290, loss = 36489.71920877
Iteration 291, loss = 36136.93924635
Iteration 292, loss = 35768.89675438
Iteration 293, loss = 35412.70296576
Iteration 294, loss = 35056.08060907
Iteration 295, loss = 34726.25712217
Iteration 296, loss = 34369.05837175
Iteration 297, loss = 34054.56356511
Iteration 298, loss = 33733.48503423
Iteration 299, loss = 33391.05819128
Iteration 300, loss = 33084.91254598
Iteration 301, loss = 32767.48390717
Iteration 302, loss = 32468.72026029
Iteration 303, loss = 32177.67610720
Iteration 304, loss = 31900.77924441
Iteration 305, loss = 31626.75093865
Iteration 306, loss = 31381.54601184
Iteration 307, loss = 31121.07272360
Iteration 308, loss = 30884.89468103
Iteration 309, loss = 30658.96270010
Iteration 310, loss = 30424.47879317
Iteration 311, loss = 30226.85517858
Iteration 312, loss = 30030.73877701
Iteration 313, loss = 29855.67645187
Iteration 314, loss = 29696.39234624
Iteration 315, loss = 29486.58164886
Iteration 316, loss = 29337.90498082
Iteration 317, loss = 29183.59240350
Iteration 318, loss = 29029.75823048
Iteration 319, loss = 28885.48264390
Iteration 320, loss = 28751.86258107
Iteration 321, loss = 28620.09214594
Iteration 322, loss = 28502.50990547
Iteration 323, loss = 28369.73689315
Iteration 324, loss = 28253.71549955
Iteration 325, loss = 28135.97071832
Iteration 326, loss = 28039.60584656
Iteration 327, loss = 27920.38649898
Iteration 328, loss = 27832.09995599
Iteration 329, loss = 27737.41175254
Iteration 330, loss = 27641.74962103
Iteration 331, loss = 27562.76678803
Iteration 332, loss = 27489.04634618
Iteration 333, loss = 27407.05115226
Iteration 334, loss = 27342.53918061
Iteration 335, loss = 27259.11326930
Iteration 336, loss = 27199.12011933
Iteration 337, loss = 27134.94902631
Iteration 338, loss = 27078.65704583
Iteration 339, loss = 27018.83986563
Iteration 340, loss = 26972.98383664
Iteration 341, loss = 26929.10994873
Iteration 342, loss = 26865.13074592
Iteration 343, loss = 26814.71563995
Iteration 344, loss = 26768.63666822
Iteration 345, loss = 26717.98275442
Iteration 346, loss = 26699.46345256
Iteration 347, loss = 26662.21408424
Iteration 348, loss = 26602.15093370
Iteration 349, loss = 26564.59670345
Iteration 350, loss = 26521.50942203
Iteration 351, loss = 26490.28861309
Iteration 352, loss = 26465.85210562
Iteration 353, loss = 26419.63604835
Iteration 354, loss = 26397.14221173
Iteration 355, loss = 26392.06786651
Iteration 356, loss = 26331.84442870
Iteration 357, loss = 26318.54172632
Iteration 358, loss = 26267.45302545
Iteration 359, loss = 26266.35461193
Iteration 360, loss = 26226.00726800
Iteration 361, loss = 26191.86339508
Iteration 362, loss = 26176.57710334
Iteration 363, loss = 26149.28731990
Iteration 364, loss = 26130.76902186
Iteration 365, loss = 26096.02338532
Iteration 366, loss = 26085.99038187
Iteration 367, loss = 26066.65490711
Iteration 368, loss = 26037.50587238
Iteration 369, loss = 26028.78521877
Iteration 370, loss = 26002.77293252
Iteration 371, loss = 25977.37485791
Iteration 372, loss = 25962.61688409
Iteration 373, loss = 25940.93165184
Iteration 374, loss = 25917.74933655
Iteration 375, loss = 25904.81747025
Iteration 376, loss = 25882.69930130
Iteration 377, loss = 25868.16660695
Iteration 378, loss = 25842.50667258
Iteration 379, loss = 25845.66676443
Iteration 380, loss = 25818.71922182
Iteration 381, loss = 25798.66254246
Iteration 382, loss = 25790.86599401
Iteration 383, loss = 25756.91365030
Iteration 384, loss = 25755.51412764
Iteration 385, loss = 25734.54628524
Iteration 386, loss = 25740.01334573
Iteration 387, loss = 25720.41558809
Iteration 388, loss = 25683.01961205
Iteration 389, loss = 25680.03983774
Iteration 390, loss = 25667.15364364
Iteration 391, loss = 25645.77064213
Iteration 392, loss = 25661.42362678
Iteration 393, loss = 25640.67770310
Iteration 394, loss = 25601.98045036
Iteration 395, loss = 25596.07622044
Iteration 396, loss = 25593.71028938
Iteration 397, loss = 25578.07502197
Iteration 398, loss = 25563.82414967
Iteration 399, loss = 25559.47659047
Iteration 400, loss = 25538.07510173
Iteration 401, loss = 25524.67474584
Iteration 402, loss = 25519.46924675
Iteration 403, loss = 25503.42268022
Iteration 404, loss = 25487.31754604
Iteration 405, loss = 25472.62101640
Iteration 406, loss = 25473.31918570
Iteration 407, loss = 25459.16920847
Iteration 408, loss = 25446.15278809
Iteration 409, loss = 25429.43440504
Iteration 410, loss = 25414.64892088
Iteration 411, loss = 25431.33535986
Iteration 412, loss = 25401.94913677
Iteration 413, loss = 25402.70891743
Iteration 414, loss = 25385.48863166
Iteration 415, loss = 25367.81896775
Iteration 416, loss = 25351.35838476
Iteration 417, loss = 25354.84920264
Iteration 418, loss = 25338.91774560
Iteration 419, loss = 25334.53316854
Iteration 420, loss = 25316.09051339
Iteration 421, loss = 25303.54634259
Iteration 422, loss = 25303.46992224
Iteration 423, loss = 25279.98376709
Iteration 424, loss = 25264.31836658
Iteration 425, loss = 25251.63284679
Iteration 426, loss = 25237.59909099
Iteration 427, loss = 25225.91955151
Iteration 428, loss = 25216.50518557
Iteration 429, loss = 25197.26884166
Iteration 430, loss = 25189.74248528
Iteration 431, loss = 25181.73925916
Iteration 432, loss = 25178.98617197
Iteration 433, loss = 25176.10710212
Iteration 434, loss = 25150.94099393
Iteration 435, loss = 25123.08497890
Iteration 436, loss = 25129.81778763
Iteration 437, loss = 25124.63889637
Iteration 438, loss = 25129.48608518
Iteration 439, loss = 25114.02971114
Iteration 440, loss = 25086.57017604
Iteration 441, loss = 25077.33303426
Iteration 442, loss = 25071.94469670
Iteration 443, loss = 25060.04967944
Iteration 444, loss = 25055.46784132
Iteration 445, loss = 25045.15503869
Iteration 446, loss = 25043.74906239
Iteration 447, loss = 25031.71248447
Iteration 448, loss = 25028.27509140
Iteration 449, loss = 25014.36974658
Iteration 450, loss = 25006.23451627
Iteration 451, loss = 24989.56960986
Iteration 452, loss = 25015.84168981
Iteration 453, loss = 24988.08912229
Iteration 454, loss = 24986.32356984
Iteration 455, loss = 24975.10418324
Iteration 456, loss = 24966.72130975
Iteration 457, loss = 24961.38855562
Iteration 458, loss = 24959.10451957
Iteration 459, loss = 24937.25482200
Iteration 460, loss = 24947.49504337
Iteration 461, loss = 24945.79704600
Iteration 462, loss = 24928.86393379
Iteration 463, loss = 24927.54894681
Iteration 464, loss = 24915.44167768
Iteration 465, loss = 24916.93931838
Iteration 466, loss = 24905.38599109
Iteration 467, loss = 24901.21679312
Iteration 468, loss = 24897.43867564
Iteration 469, loss = 24893.52846541
Iteration 470, loss = 24900.54900545
Iteration 471, loss = 24891.14944947
Iteration 472, loss = 24874.73750001
Iteration 473, loss = 24882.82225265
Iteration 474, loss = 24872.07062626
Iteration 475, loss = 24868.08914011
Iteration 476, loss = 24873.16343363
Iteration 477, loss = 24842.40875639
Iteration 478, loss = 24854.14042381
Iteration 479, loss = 24859.49559406
Iteration 480, loss = 24848.73942774
Iteration 481, loss = 24861.12970665
Iteration 482, loss = 24829.60236512
Iteration 483, loss = 24853.15693442
Iteration 484, loss = 24815.31966909
Iteration 485, loss = 24824.58559652
Iteration 486, loss = 24808.97685075
Iteration 487, loss = 24814.14725277
Iteration 488, loss = 24796.15080471
Iteration 489, loss = 24803.03048936
Iteration 490, loss = 24804.67794927
Iteration 491, loss = 24781.35605186
Iteration 492, loss = 24795.50725503
Iteration 493, loss = 24790.34498380
Iteration 494, loss = 24786.90963792
Iteration 495, loss = 24775.39058185
Iteration 496, loss = 24775.65465739
Iteration 497, loss = 24778.70834585
Iteration 498, loss = 24758.17982763
Iteration 499, loss = 24772.46157485
Iteration 500, loss = 24751.48003722
Iteration 501, loss = 24752.42236991
Iteration 502, loss = 24752.61301915
Iteration 503, loss = 24752.93303694
Iteration 504, loss = 24743.09402400
Iteration 505, loss = 24764.52697005
Iteration 506, loss = 24730.77301245
Iteration 507, loss = 24725.53149400
Iteration 508, loss = 24735.45082286
Iteration 509, loss = 24731.23478605
Iteration 510, loss = 24715.09250999
Iteration 511, loss = 24717.44109231
Iteration 512, loss = 24707.86154945
Iteration 513, loss = 24704.48013083
Iteration 514, loss = 24709.39381564
Iteration 515, loss = 24699.82159574
Iteration 516, loss = 24705.55646992
Iteration 517, loss = 24707.71652103
Iteration 518, loss = 24692.82003141
Iteration 519, loss = 24690.32449948
Iteration 520, loss = 24672.52774041
Iteration 521, loss = 24682.73395859
Iteration 522, loss = 24674.74916391
Iteration 523, loss = 24662.97515983
Iteration 524, loss = 24663.16586903
Iteration 525, loss = 24663.56214148
Iteration 526, loss = 24649.37892532
Iteration 527, loss = 24654.17531685
Iteration 528, loss = 24654.15170352
Iteration 529, loss = 24652.73394863
Iteration 530, loss = 24640.33704742
Iteration 531, loss = 24633.45377172
Iteration 532, loss = 24649.04492721
Iteration 533, loss = 24625.06783046
Iteration 534, loss = 24625.73910653
Iteration 535, loss = 24624.37588082
Iteration 536, loss = 24619.14192082
Iteration 537, loss = 24620.26362766
Iteration 538, loss = 24632.53916502
Iteration 539, loss = 24605.96135948
Iteration 540, loss = 24614.59720841
Iteration 541, loss = 24595.04854348
Iteration 542, loss = 24603.28966282
Iteration 543, loss = 24600.87222149
Iteration 544, loss = 24596.44775251
Iteration 545, loss = 24578.43876174
Iteration 546, loss = 24588.38644173
Iteration 547, loss = 24579.12147300
Iteration 548, loss = 24572.77307521
Iteration 549, loss = 24578.42428694
Iteration 550, loss = 24573.86779952
Iteration 551, loss = 24564.07683860
Iteration 552, loss = 24564.29560692
Iteration 553, loss = 24567.60312503
Iteration 554, loss = 24568.45786773
Iteration 555, loss = 24576.67743953
Iteration 556, loss = 24548.85681108
Iteration 557, loss = 24552.85180165
Iteration 558, loss = 24544.49539712
Iteration 559, loss = 24544.72591151
Iteration 560, loss = 24537.94094543
Iteration 561, loss = 24542.97820206
Iteration 562, loss = 24539.98081315
Iteration 563, loss = 24538.48026394
Iteration 564, loss = 24528.49074599
Iteration 565, loss = 24529.25375825
Iteration 566, loss = 24522.77351568
Iteration 567, loss = 24520.70931774
Iteration 568, loss = 24529.59991830
Iteration 569, loss = 24513.06054038
Iteration 570, loss = 24519.65812405
Iteration 571, loss = 24517.33141606
Iteration 572, loss = 24514.10723073
Iteration 573, loss = 24495.04849222
Iteration 574, loss = 24501.43615701
Iteration 575, loss = 24508.39513098
Iteration 576, loss = 24487.47157624
Iteration 577, loss = 24486.84135226
Iteration 578, loss = 24482.00469192
Iteration 579, loss = 24487.23524825
Iteration 580, loss = 24509.68334100
Iteration 581, loss = 24482.09610864
Iteration 582, loss = 24474.17824055
Iteration 583, loss = 24500.01283933
Iteration 584, loss = 24474.41134049
Iteration 585, loss = 24464.21506712
Iteration 586, loss = 24469.75139670
Iteration 587, loss = 24477.13223399
Iteration 588, loss = 24487.74369932
Iteration 589, loss = 24465.75858094
Iteration 590, loss = 24458.39779826
Iteration 591, loss = 24454.45965897
Iteration 592, loss = 24456.58712460
Iteration 593, loss = 24453.28449678
Iteration 594, loss = 24444.53968347
Iteration 595, loss = 24449.33466268
Iteration 596, loss = 24453.03888270
Iteration 597, loss = 24446.06330575
Iteration 598, loss = 24434.34711606
Iteration 599, loss = 24430.48489020
Iteration 600, loss = 24426.29723571
Iteration 601, loss = 24422.85489651
Iteration 602, loss = 24420.33605582
Iteration 603, loss = 24421.81965737
Iteration 604, loss = 24418.11297801
Iteration 605, loss = 24439.74662482
Iteration 606, loss = 24411.94411690
Iteration 607, loss = 24443.89850876
Iteration 608, loss = 24418.42686798
Iteration 609, loss = 24405.49024546
Iteration 610, loss = 24397.61974784
Iteration 611, loss = 24416.15072869
Iteration 612, loss = 24417.09344003
Iteration 613, loss = 24425.47498458
Iteration 614, loss = 24399.30361404
Iteration 615, loss = 24395.36077099
Iteration 616, loss = 24398.26564212
Iteration 617, loss = 24382.72832986
Iteration 618, loss = 24383.62454332
Iteration 619, loss = 24382.09088382
Iteration 620, loss = 24377.52407640
Iteration 621, loss = 24400.70257173
Iteration 622, loss = 24385.50494274
Iteration 623, loss = 24383.26382535
Iteration 624, loss = 24364.15551974
Iteration 625, loss = 24385.87933093
Iteration 626, loss = 24375.94471865
Iteration 627, loss = 24371.92008165
Iteration 628, loss = 24362.47154918
Iteration 629, loss = 24373.15092142
Iteration 630, loss = 24350.82401263
Iteration 631, loss = 24356.44013687
Iteration 632, loss = 24353.96799793
Iteration 633, loss = 24346.15432320
Iteration 634, loss = 24343.48402814
Iteration 635, loss = 24352.55854700
Iteration 636, loss = 24339.61107613
Iteration 637, loss = 24333.49717213
Iteration 638, loss = 24345.17744163
Iteration 639, loss = 24338.52944259
Iteration 640, loss = 24324.20337676
Iteration 641, loss = 24332.80841031
Iteration 642, loss = 24320.99018177
Iteration 643, loss = 24328.22755194
Iteration 644, loss = 24338.19856597
Iteration 645, loss = 24333.38752618
Iteration 646, loss = 24324.95150420
Iteration 647, loss = 24313.44043688
Iteration 648, loss = 24309.68342451
Iteration 649, loss = 24315.05986216
Iteration 650, loss = 24309.46841514
Iteration 651, loss = 24300.43435648
Iteration 652, loss = 24305.45001406
Iteration 653, loss = 24294.02444061
Iteration 654, loss = 24296.50402190
Iteration 655, loss = 24297.91846515
Iteration 656, loss = 24299.25117126
Iteration 657, loss = 24290.31122614
Iteration 658, loss = 24293.36795048
Iteration 659, loss = 24294.18097784
Iteration 660, loss = 24278.85566510
Iteration 661, loss = 24292.50939055
Iteration 662, loss = 24281.65565269
Iteration 663, loss = 24286.13111059
Iteration 664, loss = 24285.49793739
Iteration 665, loss = 24272.29662124
Iteration 666, loss = 24273.27856016
Iteration 667, loss = 24281.50006015
Iteration 668, loss = 24262.93010277
Iteration 669, loss = 24255.49501669
Iteration 670, loss = 24258.45381553
Iteration 671, loss = 24263.38034862
Iteration 672, loss = 24277.20412829
Iteration 673, loss = 24266.47708687
Iteration 674, loss = 24255.08152905
Iteration 675, loss = 24249.88367857
Iteration 676, loss = 24258.80067631
Iteration 677, loss = 24260.60577685
Iteration 678, loss = 24253.12733942
Iteration 679, loss = 24249.22479266
Iteration 680, loss = 24250.75221294
Iteration 681, loss = 24268.67519172
Iteration 682, loss = 24242.30102274
Iteration 683, loss = 24232.74709564
Iteration 684, loss = 24234.30135894
Iteration 685, loss = 24236.10228798
Iteration 686, loss = 24241.45493106
Iteration 687, loss = 24234.76403654
Iteration 688, loss = 24237.45876326
Iteration 689, loss = 24234.69116786
Iteration 690, loss = 24237.42517739
Iteration 691, loss = 24248.45180880
Iteration 692, loss = 24217.78042607
Iteration 693, loss = 24236.61961957
Iteration 694, loss = 24221.89026044
Iteration 695, loss = 24214.00370156
Iteration 696, loss = 24214.76407566
Iteration 697, loss = 24214.70463534
Iteration 698, loss = 24221.84075692
Iteration 699, loss = 24214.34352557
Iteration 700, loss = 24214.78474311
Iteration 701, loss = 24209.94189639
Iteration 702, loss = 24220.66468754
Iteration 703, loss = 24224.35258935
Iteration 704, loss = 24217.05477129
Iteration 705, loss = 24215.59769598
Iteration 706, loss = 24193.88018822
Iteration 707, loss = 24211.33374235
Iteration 708, loss = 24210.05493304
Iteration 709, loss = 24193.21936380
Iteration 710, loss = 24200.56364884
Iteration 711, loss = 24187.86419332
Iteration 712, loss = 24186.58315272
Iteration 713, loss = 24192.30441078
Iteration 714, loss = 24185.53277905
Iteration 715, loss = 24181.82707390
Iteration 716, loss = 24187.50960939
Iteration 717, loss = 24196.40255890
Iteration 718, loss = 24192.04349786
Iteration 719, loss = 24192.93738818
Iteration 720, loss = 24204.55877067
Iteration 721, loss = 24171.89514635
Iteration 722, loss = 24173.12585907
Iteration 723, loss = 24166.35294916
Iteration 724, loss = 24159.67719457
Iteration 725, loss = 24178.33259962
Iteration 726, loss = 24164.01868684
Iteration 727, loss = 24189.24462797
Iteration 728, loss = 24164.93386175
Iteration 729, loss = 24159.17418991
Iteration 730, loss = 24162.83207300
Iteration 731, loss = 24163.59464797
Iteration 732, loss = 24177.80823364
Iteration 733, loss = 24172.29781918
Iteration 734, loss = 24160.95361360
Iteration 735, loss = 24162.81582655
Iteration 736, loss = 24159.08744949
Iteration 737, loss = 24168.42307540
Iteration 738, loss = 24153.63506594
Iteration 739, loss = 24170.14393718
Iteration 740, loss = 24143.07013953
Iteration 741, loss = 24150.20109143
Iteration 742, loss = 24152.55860206
Iteration 743, loss = 24145.73132956
Iteration 744, loss = 24159.41672939
Iteration 745, loss = 24148.02461678
Iteration 746, loss = 24144.03961321
Iteration 747, loss = 24142.01925275
Iteration 748, loss = 24153.03574489
Iteration 749, loss = 24137.46738280
Iteration 750, loss = 24147.57088591
Iteration 751, loss = 24160.08666744
Iteration 752, loss = 24138.90480243
Iteration 753, loss = 24137.34157152
Iteration 754, loss = 24137.79135628
Iteration 755, loss = 24141.34413626
Iteration 756, loss = 24133.88737238
Iteration 757, loss = 24127.14746650
Iteration 758, loss = 24144.46314762
Iteration 759, loss = 24133.84480961
Iteration 760, loss = 24135.87466267
Iteration 761, loss = 24133.99142754
Iteration 762, loss = 24130.46181088
Iteration 763, loss = 24120.50462358
Iteration 764, loss = 24130.13009435
Iteration 765, loss = 24116.84375611
Iteration 766, loss = 24130.78600442
Iteration 767, loss = 24146.47162761
Iteration 768, loss = 24123.44461690
Iteration 769, loss = 24123.55694354
Iteration 770, loss = 24134.25181701
Iteration 771, loss = 24123.72620876
Iteration 772, loss = 24114.33596837
Iteration 773, loss = 24129.76240435
Iteration 774, loss = 24140.08622093
Iteration 775, loss = 24127.96246910
Iteration 776, loss = 24110.13843567
Iteration 777, loss = 24115.25072331
Iteration 778, loss = 24126.96528499
Iteration 779, loss = 24114.88149357
Iteration 780, loss = 24109.56383703
Iteration 781, loss = 24103.81933425
Iteration 782, loss = 24116.39125317
Iteration 783, loss = 24121.58377503
Iteration 784, loss = 24107.81965645
Iteration 785, loss = 24114.25928480
Iteration 786, loss = 24117.63501283
Iteration 787, loss = 24105.67808069
Iteration 788, loss = 24105.53329327
Iteration 789, loss = 24106.24584247
Iteration 790, loss = 24100.91344471
Iteration 791, loss = 24108.79807854
Iteration 792, loss = 24118.75615034
Iteration 793, loss = 24096.65494942
Iteration 794, loss = 24107.77719369
Iteration 795, loss = 24103.77591660
Iteration 796, loss = 24108.32438056
Iteration 797, loss = 24122.27252679
Iteration 798, loss = 24098.36596943
Iteration 799, loss = 24095.93209765
Iteration 800, loss = 24104.72964076
Iteration 801, loss = 24102.23677993
Iteration 802, loss = 24094.32903566
Iteration 803, loss = 24099.07211507
Iteration 804, loss = 24104.75197548
Iteration 805, loss = 24100.55388481
Iteration 806, loss = 24103.91752554
Iteration 807, loss = 24089.91038593
Iteration 808, loss = 24101.01603622
Iteration 809, loss = 24088.26480056
Iteration 810, loss = 24102.80534088
Iteration 811, loss = 24084.70908784
Iteration 812, loss = 24092.58352411
Iteration 813, loss = 24094.22036347
Iteration 814, loss = 24100.42401689
Iteration 815, loss = 24099.83269656
Iteration 816, loss = 24083.10801264
Iteration 817, loss = 24089.10580237
Iteration 818, loss = 24096.19923371
Iteration 819, loss = 24080.97851890
Iteration 820, loss = 24096.30158764
Iteration 821, loss = 24076.61520751
Iteration 822, loss = 24094.22490309
Iteration 823, loss = 24081.43879402
Iteration 824, loss = 24094.30087390
Iteration 825, loss = 24075.56764537
Iteration 826, loss = 24089.02282661
Iteration 827, loss = 24076.29358873
Iteration 828, loss = 24089.65735564
Iteration 829, loss = 24081.91563205
Iteration 830, loss = 24099.90305851
Iteration 831, loss = 24069.82654357
Iteration 832, loss = 24075.27474529
Iteration 833, loss = 24081.09808310
Iteration 834, loss = 24072.49909765
Iteration 835, loss = 24065.17063459
Iteration 836, loss = 24082.75877311
Iteration 837, loss = 24071.94559666
Iteration 838, loss = 24081.62551839
Iteration 839, loss = 24069.98628603
Iteration 840, loss = 24067.00754319
Iteration 841, loss = 24064.65258967
Iteration 842, loss = 24079.21989564
Iteration 843, loss = 24064.05964430
Iteration 844, loss = 24069.95438478
Iteration 845, loss = 24077.45275966
Iteration 846, loss = 24050.77918750
Iteration 847, loss = 24078.54426970
Iteration 848, loss = 24076.85277752
Iteration 849, loss = 24076.51088421
Iteration 850, loss = 24064.18582990
Iteration 851, loss = 24075.42785565
Iteration 852, loss = 24063.17166332
Iteration 853, loss = 24060.48874850
Iteration 854, loss = 24059.49998649
Iteration 855, loss = 24059.33152804
Iteration 856, loss = 24054.28714868
Iteration 857, loss = 24063.43881214
Training loss did not improve more than tol=0.000010 for 10 consecutive epochs. Stopping.
train R2 score= 0.8757665807741901
test R2 score= 0.8665886639074808
y_train_predict= [2960.689 2995.327 2795.9   ... 3594.133 3646.591 3646.591]
y_test_predict= [3133.878 3189.298 3223.936 ... 3739.309 3791.767 3756.795]
original x_train= [[[4.8 0.  1.  ... 0.  0.  0. ]
  [4.3 0.  1.  ... 0.  0.  0. ]
  [4.3 0.  0.  ... 0.  0.  0. ]
  ...
  [3.3 1.  0.  ... 0.  0.  1. ]
  [2.7 1.  0.  ... 0.  0.  1. ]
  [2.7 1.  0.  ... 0.  0.  1. ]]]
original x_test= [[[2.3 0.  1.  ... 0.  0.  0. ]
  [1.5 0.  1.  ... 0.  0.  0. ]
  [1.  0.  1.  ... 0.  0.  0. ]
  ...
  [4.8 0.  0.  ... 0.  0.  1. ]
  [4.2 0.  0.  ... 0.  0.  1. ]
  [4.6 0.  0.  ... 0.  0.  1. ]]]

さて,結果を見てみましょう.ニューラルネットワークモデルの予測精度(決定係数$R^2$)は,学習用データは0.88を,検証用データは0.87となりました.気温,曜日,時間を説明変数とすることで,精度良く電力消費量を予測できると言えるでしょう.

さて,2.2.5節に戻って,1)説明変数の種類を減らしたり増やしたりすると予測精度がどう変化するのか調べてみましょう.予測的分析では,試行錯誤的に説明変数の組み合わせを変化させることで,モデルの精度を比較しながら因果関係を分析していきます.

また,ニューラルネットワークモデルの2)ニューロンの数を増やしたら予測精度がどう変化するか,3)中間層の数を増やしたら予測精度がどう変化するか調べてみましょう.

さらには,ニューラルネットワークモデル以外の線形重回帰モデルLinearRegression()を使ったら予測精度はどう変化するか調べてみましょう.線形重回帰モデルとは,$y=a x_1 + b x_2 + c x_3 + \cdots + d$のように複数の説明変数($x_1$,$x_2$,$x_3$,…)により目的変数($y$)を予測する統計的モデルです.重回帰モデルでは上手く行かず,ニューラルネットワークでは上手くいく理由を考察してみましょう.

そして,電力消費量予測の精度をどこまで向上させることができるでしょうか?皆さんも,自分自身で新たな説明変数を追加するなどして,決定係数$R^2>$0.99の予測モデル開発に挑戦してみましょう.

統計的モデルの最適化計算をやり直す場合には,2.2.5節のセルまで戻る必要がありますので気をつけて下さい.

2.2.8 時系列図と散布図の作成

次に,時系列図を作成して調べてみましょう.予測データと正解データを適切な形にグラフ化することで,その特徴がみえてきます.数値だけでは分かりにくい例外的な特徴も,グラフにすることによって抽出することできる.ここでは,pythonのmatplotlibを使って時系列図を作成します.表計算ソフトExcelで例えるならば,3つの列データを選択して,ツールバーの中の折れ線グラフボタンを押すプロセスです.

次のセルは,時系列図を作成する関数timeseries(2行目)を定義し,それを呼び出すことで,訓練データの中の正解データy_trainと予測データy_train_predictを縦軸に,データフレームdataのインデックスから抽出したdatetime_xを横軸に設定して作図します(36行目).datetime_xのスライスの中のdt1は,全体の行(len(data))を2で割った値が設定されています(29行目).関数timeseriesの中では,まず,横軸のデータX(4行目)と縦軸の正解データY1(6行目)と予測データY2(8行目)を割り当てています.2つの時系列を1枚の図に収めます.図のサイズを10行目で設定しています.12行目と14行目で正解データY1(青線)と予測データY2(赤線)の時系列図を描画し,15行目~20行目で図のタイトルや軸ラベルの設定をし,22行目で図の凡例を設定しています.

同じ要領で関数timeseriesを利用することで,38行目では検証データの中の正解データy_testと予測データy_test_predictの時系列を作成しています.スライスの中では,真ん中のデータ行dt1=int(len(data)/2)以降を取り出すように範囲が指定されています(29行目).

40行目では2017年1月の1ヶ月分(744時間分)の訓練データの中の正解データy_trainと予測データy_train_predictの時系列図を作成しています.スライスの中のdtm1=744が訓練用データの開始1ヶ月分のデータの行数を示しています.42行目では2018年1月の1ヶ月分(744時間分)の検証データの中の正解データy_testと予測データy_test_predictの時系列を作成しています.40行目と42行目のスライスの中では,検証用データの開始1ヶ月分のデータの範囲が指定されています.

それでは以下のセルを実行してみましょう.

In [8]:
# 予測モデルの結果(testおよびtrain)から時系列図を作成します
def timeseries(datetime_x, y1, y2): 
    # dfのインデックス(時間)をXとする
    X=datetime_x
    # dfのname1列を指定してデータを取り出し,numpy配列で値をY1に与える.
    Y1=y1
    # dfのname1列を指定してデータを取り出し,numpy配列で値をY2に与える.
    Y2=y2
    # 時系列図の大きさを指定
    plt.figure(figsize=(20, 10))
    # y_obsvの時系列図
    plt.plot(X,Y1,color='blue',label='observed')
    # y_frcstの時系列図 
    plt.plot(X,Y2,color='red',label='predicted')
    # グラフのタイトル
    plt.title("Timeseries")
    # x軸のラベル
    plt.xlabel('Time')
    # y軸(左側の第1軸)のラベル
    plt.ylabel('Electric Power [$10^6$kW]')
    # 凡例(左上に置く) 
    plt.legend(loc='upper left')
    return

# datetime_xとする
datetime_x=list(data.index)
# 2018年1月1日のデータの行番号を取得する
# 学習用データと検証用データに別けるために,全体の行(len(data))を2で割る
dt1=int(len(data)/2)
# 学習用データの1ヶ月分の行
dtm1=744
# 2018年2月1日のデータの行番号を取得する
dtm2=dt1+dtm1
# データを用いて時系列図を作成します.
# 訓練データ(正解データと予測データ)の1年分の時系列図 
timeseries(datetime_x[:dt1],y_train,y_train_predict)
# 検証データ(正解データと予測データ)の1年分の時系列図 
timeseries(datetime_x[dt1:],y_test,y_test_predict)
# 訓練データ(正解データと予測データ)の1ヶ月分の時系列図 
timeseries(datetime_x[:dtm1], y_train[:dtm1], y_train_predict[:dtm1])
# 検証データ(正解データと予測データ)の1ヶ月分の時系列図 
timeseries(datetime_x[dt1:dtm2],y_test[:dtm1], y_test_predict[:dtm1])

ニューラルネットワークにより予測された電力消費量の時系列図の結果を見てみると,正解データと予測データは良い一致を示していることが分かります.「気温」,「曜日」,「時刻」というたった3つのパラメータで,季節スケールの変化,週スケールの変化,日スケールの変化を適切に表現できていることが分かります.前回の診断的分析の議論の結果ともよく整合していることが分かります.非線形的な変化を表現できるニューラルネットワークを用いることで,「夏季」「12~15時」「平日」といった条件別けをしてそれぞれに予測モデルを作らなくても一気通貫で予測できることが明らかとなりました.

次に,横軸に正解データの電力消費量と予測データの電力消費量との間を取ることで,各時刻のデータをプロットした散布図matplotlibで作成してみましょう.表計算ソフトExcelの場合には,2つの列データを選択して,ツールバーの中の散布図ボタンを押すプロセスです.

予測モデルの精度が高いときには,散布図上では,一方が上がれば一方が上がる関係(正の相関)となるはずです.予測モデルの精度が低いときには,このように一本の線に乗らず全体的に散らばったような分布となるはずです.

次のセルは,散布図を作成する関数scatter(2行目)を定義し,それを呼び出すことで訓練データの中の正解データy_train予測データy_train_predictedから散布図を作成しています(23行目).関数scatterの中では,まず,図の大きさを指定し(6行目)とx(ここではy_train)とy(ここではy_train_predicted)の散布図を作成しています(8行目).

比較となる,$y=x$の線を引いています(10行目).散布図の点がこの線の上に乗ると高精度な予測モデルであると言えます.12行目では,先に計算した予測モデルの決定係数$R^2$を文字型scoreとして図の左上に置いています.13行目から18行目は,グラフのタイトル,横軸Xのラベル,縦軸Yのラベルを設定しています.

同様の要領で,25行目では,関数scatterを利用して,検証データの中の正解データy_test予測データy_test_predictedから散布図を作成しています.

それでは以下のセルを実行してみましょう.

In [9]:
# 予測モデルの結果(testおよびtrain)から散布図を作成します
def scatter(x, y,score): 
    # 文字列"R2=score"
    score="R2="+str(score)
    # 散布図の大きさを指定
    plt.figure(figsize=(8, 8))
    # 散布図のプロット
    plt.plot(x, y, 'o')
    # Y=Xの線を引く
    plt.plot(x, x)
    # 文字列R2=r2)"を図の左上に置く
    plt.text(np.nanmin(x), np.nanmax(x), score)
    # グラフのタイトル
    plt.title("Scatter diagram")
    # x軸のラベル
    plt.xlabel('Observed Electric Power [$10^6$kW]')
    # y軸のラベル
    plt.ylabel('Predicted Electric Power [$10^6$kW]')    
    return

# データを用いて散布図を作成します.
# 訓練データ(正解データと予測データ)の1年分の時系列
scatter(y_train, y_train_predict,score_train)
# 検証データ(正解データと予測データ)の1年分の時系列
scatter(y_test, y_test_predict,score_test)

電力消費量の正解(横軸)と予測(縦軸)の散布図を見ると,$y=x$の直線の周りに分布し,決定係数$R^2$のスコアも0.9に近く,ニューラルネットワークの高いパフォーマンスが確認できました.しかし,実用面を考えるとより精度が高くなることが望まれます.

是非,皆さんも,説明変数を変更したり,ニューラルネットワークのハイパーパラメータを変更したりすることで,より精度の高い予測モデルの開発を目指して下さい.

2.3 まとめ

● 予測モデルの開発のためには,データ分析により説明変数と目的変数の設計を行います.
● データを適切に要約・整理します.
● データを適切に標準化します.
● 線形重回帰モデルやニューラルネットワークモデルなどの統計的モデルを利用します.
● 統計的モデルにより重み係数やバイアスの最適化を行います.
● 決定係数$R^2$により予測モデルの精度を検証します.
● 時系列図や散布図を確認することで予測モデルの特性を考察し,更なる高精度化に向けて検討します.

高精度な気象予測データと皆さん自身で開発した予測モデルを組み合わせることによって,これまで誰もできなかった新しいビジネスが実現できるようになるかもしれませんね.