次に,過去の気象データと過去の電力消費量データを分析することで統計的モデル(ニューラルネットワーク)を構築して将来の電力消費量を予測します.このような分析を予測的分析(Predictive Analytics)と言います.予測的分析により,例えば,商品ごとの需要を予測分析することによって,需要が上がるタイミングで商品を発注し,そして需要が下がるタイミングで値下げを行い,企業の潜在的な商機やリスクを特定することが可能となります.前回のアメダス気象データ分析チャレンジ!(Python版)診断的分析の理解を前提としていますので,必要があれば前回の資料を復習しながら読み進めてください.
今回利用する ニューラルネットワーク は,教師あり学習 の一種です.教師あり学習では,教師データ(正解データ) による学習(訓練)が必要となります.例えば,スパムメールフィルタでは,事前にそのメールであるかどうかの正解を教えてあげる必要があります.教師あり学習には 「分類(classification)問題」 と 「回帰(regression)問題」 があります.「分類問題」とは,例えば「気温(24℃)」という入力データから「アイスが売れる(1)」か「アイスが売れない(0)」といった2クラスや「今日売れそうな商品(ビール,コーラ,コーヒー,お茶・・・)」といった多クラスを分類するような問題を指します.また,「回帰問題」とは,例えば「気温(24℃)」という入力データから「アイスが何個売れる(100個)」という出力値を予測するような問題を指します.
前回の診断的分析の中でも気温と電力消費量との間の散布図を作成し,線形回帰モデル(ライブラリscikit-learn(sklearn)
)により回帰直線を引きました.これも回帰問題です.線形回帰モデルでは説明変数と目的変数との間に線形の関係があることを仮定されていますが,非線形な関係であってもその関係を表現できることができるのがニューラルネットワーク(非線形回帰モデル) です.
そもそもモデルとは,入力値(説明変数)$x$(例えば,気温)を与えることで,出力値(目的変数)$y=f(x)$(例えば,電力消費量)を得るような関数$f$ を指します.
図:モデルとは?
ニューラルネットワークは,今流行の人工知能(AI)の基本的な技術となっていますので,まずはニューラルネットワークの基本的な仕組みから理解していきましょう.
ニューラルネットワークの歴史を簡単に振り返ると,まず,1943年に 「ニューロン」 という人間の脳を模した数理モデルが提唱され,1957年にデータによる学習が可能な 「パーセプトロン」 というニューラルネットワークが開発されます.パーセプトロンは1960年代に 第1次ブーム を引き起こしますが,XOR問題という弱点が指摘され 冬の時代 を迎えます.その後,1986年に誤差逆伝播法が提案されると,再び 第2次ブーム を迎えますが,コンピュータの能力が十分でなかったり,学習に必要なデータがなかったり,勾配消失などの問題があったりと,再び 冬の時代 となります.その後,2006年になり,トロント大学のヒントン教授によるランプ関数ReLU(Rectified Linear Unit)と呼ばれる活性化関数とドロップアウト正則化などの技術革新に加えて,コンピュータ性能の爆発的な向上により大規模なニューラルネットワークと大量のデータを処理できるようになり現在の 第3次ブーム を迎えることになります.現在も続く第3次ブームで,多層構造をもち「畳み込みネットワーク」を有するニューラルネットワークが広く普及するようになり,ディープラーニング(深層学習) と呼ばれるようになりました.ディープラーニングは音声認識や画像診断の分野で広く利用されるようになり産業応用分野でも急速に利用が広がっています.
図:ニューラルネットワークの歴史
ニューラルネットワークでは,人間の脳内にある 神経細胞(ニューロン) とそのつながり,つまり神経回路網を人工ニューロンという数式的なモデルで表現しています.ニューロンが他のニューロンと繋がっている状態を シナプス結合 といい,電気信号が伝わることによって情報的なつながり ネットワーク を形成しています.シナプスとシナプスの間の情報伝達効率は,前後のニューロンの活動状態によって決まり保持されます.このようなシナプスの情報伝達効率が生物の記憶や学習において重要な役割を果たしています.人間の脳では,1個のニューロンにおよそ1000個のニューロンが結合し,全部で1000億個のニューロンが存在していると言われており,人間の複雑な知性を司っていると考えられています.
図:ニューロンの構造
神経細胞(ニューロン)をコンピュータ上でモデル化したものを人工ニューロンやパーセプトロンともいいます.まず,中央のニューロン(破線の円)に着目すると,まず,入力側($u$)では$m$個のニューロンから $u = b + w_1 \times x_1 + w_2 \times x_2 + \cdots + w_m \times x_m = b + \sum_{i=1}^{m} w_m x_m $ という値を受け取っています.ここで,$x_1、x_2、x_3$ が前のニューロンから伝わってくる電気信号に当たり,$w_1、w_2、w_3$がニューロン間の情報伝達効率に当たります.多数のニューロンから伝達される入力$x$ と 重み$w$の線形結合に バイアス$b$ と呼ばれる定数を足すことで,ニューロンに入る電気信号の強さが決まります.また,中央のニューロンの出力側($o$)に注目すると,ニューロンが蓄積した電気量$y=f(u)$がある一定の閾値aを超えるかどうかによって,1または0どちらかの値が出力されています.1が出力される場合はニューロンの発火を,0が出力される場合にはニューロンの非発火を意味します.ここで,$f$は活性化関数と呼ばれる関数であり,蓄積した電気量をもとに発火(1)するか非発火(0)かを判断する関数です.
図:パーセプトロンの構造
このような単純な計算をしているニューロンが複雑に組み合わさりネットワークを組むことで,次々と電気信号が下流側へと伝達され,高度な「分類」や「回帰」が可能となります.これが ニューラルネットワーク です.ニューラルネットワークを理解する上で重要なのは,「入力層」「中間層(隠れ層)」「出力層」 です.そして,ニューロン同士をつなぐ 「重み」 です.ここで,「重み」 は,情報と情報のつながりの強さを表すパラメータです.この「重み」の大きさは,「入力層」と「出力層」に教師データを入力して 学習 することによって設定されます.情報(説明変数) が入力されるニューロンは 入力層 といいます.複数の説明変数があれば層は複数ニューロンで構成されます.入力層から情報を受け取った 「中間層(隠れ層)」 のニューロンは受け取った情報を自分なりに評価して,伝言ゲームのようにして,その結果を次の層のニューロンに伝達します.それを繰り返した末に,最終的に,情報(目的変数) が 「出力層」 に伝わります.このような,入力層から出力層へと一方向的に情報が伝わるニューラルネットワークのことを,「フィードフォワード(順伝播型)ニューラルネットワーク」 という.
下のフィードフォワードニューラルネットワークの例では,傘の画像の1つ1つのピクセルの輝度が配置された入力層から入っていきます.出力層には,「傘」「ゴルフクラブ」「杖」のニューロンが配置されています.入力層から入った情報は,中間層を伝って,最終的には,出力層の中の「傘」のニューロンに大きな電気信号(情報)が伝わり,画像が80%の確率で「傘」であると認識できたわけです(分類問題の場合).
図:フィードフォワードニューラルネットワークの構造
さて,精度の高いニューラルネットワークを作る上で,ニューロン同士をつなぐ多数の重み$w$の 学習 が重要となってきます(正確にはバイアス$b$も学習の対象となりますが,重み$w$の一種として考えることとしましょう).重み$w$の学習とは,様々な画像に対応する出力値と正解値との間の 誤差が最小 となるに重み$w$をチューニングしていく過程のことを指します.順伝播によって得られた出力値と予め用意された正解値との間で評価された誤差$E$(下の例では2乗和誤差)から,層を逆方向に遡るようにして重み$w$の値を少しずつ 更新 していきます.伝言ゲームの答え合わせをするようなイメージです.たくさんの教師データによりこのような更新の手続きを繰り返すことで次第にニューラルネットワークの誤差が最小になるように 最適化 されていきます.
図:ニューラルネットワークの学習
さて,このような誤差$E$が最小になるように重み$w$をチューニングする方法について考えてみましょう.一般的に,誤差$E$を重み$w$の関数$E(w)$として考え,この誤差関数$E(w)$が最小 となるような$w$の近似値を求めるというアプローチをとります.この近似値を求める際に使われる手法が,「勾配降下法」 と 「誤差逆伝播法」 です.ここで,勾配降下法 とは,重み$w$の初期値における誤差関数の勾配ベクトル(傾斜の向きと大きさ)を求め,その勾配ベクトルの成分が正の場合には重み$w$を小さくするように更新し,勾配ベクトルの成分が負の場合には重み$w$を大きくするように値の更新を行います.次に,1回目の更新後の重み$w$における勾配ベクトルを求めて,同様に勾配の正負に応じてさらに重み$w$の更新を行います.この作業を,重み$w$の更新ができなくなるまで,すなわち,誤差$E$が最小となるまで繰り返していくのです.
下の例では,2つの重み$w_1$と$w_2$の勾配降下法による更新のプロセスを模式化しています.傾斜のある地面に沿ってボールを転がしていくようなイメージで,誤差$E$(傾斜のある地面)がより低くなるような重み$w$の場所(ボールが止まる場所)を探すのです.
勾配降下法による 最適化アルゴリズム には,SGD法,Momentum法,AdaGrad法,Adam法など様々なものがあります.
図:勾配降下法の仕組み.
「勾配降下法」によって重み$w$の最適化を行うためには,上の図のように重み$w$の勾配に関する情報が必要となります.この勾配を効率的に計算する手法が 「誤差逆伝播法」(Backpropagation) と呼ばれるものです.出力側に近い勾配の値から順番に遡って入力側の勾配を求めていくことから,順伝播(フィードフォワード)の対義語として 「逆伝播」 と呼んでいます.
誤差関数$E$は多数の重み$w$で表現される合成関数であるため,多変数関数の連鎖律(チェーンルール)により,構成する各関数の導関数(勾配)の積によって表現することができます.そのため,層の深いネットワークの場合,入力層に近いところの重み$w$の勾配$\frac{\partial E}{\partial w}$を直接求めようとすると,多数の勾配の積を求めなくてはならず計算量が膨大となります.しかしながら,「誤差逆伝播法」では,出力層から入力層に遡るように,重み$w$の勾配$\frac{\partial E}{\partial w}$や入力$u$の勾配$\frac{\partial E}{\partial u}$をメモリに格納しながら効率的に計算することができます.
ただし,層の深いネットワークの場合,誤差逆伝播法によって重み$w$の勾配の値が0になってしまう 勾配消失 という問題があり,ニューラルネットワークの冬の時代をもたらす原因となっていましたが,前述した2006年のヒントン教授による技術革新(活性化関数ReLUとドロップアウト正則化)で解決され,今日のディープラーニングのブームのきっかけとなっています.
図:誤差逆伝播法の仕組み.
ニューラルネットワークにおける活性化関数(Activation function)は,ニューロンを興奮させるための関数であり,ニューロンの興奮状態を発火(1)か非発火(0)を表す信号に変換します.この活性化関数は,非線形変換 である必要があり,これによりニューラルネットワークの複雑な表現力を可能にします(線形変換は何回重ねてもても線形の変化しか表現できない).
使われる活性化関数は時代とともに変化しています.初期には,「ステップ関数」 という活性化関数が用いられていましたが,「誤差逆伝播法」が登場してからは 「シグモイド関数(ロジスティック関数)」 や 「tanh関数」 が使われるようになりました.また,最近のディープニューラルネットワークでは勾配消失の問題を解決させた 「ReLU」 がよく使われるようになりました.
また,出力層で使われる活性化関数は,二クラスの分類問題 の場合は 「シグモイド関数(ロジスティック関数)」,多クラスの分類問題 の場合は 「ソフトマックス関数」,回帰問題の場合は 「恒等関数」が挙げられます.
図:様々な活性化関数.
以下,ニューラルネットワークによく使われるシグモイド関数,tanh関数,ReLU関数,恒等関数のグラフをPythonのプログラミングによって示しています.恒等関数を除くいずれも$z=0$付近で 発火 している様子が見て取れます.
ここでは機械学習ライブラリsklearn
(scikit-learn)の中の統計モデルを用いて予測モデルを構築してみましょう.
ライブラリsklearn
(scikit-learn)は,基本的な機械学習のアルゴリズムやデータセットを備えており,単純で定型的なコードで実装できるため手軽に機械学習を始めることができます.sklearn
(scikit-learn)には機械学習に関する6つのパッケージ:1) 回帰問題,2) 分類問題,3) クラスタリング,4) 次元削減,5) モデル評価,6) データ前処理があります.ここでは,「回帰問題」 のための ニューラルネットワークモデルのクラスneural_network.MLPRegressor
を統計モデルとして導入します.また, sklearn
には,「分類問題」 に対しても,ニューラルネットワークを使うこともできます.その場合は,クラスneural_network.MLPClassifier
を導入します.sklearn
(scikit-learn)のニューラルネットワークの詳細については,scikit-learnの公式ホームページ https://scikit-learn.org/stable/modules/neural_networks_supervised.html を参照してください.
neural_network.MLPRegressor()
クラスでは,多数の引数を設定し,試行錯誤でユーザーがハイパーパラメータの最適値を見つけて設定する必要があります.以下は,ハイパーパラメータのデフォルトの設定値を示します.
class sklearn.neural_network.MLPRegressor(hidden_layer_sizes=(100, ), activation='relu', *, solver='adam', alpha=0.0001, batch_size='auto', learning_rate='constant', learning_rate_init=0.001, power_t=0.5, max_iter=200, shuffle=True, random_state=None, tol=0.0001, verbose=False, warm_start=False, momentum=0.9, nesterovs_momentum=True, early_stopping=False, validation_fraction=0.1, beta_1=0.9, beta_2=0.999, epsilon=1e-08, n_iter_no_change=10, max_fun=15000)
以下,主なハイパーパラメータについて説明します.
●solver
は,学習の最適化アルゴリズムを選択します.ここでの最適化アルゴリズムに基づいて予測誤差が最小になるように重み付けを決定します.ここの選択を誤ると学習速度が遅くなったり,最終的な学習結果が最適な場所(最小値)に行き着かない可能性があります.solver=lbfgs
は,準ニュートン法を省メモリにて実現した手法です.1000以下の小さいデータセットの場合に高速で学習できます.solver=sgd
は,確率的勾配降下法(Stochastic Gradient Descent)であり,訓練データをランダムにシャッフルしてミニバッチとして抽出し,重みを更新していきます.この手法は確率的に局所解にはまりにくくなるという長所やオンライン学習で使いやすいという長所を持っていますが,短所として学習率の設定が難しいことが挙げられます.solver=adam
は,同じく確率的勾配降下法をベースとする手法(ADAptive Momentum estimation)で,「AdaGard法」(過去の学習によって更新されてこなかったパラメータを優先的に更新する)と「Momentum法」(更新量の急変を抑え滑らかに更新する)を組み合わせた方法で高い性能で最適化できる手法です.このsolver=adam
がデフォルト設定となっています.
●activate
は,活性化関数を選択します.activate=identity
は,特に活性化関数を何も設定しません.つまり,$h(x)=x$であり,入力された値をそのまま出力として返します.activate=logistic
は,ロジスティック関数(シグモイド関数)です.式では,$h(x)=\frac{1}{1+e^{-x}}$となります.activate=tanh
は,ハイパボリックタンジェントを活性化関数として使います.tanh
はロジスティック関数と同じ特徴を持ちますが,tanh
では出力範囲が-1から1に変換されています.$h(x)= \tanh x = \frac{1-e^{-2x}}{1+e^{-2x}} $となります.activate=relu
は,ランプ関数ReLU(Rectified Linear Unit)と呼ばれます.式で表すと,$h(x)=\max(0,x)$となります.ReLU関数の導関数はは$x>0$で常に1となるので,ニューラルネットワークの学習で長年問題となっていた「勾配消失」という問題が生じないという利点があります.このactivate=relu
がデフォルト設定となっています.
●hidden_layer_sizes
は,中間層の数とニューロンの数を指定します.例えば,中間層が2層10ニューロンずつ配置する場合は(10,10,)のように指定します.中間層のニューロンの数を入力層の次元よりも小さくすると次元圧縮のように機能し,主成分分析のようなことが可能になります.逆に入力層の次元よりも大きくするとスパースなデータ表現を得ることができます.中間層の数を増やし,ニューロンの数を増やすことで,複雑な非線形回帰問題を解くことができます.ただし,ニューロンの数が多く,訓練データが少なすぎると,この少ない既知データに最適化されすぎてしまい,未知のデータへの対応ができなくなってしまう 「過学習」 という問題が生じてしまいます.過学習とは,学習データによる予測精度は高く,検証用データによる予測データは低い状態を指します.過学習が生じた場合には,より多くの訓練データを用意したり,次の 「L2正則化」 と呼ばれるテクニックを使います.
●alpha
は,L2正則化のペナルティ項の係数 を表しています.正則化とは,過学習を防ぐために,誤差関数に重みを付けた制約(重み$w$を小さくするような成約)を与える項(ペナルティ項)を加えます.alpha
は,一般的に0.01~0.00001の間の値を指定して,この値を大きくするほど,学習時により小さな重みが重要視されるようになります.過学習が生じているときにはこのパラメータを上げて試してみましょう.
●batch_size='auto'
は,最適化ソルバー'sgd'や'adam'のためのミニバッチ学習のサイズを設定します.ソルバーが'lbfgs'の場合,分類器はミニバッチを使用しません.auto "に設定されている場合,batch_size=min(200, n_samples)となります(サンプル数と200の小さい方の数字という意味).
●max_iter
は,学習の反復の最大回数を設定します.常に最大回数学習まで計算するわけではなく,途中で学習が完了したと判断された場合はこれよりも早く終了します.
●tol
は,学習の収束値を設定します.10エポック(学習回数のこと)連続でこの値よりもlossやscoreが向上しなかった場合には,学習が収束したと判断して,学習を終了します.大きくしすぎると学習途中で終了してしまいます.
●n_iter_no_change
は,このエポック数で連続して学習が改善しない場合には計算が収束したものとして計算を終了します.デフォルト値は10です.
●random_state
は,乱数発生器の番号をしてします.重み付けやバイアスの初期値やバッチサンプリング(ソルバーsgd
とadam
のみ)で乱数を発生させる場合に,指定した番号の乱数を発生させます.これをNone(デフォルト)とすることで,毎回異なる結果となりますが,乱数発生器の番号を指定することで再現性のある同じ結果が得られます.
●verbose
は,進行状況のメッセージを標準出力として出力するかしないかを指定します.結果がおかしい場合には,True
としてください.デフォルトは,False
です.
その他のハイパーパラメータの設定方法については,scikit-learnの公式ホームページ https://scikit-learn.org/stable/modules/generated/sklearn.neural_network.MLPRegressor.html を参照してください.
まず,最も簡単なニューラルネットワークを作ってみましょう.以下に示すプログラムのように,sklearn
(scikit-learn)を使えば実質10行のプログラムで機械学習ができてしまいます.
ここでは,まず,「1」という入力を「1,2,3」と与えると(説明変数は3つ),その値を2倍して「2, 4, 6」と出力してくれるモデルを作成してみましょう.
以下のプログラムでは,11行目のx_train=np.array([[1],[2],[3]], dtype=np.float)
では,訓練データの 説明変数$x$ を 「1」「2」「3」と3つ設定 しています.それに対応する 目的変数$y$ は,13行目のy_train=np.array([[2],[4],[6]], dtype=np.float)
で,「2」「4」「6」と3つ設定 しています.
15行目のx_test=np.array([[4],[5],[6]], dtype=np.float)
では,検証データの 説明変数$x$を,未知の「4」「5」「6」を3つ設定 しています.そして,対応する 目的変数$y$つまり正解値 は,17行目のy_test=np.array([[8],[10],[12]], dtype=np.float)
で,説明変数を2倍した「8」「10」「12」と3つ設定しています.
これらのデータを用いて,20行目と22行目でmodel1
インスタンスを作成し, 確率的勾配降下法sgd
による最適化計算を行っています.
24行目のy_test_predict=model1.predict(x_test)
で予測された結果が,どの程度正しいのかを27行目のmse=mean_squared_error(y_test, y_test_predict)
で平均二乗誤差を計算することで判定します.
早速,以下のプログラムを動かしてみましょう.入力値を2倍にするモデルができているでしょうか??
# numpyをnpという別名でインポートします.
import numpy as np
# sklearn(scikit-learn)のニューラルネットワークのライブラリをインポートします
from sklearn import neural_network
# sklearn(scikit-learn)の平均二乗誤差のライブラリをインポートします
from sklearn.metrics import mean_squared_error
%precision 3
%matplotlib inline
# 訓練データの説明変数の設定
x_train=np.array([[1],[2],[3]], dtype=np.float64)
# 訓練データの目的変数の設定
y_train=np.array([[2],[4],[6]], dtype=np.float64)
# 検証データの説明変数の設定
x_test=np.array([[4],[5],[6]], dtype=np.float64)
# 検証データの目的変数(正解)の設定
y_test=np.array([[8],[10],[12]], dtype=np.float64)
# neural_networkクラスからmodel1インスタンスを作成
model1 = neural_network.MLPRegressor(solver="sgd",activation="identity", hidden_layer_sizes=(1),max_iter=10000,tol=1e-6,n_iter_no_change=50,random_state=5,verbose=False)
# 予測モデルの最適化
model1.fit(x_train, y_train.ravel())
# モデルによる予測データの作成(検証用データ)
y_test_predict=model1.predict(x_test)
print('y_test=',y_test.ravel(),'y_test_predict=',y_test_predict)
# 平均二乗誤差の評価
mse=mean_squared_error(y_test, y_test_predict)
print('Mean square error=',mse)
y_test= [ 8. 10. 12.] y_test_predict= [ 7.983 9.972 11.962] Mean square error= 0.0008503844159932146
予測された結果を見ると,説明変数「4」「5」「6」 に対して,いずれの 予測値も「7.983」「9.972」「11.962」 となっています.今回は,正解値「8」「10」「12」とほぼ一致する結果が得られています.平均二乗誤差MSEも小さく適切にモデル化できていると言えるでしょう.みなさんも,入力する訓練データや検証データの数値を変えて試してみてください.
さて,実際の気象データと電力消費データからsklearn(scikit-learn)のニューラルネットワークにより統計モデルを作成してみましょう.
気象データとビジネスデータにより統計的モデルを作ってみましょう.このプログラムでは,気象データとして,気象庁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版)診断的分析」の講義内容をご覧下さい.
気象庁のホームページのリンク: http://www.data.jma.go.jp/gmd/risk/obsdl/
今回のチャレンジ!では事前に用意したcsv形式のファイル(amedas.csv)を使用しますので,データをダウンロードしていただく必要はありませんが,異なる地点や期間,異なる気象要素を各自でダウンロードして分析してみて下さい.
また,このプログラムでは,東京電力の電力消費量データを使用します.気象庁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版)診断的分析」の講義内容をご覧下さい.
東京電力でんき予報のリンク: http://www.tepco.co.jp/forecast/index-j.html
今回のチャレンジ!では事前に用意したcsv形式のファイル(tepco.csv)を使用しますので,データをダウンロードしていただく必要はありませんが,異なる地域の電力消費量データやその他のオープンデータを各自でダウンロードして分析してみて下さい.
このプログラムは,本日配布した資料一式の中のpythonプログラムpredictive_analytics.py
を,Jupyter Notebook用に再収録したものです.どちらも同じ内容のpythonプログラムとなっています.このプログラムは,以下のような構成となっています.
1. ライブラリのインポート
2. 気象庁AMeDASの気象データ(csv形式)を読む
3. 東京電力の電力消費量データ(csv形式)を読む
4. 2つのデータを結合して整理する
5. 学習データと検証データの作成
6. データの標準化
7. 統計的モデルの設定と最適化
8. 時系列図と散布図の作成
早速,1つずつ手を動かしてみましょう.
以下のセルから,プログラムが始まります.上のプルダウンメニューが「Code」となっているセルがpythonプログラムに相当します.セルにカーソルを合わせて上にある「Run」ボタンを押して実行して下さい.以降のプログラムは上のプログラムを前提にして組み立てられてますので,上から順番に実行する必要がありますので注意して下さい.
上の「View」タブの中の「Toggle Line Numbers」をクリックすると,セルに行番号が表示されるので表示させて読み進めて下さい.
まずは,必要なライブラリをインポートします.2行目で数値計算ライブラリnumpy
をインポートします.import numpy as np
とすることで,プログラム中ではnpと省略してインポートすることができます.作図のための,matplotlib
とseaborn
も4行目と6行目でインポートします.また,8行目のdatetime
は日付(日時や時間)の処理ができる便利なライブラリで,気象データ分析には必須となります.また,10行目でデータの前処理や整理に不可欠なpandas
もインポートします.機械学習に必要となる標準化などの前処理のためのライブラリpreprocessing
を16行目でインポートします.最後に,18行目の線形重回帰モデル(linear_model
)と20行目のニューラルネットワーク(neural_network
)を作成するのに必要なsklearn
(scikit-learn)もインポートします.
Jupyter Notebookの中だけで有効なマジックコマンド%precison 3
(小数点3桁まで)と%matplotlib inline
(インラインで画像を表示)も忘れずにつけましょう.
# 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
まずは,前回の診断的分析と同様に,気象庁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(
filename:ファイル名(ここでは,amedas.csv)
encoding:文字エーコーディング(ここでは,Shift_JIS)
skiprows:上から何行分を読み飛ばすか(ここでは,5行)
names:列のラベル(ここでは,A列~D列を['date_time', 'T', 'dummy1', 'dummy2', 'Td', 'dummy3', 'dummy4']とつけます)
parse_dates:datetime型で読み込む列(ここでは,'date_time'の列をdatetime型で読み込む)
index_col:インデックスとする列(ここでは,'datetime')
)
この様にして読み込んだデータは,pandasのデータフレーム型で割り当てられて,関数の中のdf
という変数から,メインプログラムではamedas
という変数に戻されて格納されます(29行目).その直後の31行目では,
amedas=amedas.drop(['dummy1', 'dummy2', 'dummy3', 'dummy4'], axis=1)
のメソッドdrop
により,ラベルが'dummy1'~'dummy4'の列が削除されています.axis=1
とすると列の方向で削除されるという意味です.axis=0
とすれば行の方向で削除されます.
それでは実行してみましょう.
# 気象庁アメダスの気温の時系列データを読み込んで,
# 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
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')からなる表が現れれば成功です.dummy1
~dummy4
も消えていることを確認して下さい.これは,amedas
というデータフレームの内容を省略して表示しています.2017年1月1日01時00分00秒~2019年1月1日0時0分0秒までの1時間毎の気温データ(合計17520行)を読み込めました.
次に,同様に,東京電力でんき予報の電力消費量データ(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(
filename:ファイル名(ここでは,tepco.csv)
encoding:文字エーコーディング(ここでは,Shift_JIS)
skiprows:上から何行分を読み飛ばすか(ここでは,3行)
names:列のラベル(ここでは,A列~C列を['date','time','POWER']とする)
parse_dates:datetime型で読み込む列(ここでは,'date'と'time'の2つの列をまとめてdatetime型で読み込む)
index_col:インデックスとする列(ここでは,'datetime')
)
この様にして読み込んだデータは,pandasのデータフレーム型で割り当てられて,関数の中のdf
という変数から,メインプログラムではtepco
という変数に戻されて格納されます(29行目).
それでは実行してみましょう.
# 東京電力の電力消費量の時系列データを読み込んで,
# 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
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行分)を読み込めました.
次に,気象データの時系列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~23行目では,得られた各時刻のHOURから3時ごとの時間区分(キー'HOUR1','HOUR2',…,'HOUR7')に分類し,該当すれば1,該当しなければ0としています.21時以降の場合は全ての時間区分で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)
その後,24行目では,不要になった'HOUR'の列をデータフレームhour
から.drop()
メソッドで削除し,26行目では,pd.merge()
メソッドによりdf
にhour
を結合しています.
次に,31行目では,データフレームdf
のインデックスから .day
により 日にち を取り出し,キーとして'DAY'を割り当て,day
というデータフレームに格納しています.
day=pd.DataFrame({'DAY': df.index.day}, index=df.index)
さらに,32~33行目では,得られた各時刻のDAYから上旬(キー'DAY1'),中旬(キー'DAY2')に分類し,該当すれば1,該当しなければ0としています.下旬については,'DAY1'と'DAY2'のキーが0となることで表現されます.
day['DAY1']=( ( day['DAY'] >= 1 ) & ( day['DAY'] < 11 ) ).astype(int)
day['DAY2']=( ( day['DAY'] >= 11 ) & ( day['DAY'] < 21 ) ).astype(int)
その後,34行目では,不要になった'DAY'の列をデータフレームday
から.drop()
メソッドで削除し,36行目では,pd.merge()
メソッドによりdf
にday
を結合しています.
また,39行目では,データフレームdf
のインデックスから .weekday
により 曜日 を取り出し,キーとしてWEEK
を割り当て,week
というデータフレームに格納しています. .weekday
は,月曜日は0,火曜日は1,水曜日は2,木曜日は3,金曜日は4,土曜日は5,日曜日は6と返します.
week=pd.DataFrame({'WEEK': df.index.weekday}, index=df.index)
さらに,40行目では,得られた各時刻のWEEKから平日(月0,火1,水2,木3,金4)と休日(土5,日6)に分類して,平日ならば0,休日ならば1となるような修正を行います.
week['WEEK']=( week['WEEK'] >= 5 ).astype(int)
その後,42行目では,pd.merge()
メソッドによりdf
にweek
を結合しています.
最後に,47行目では,データフレームdf
のインデックスから .month
により 月 を取り出し,キーとしてMONTH
を割り当て,month
というデータフレームに格納しています.
month=pd.DataFrame({'MONTH': df.index.month}, index=df.index)
さらに,48~50行目では,得られた各時刻のMONTHから12月・1月・2月(キー'DJF'),3月・4月・5月(キー'MAM'),6月・7月・8月(キー'JJA')に分類して,それぞれ該当すれば1,該当しなければ0とします.9月・10月・11月については,全てのキー'DJF','MAM','JJA'が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)
その後,51行目では,不要になった'MONTH'の列をデータフレームmonth
から.drop()
メソッドで削除し,53行目では,pd.merge()
メソッドによりdf
にmonth
を結合しています.
また,54~58行目では,目的変数の電力消費量('POWER')のデータから1日前の発電量(power24, キー:'POWER24') を説明変数として加えています.電力消費量を予測するモデルを作ろうとしているのに,電力消費量を説明変数として加えるのはおかしいのではないかと感じる人もいるかもしれませんが,1日前の過去の電力消費量は予測時点では既知のデータとなっていると仮定して,過去の電力消費量をも目的変数としているわけです.同じ要領で,59~63行目では1時間前の発電量(power1,キー:'POWER1'),64~68行目では2時間前の発電量(power2,キー:'POWER2'),69~73行目では3時間前の発電量(power3,キー:'POWER3') といった説明変数をデータフレームdf
に追加しています.ちなみに,「1時間前の発電量」 は,目的変数であるキー'POWER'の列データを1つ下にずらすことで作成できます.1つ下にずらすには,pandasのshift()
メソッド を使います(61行目).つまり,
df['POWER'].shift(+1)
とすることで,'POWER'列のデータを1つ下にずらした列データとなります.
最後に,再び,75行目ではdropna()
メソッドにより,データフレーム中の欠測値(NaN)のある行を削除します(1つ下にずらすことで上端で欠測が生じるため).
それでは,以下のセルを実行してみましょう.
# 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, …, HOUR7に分ける.
# 該当すれば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=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に分ける.
# 該当すれば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=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に分ける.
# 該当すれば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=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
POWER | T | T-Td | HOUR1 | HOUR2 | HOUR3 | HOUR4 | HOUR5 | HOUR6 | HOUR7 | DAY1 | DAY2 | WEEK | DJF | MAM | JJA | POWER24 | POWER1 | POWER2 | POWER3 | |
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
datetime | ||||||||||||||||||||
2017-01-02 01:00:00 | 2383.0 | 4.8 | 6.2 | 1 | 0 | 0 | 0 | 0 | 0 | 0 | 1 | 0 | 0 | 1 | 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 | 1 | 0 | 0 | 1 | 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 | 1 | 0 | 0 | 1 | 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 | 1 | 0 | 0 | 1 | 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 | 1 | 0 | 0 | 1 | 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 | 0 | 0 | 0 | 1 | 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 | 0 | 0 | 0 | 1 | 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 | 0 | 0 | 0 | 1 | 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 | 0 | 0 | 0 | 1 | 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 | 0 | 0 | 0 | 1 | 0 | 0 | 3211.0 | 3252.0 | 3376.0 | 3473.0 |
17464 rows × 20 columns
データフレームdata
に,予測モデルの作成に必要となる目的変数('POWER') と多数の 説明変数('T', 'T-Td', 'HOUR1', 'HOUR2', 'HOUR3', 'HOUR4', 'HOUR5', 'HOUR6', 'HOUR7', 'DAY1', 'DAY2', 'WEEK', 'DJF', 'MAM', 'JJA', 'POWER24','POWER1','POWER2','POWER3') を1つに集約することができました.次の節では,これらの説明変数から予測モデルの学習用データや検証用データを選択する方法について説明します.
次に,作成されたデータフレームdata
から予測モデルを作成する際に使用する説明変数と目的変数を設定し,また,前半1年間分のデータを学習用データ(x_train,y_train)として,後半1年間分のデータを検証用データ(x_test,y_test)としてnp.ndarray
のfloat型で割り当てます.まず,22~42行目のリストx_cols
で予測モデルの説明変数を設定しています.ここで1つを設定すれば単回帰モデルになりますが,複数を設定すれば重回帰モデルとなります.リストx_cols
は.extend()
メソッドにより,必要な説明変数が追加されていきます.デフォルトでは,3種類(計9個)の説明変数
x_cols=['T','WEEK', 'HOUR1', 'HOUR2', 'HOUR3', 'HOUR4', 'HOUR5', 'HOUR6', 'HOUR7']
が設定されていることになります.前回の診断的分析による知見に基づいて,気温('T') と 曜日('WEEK') と 時刻( 'HOUR1', 'HOUR2', 'HOUR3', 'HOUR4', 'HOUR5', 'HOUR6', 'HOUR7') が指定されています.次に,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行目では,データフレームdf
(data
のこと)から先に設定したリストx_cols
に該当する列の値(.values
)を抽出し,np.ndarray
のfloat(浮動小数点)型x
として割り当てています.また,9行目では,データフレームdf
(data
のこと)からリストy_cols
に該当する列の値(.values
メソッド)を抽出し,同様にnp.ndarray
のfloat(浮動小数点)型y
として割り当てています.ただし,後の予測モデル作成の際に,変数y
は行データである必要があるため,.ravel()
メソッドにより1次元(ベクトル) に変換しています.
抽出された説明変数x
と目的変数y
に対して,numpy
の スライスの機能 を使って,学習用データ(前半)と検証用データ(後半)の境目の行数split
(dt1
のこと)を指定することにより,11行目と13行目のように,
x_train= x[:split]
y_train= y[:split]
とすることで,インデックスが先頭行~split-1
行までが学習用データとなります.一方で,15行目と17行目のように,
x_test= x[split:]
y_test= y[split:]
とすることで,インデックスがsplit
行~最終行までが検証用データとなります.
それでは以下のセルを実行してみましょう.
# データを学習用(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'])
# 日にち
#x_cols.extend(['DAY1','DAY2'])
# 月
#x_cols.extend(['DJF','MAM','JJA'])
# 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, 9) y_train.shape= (8732,) x_test.shape= (8732, 9) y_test.shape= (8732,)
割り当てられたデータのサイズを.shape
メソッドで確認してみましょう.まず,説明変数x
については,学習用データx_train
と検証用データx_test
に対しては,デフォルト設定で8732行10列のデータが割り当てられているはずです(ただし,x_cols
の設定を変えれば列数が変わります).つまりこの場合,入力層は10ニューロンです.また,目的変数y
については,学習用データy_train
と検証用データy_test
に対しては,8732行のデータが割り当てられていれば問題ありません.つまり,出力層は1ニューロンです.
次に,データの標準化を行います.機械学習により予測モデルを作成する場合,入力するデータのスケール変換が必要となります.異なる桁数のデータを扱う場合にはスケール変換はほぼ必須となります.ニューラルネットワークやサポートベクターマシンではこの種の前処理をしないと学習の進みが遅いという問題もあります.そこで,機械学習ライブラリ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に対してそれぞれ標準化を行います.
それでは以下のセルを実行してみましょう.
# 説明変数に対して,平均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. 0. ] [2.7 1. 0. ... 0. 0. 0. ] [2.7 1. 0. ... 0. 0. 0. ]] 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 -0.378] [-1.592 1.581 -0.378 ... -0.378 -0.378 -0.378] [-1.592 1.581 -0.378 ... -0.378 -0.378 -0.378]] 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. 0. ] [4.2 0. 0. ... 0. 0. 0. ] [4.6 0. 0. ... 0. 0. 0. ]] 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 -0.378] [-1.411 -0.633 -0.378 ... -0.378 -0.378 -0.378] [-1.363 -0.633 -0.378 ... -0.378 -0.378 -0.378]]
学習用データの説明変数x_train
と検証用データの説明変数x_test
に対して,平均0,標準偏差1の標準化ができていることが確認できます.
これで入力データの準備が完了しました.ここでは機械学習ライブラリsklearn
(scikit-learn)の中の統計的モデルを用いて予測モデルを構築してみましょう.ここの説明では,ニューラルネットワークモデルのneural_network.MLPRegressor
クラス を統計的モデルとして導入します.neural_network.MLPRegressor
クラスの詳細については,scikit-learnの公式ホームページ( https://scikit-learn.org/stable/modules/generated/sklearn.neural_network.MLPRegressor.html )を御覧ください.
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
は,乱数発生器の番号をしてします.重み付けやバイアスの初期値やバッチサンプリング(ソルバーsgd
とadam
のみ)で乱数を発生させる場合に,指定した番号の乱数を発生させます.これを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])
それでは以下のセルを実行してみましょう.
# 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 = 5532919.08330685 Iteration 2, loss = 5529359.96105019 Iteration 3, loss = 5524792.70771636 Iteration 4, loss = 5518242.22615091 Iteration 5, loss = 5508785.35564959 Iteration 6, loss = 5495166.10625166 Iteration 7, loss = 5476005.75383042 Iteration 8, loss = 5450467.69820793 Iteration 9, loss = 5417429.30722369 Iteration 10, loss = 5376057.52497189 Iteration 11, loss = 5325618.89885591 Iteration 12, loss = 5265419.93560834 Iteration 13, loss = 5194940.98340013 Iteration 14, loss = 5113765.69265194 Iteration 15, loss = 5021374.33048579 Iteration 16, loss = 4917882.43892466 Iteration 17, loss = 4803575.05095018 Iteration 18, loss = 4678694.78156944 Iteration 19, loss = 4543367.99102858 Iteration 20, loss = 4398320.21683405 Iteration 21, loss = 4244085.34072051 Iteration 22, loss = 4081641.41302536 Iteration 23, loss = 3911551.00916618 Iteration 24, loss = 3735048.17126892 Iteration 25, loss = 3553067.41040913 Iteration 26, loss = 3366843.88383689 Iteration 27, loss = 3177749.29759184 Iteration 28, loss = 2986766.74542442 Iteration 29, loss = 2795388.78615781 Iteration 30, loss = 2605303.94016361 Iteration 31, loss = 2417343.63809424 Iteration 32, loss = 2233104.63931920 Iteration 33, loss = 2053833.08817230 Iteration 34, loss = 1880770.04557951 Iteration 35, loss = 1714904.25995904 Iteration 36, loss = 1557075.15679835 Iteration 37, loss = 1408157.39479365 Iteration 38, loss = 1268717.50519734 Iteration 39, loss = 1139461.74080442 Iteration 40, loss = 1020421.10250980 Iteration 41, loss = 912014.50705961 Iteration 42, loss = 813724.75188318 Iteration 43, loss = 725647.96804256 Iteration 44, loss = 647166.14059462 Iteration 45, loss = 577859.38817698 Iteration 46, loss = 516906.15644301 Iteration 47, loss = 463878.66451350 Iteration 48, loss = 417729.12992215 Iteration 49, loss = 377784.69727125 Iteration 50, loss = 343297.03350102 Iteration 51, loss = 313601.57415822 Iteration 52, loss = 287858.07025354 Iteration 53, loss = 265540.59662090 Iteration 54, loss = 246260.40507967 Iteration 55, loss = 229412.51841743 Iteration 56, loss = 214619.30961103 Iteration 57, loss = 201605.35009878 Iteration 58, loss = 190086.22138031 Iteration 59, loss = 179845.05895510 Iteration 60, loss = 170638.20066047 Iteration 61, loss = 162426.28988224 Iteration 62, loss = 155037.78922936 Iteration 63, loss = 148369.87482445 Iteration 64, loss = 142379.55677498 Iteration 65, loss = 136928.79009197 Iteration 66, loss = 132010.20705783 Iteration 67, loss = 127550.93042342 Iteration 68, loss = 123501.57291468 Iteration 69, loss = 119809.67625004 Iteration 70, loss = 116454.40161834 Iteration 71, loss = 113359.67445982 Iteration 72, loss = 110530.19315226 Iteration 73, loss = 107954.63794883 Iteration 74, loss = 105564.29661863 Iteration 75, loss = 103384.54106343 Iteration 76, loss = 101357.42674233 Iteration 77, loss = 99493.34439815 Iteration 78, loss = 97762.33068593 Iteration 79, loss = 96144.67755915 Iteration 80, loss = 94642.53294076 Iteration 81, loss = 93246.45071421 Iteration 82, loss = 91935.13340959 Iteration 83, loss = 90725.04204245 Iteration 84, loss = 89581.62075478 Iteration 85, loss = 88504.93152805 Iteration 86, loss = 87511.60676952 Iteration 87, loss = 86568.97982654 Iteration 88, loss = 85683.63899679 Iteration 89, loss = 84859.11625305 Iteration 90, loss = 84079.01270757 Iteration 91, loss = 83361.48747693 Iteration 92, loss = 82681.20814304 Iteration 93, loss = 82060.56222616 Iteration 94, loss = 81455.91247744 Iteration 95, loss = 80897.12451723 Iteration 96, loss = 80379.45071255 Iteration 97, loss = 79890.38362519 Iteration 98, loss = 79429.87684017 Iteration 99, loss = 78992.82204209 Iteration 100, loss = 78594.85318910 Iteration 101, loss = 78205.44599328 Iteration 102, loss = 77851.85093742 Iteration 103, loss = 77496.44665039 Iteration 104, loss = 77173.41196826 Iteration 105, loss = 76861.97242551 Iteration 106, loss = 76562.01319953 Iteration 107, loss = 76284.07555360 Iteration 108, loss = 76016.98458365 Iteration 109, loss = 75755.99942901 Iteration 110, loss = 75510.23542053 Iteration 111, loss = 75261.31333515 Iteration 112, loss = 75032.52413789 Iteration 113, loss = 74810.80309215 Iteration 114, loss = 74589.71487858 Iteration 115, loss = 74384.99718221 Iteration 116, loss = 74172.82222838 Iteration 117, loss = 73976.09392103 Iteration 118, loss = 73773.08231003 Iteration 119, loss = 73583.43452818 Iteration 120, loss = 73388.86691570 Iteration 121, loss = 73210.29635380 Iteration 122, loss = 73015.03422731 Iteration 123, loss = 72847.54926246 Iteration 124, loss = 72655.10369008 Iteration 125, loss = 72488.30747360 Iteration 126, loss = 72311.19018171 Iteration 127, loss = 72137.74679091 Iteration 128, loss = 71965.10584761 Iteration 129, loss = 71793.20334939 Iteration 130, loss = 71628.89577282 Iteration 131, loss = 71449.08270819 Iteration 132, loss = 71272.44437013 Iteration 133, loss = 71110.79524938 Iteration 134, loss = 70926.11595369 Iteration 135, loss = 70751.60846860 Iteration 136, loss = 70587.91994187 Iteration 137, loss = 70421.34308709 Iteration 138, loss = 70250.08870132 Iteration 139, loss = 70080.20265793 Iteration 140, loss = 69911.21127524 Iteration 141, loss = 69740.86641951 Iteration 142, loss = 69568.90514468 Iteration 143, loss = 69402.39449148 Iteration 144, loss = 69227.13438194 Iteration 145, loss = 69056.86940159 Iteration 146, loss = 68886.18997217 Iteration 147, loss = 68713.85936791 Iteration 148, loss = 68547.11199463 Iteration 149, loss = 68374.00383762 Iteration 150, loss = 68203.41572920 Iteration 151, loss = 68024.12399113 Iteration 152, loss = 67865.71398135 Iteration 153, loss = 67675.58219033 Iteration 154, loss = 67495.62453504 Iteration 155, loss = 67315.67505864 Iteration 156, loss = 67127.41751634 Iteration 157, loss = 66945.40456252 Iteration 158, loss = 66755.23465526 Iteration 159, loss = 66553.03487106 Iteration 160, loss = 66357.05769551 Iteration 161, loss = 66159.65654360 Iteration 162, loss = 65954.54034628 Iteration 163, loss = 65755.19542056 Iteration 164, loss = 65542.80107526 Iteration 165, loss = 65340.69822600 Iteration 166, loss = 65135.16004209 Iteration 167, loss = 64928.00355738 Iteration 168, loss = 64714.04154198 Iteration 169, loss = 64495.12987595 Iteration 170, loss = 64271.49303037 Iteration 171, loss = 64034.65318664 Iteration 172, loss = 63800.02386395 Iteration 173, loss = 63572.35923664 Iteration 174, loss = 63315.99059606 Iteration 175, loss = 63054.79515447 Iteration 176, loss = 62784.45864266 Iteration 177, loss = 62522.56130912 Iteration 178, loss = 62226.08393455 Iteration 179, loss = 61927.93979346 Iteration 180, loss = 61586.70791745 Iteration 181, loss = 61226.61834977 Iteration 182, loss = 60844.50190123 Iteration 183, loss = 60412.77830728 Iteration 184, loss = 59972.98534799 Iteration 185, loss = 59497.86658241 Iteration 186, loss = 58994.11839739 Iteration 187, loss = 58489.15638225 Iteration 188, loss = 57979.40391631 Iteration 189, loss = 57449.17016103 Iteration 190, loss = 56897.39484729 Iteration 191, loss = 56373.47302163 Iteration 192, loss = 55815.99307399 Iteration 193, loss = 55273.82615542 Iteration 194, loss = 54755.25072134 Iteration 195, loss = 54221.69906743 Iteration 196, loss = 53695.20895733 Iteration 197, loss = 53175.13832495 Iteration 198, loss = 52653.74176266 Iteration 199, loss = 52158.15909283 Iteration 200, loss = 51658.36705240 Iteration 201, loss = 51176.85999425 Iteration 202, loss = 50693.75381584 Iteration 203, loss = 50213.35157013 Iteration 204, loss = 49758.71574035 Iteration 205, loss = 49275.59715204 Iteration 206, loss = 48797.21768140 Iteration 207, loss = 48321.05403286 Iteration 208, loss = 47862.62463556 Iteration 209, loss = 47378.99233548 Iteration 210, loss = 46918.18487205 Iteration 211, loss = 46437.58841922 Iteration 212, loss = 45963.09128156 Iteration 213, loss = 45491.08623235 Iteration 214, loss = 45021.04515553 Iteration 215, loss = 44582.80544388 Iteration 216, loss = 44086.86627466 Iteration 217, loss = 43631.70789150 Iteration 218, loss = 43163.67174400 Iteration 219, loss = 42719.20796146 Iteration 220, loss = 42283.60386767 Iteration 221, loss = 41838.22014509 Iteration 222, loss = 41407.07240488 Iteration 223, loss = 40967.90220864 Iteration 224, loss = 40547.93950002 Iteration 225, loss = 40140.07362731 Iteration 226, loss = 39746.62760813 Iteration 227, loss = 39348.48583939 Iteration 228, loss = 38976.54780758 Iteration 229, loss = 38607.41202939 Iteration 230, loss = 38247.41386116 Iteration 231, loss = 37908.62424926 Iteration 232, loss = 37587.43722724 Iteration 233, loss = 37265.79691931 Iteration 234, loss = 36956.86519445 Iteration 235, loss = 36659.64408886 Iteration 236, loss = 36366.97051283 Iteration 237, loss = 36079.94276561 Iteration 238, loss = 35809.96960600 Iteration 239, loss = 35530.87117107 Iteration 240, loss = 35261.04218012 Iteration 241, loss = 35018.03139501 Iteration 242, loss = 34774.58025797 Iteration 243, loss = 34534.17305485 Iteration 244, loss = 34317.84622402 Iteration 245, loss = 34084.08755551 Iteration 246, loss = 33881.07611608 Iteration 247, loss = 33654.37246069 Iteration 248, loss = 33452.24653156 Iteration 249, loss = 33262.85033777 Iteration 250, loss = 33055.06628205 Iteration 251, loss = 32855.38531934 Iteration 252, loss = 32685.30311549 Iteration 253, loss = 32491.29549606 Iteration 254, loss = 32314.95490694 Iteration 255, loss = 32155.15336738 Iteration 256, loss = 32002.97597335 Iteration 257, loss = 31858.91677775 Iteration 258, loss = 31717.79776396 Iteration 259, loss = 31589.36967139 Iteration 260, loss = 31464.31086841 Iteration 261, loss = 31335.07832346 Iteration 262, loss = 31227.32242643 Iteration 263, loss = 31107.50444871 Iteration 264, loss = 31011.79033443 Iteration 265, loss = 30910.16328398 Iteration 266, loss = 30820.92962729 Iteration 267, loss = 30734.61901342 Iteration 268, loss = 30641.66260029 Iteration 269, loss = 30565.57187089 Iteration 270, loss = 30471.13670878 Iteration 271, loss = 30409.89300169 Iteration 272, loss = 30321.82974929 Iteration 273, loss = 30248.08537207 Iteration 274, loss = 30175.90107021 Iteration 275, loss = 30113.97636386 Iteration 276, loss = 30051.86261326 Iteration 277, loss = 29982.16353986 Iteration 278, loss = 29920.04525515 Iteration 279, loss = 29853.51934298 Iteration 280, loss = 29801.94970064 Iteration 281, loss = 29755.45341511 Iteration 282, loss = 29690.64054827 Iteration 283, loss = 29634.48661996 Iteration 284, loss = 29599.02140045 Iteration 285, loss = 29535.94752619 Iteration 286, loss = 29496.24298959 Iteration 287, loss = 29422.98059970 Iteration 288, loss = 29366.62103347 Iteration 289, loss = 29286.56629441 Iteration 290, loss = 29219.18819457 Iteration 291, loss = 29139.19614106 Iteration 292, loss = 29072.48610444 Iteration 293, loss = 28979.27671196 Iteration 294, loss = 28922.57437719 Iteration 295, loss = 28883.85512832 Iteration 296, loss = 28816.08022224 Iteration 297, loss = 28768.24530470 Iteration 298, loss = 28726.59497838 Iteration 299, loss = 28682.80430879 Iteration 300, loss = 28647.19729557 Iteration 301, loss = 28603.95970532 Iteration 302, loss = 28560.93092236 Iteration 303, loss = 28524.36962064 Iteration 304, loss = 28496.41815214 Iteration 305, loss = 28453.14321607 Iteration 306, loss = 28433.98559588 Iteration 307, loss = 28403.88425993 Iteration 308, loss = 28380.59574798 Iteration 309, loss = 28325.85414556 Iteration 310, loss = 28305.43825068 Iteration 311, loss = 28258.31709778 Iteration 312, loss = 28230.79710651 Iteration 313, loss = 28195.27992645 Iteration 314, loss = 28176.43549292 Iteration 315, loss = 28132.34791860 Iteration 316, loss = 28109.31021625 Iteration 317, loss = 28089.22050505 Iteration 318, loss = 28065.54996589 Iteration 319, loss = 28011.91611263 Iteration 320, loss = 27995.66335932 Iteration 321, loss = 27963.76428528 Iteration 322, loss = 27929.74703919 Iteration 323, loss = 27909.21819138 Iteration 324, loss = 27872.30438857 Iteration 325, loss = 27844.96256316 Iteration 326, loss = 27822.93923202 Iteration 327, loss = 27796.49914549 Iteration 328, loss = 27768.18860960 Iteration 329, loss = 27771.10597529 Iteration 330, loss = 27723.71522393 Iteration 331, loss = 27701.60691083 Iteration 332, loss = 27684.87074307 Iteration 333, loss = 27664.78694075 Iteration 334, loss = 27639.18842203 Iteration 335, loss = 27620.35272130 Iteration 336, loss = 27578.76154651 Iteration 337, loss = 27582.81926754 Iteration 338, loss = 27543.70315645 Iteration 339, loss = 27518.88703642 Iteration 340, loss = 27494.72372665 Iteration 341, loss = 27481.85824716 Iteration 342, loss = 27458.25620499 Iteration 343, loss = 27438.54634333 Iteration 344, loss = 27448.77728428 Iteration 345, loss = 27385.56989415 Iteration 346, loss = 27361.22672931 Iteration 347, loss = 27344.70485964 Iteration 348, loss = 27313.66209602 Iteration 349, loss = 27291.18942876 Iteration 350, loss = 27282.42351257 Iteration 351, loss = 27259.42911991 Iteration 352, loss = 27241.33179277 Iteration 353, loss = 27225.27128686 Iteration 354, loss = 27203.57894630 Iteration 355, loss = 27169.76118257 Iteration 356, loss = 27154.84214175 Iteration 357, loss = 27137.35977692 Iteration 358, loss = 27127.81937691 Iteration 359, loss = 27091.47133020 Iteration 360, loss = 27068.32917412 Iteration 361, loss = 27061.21369138 Iteration 362, loss = 27034.76438147 Iteration 363, loss = 27006.90746040 Iteration 364, loss = 26993.05973689 Iteration 365, loss = 26977.06576795 Iteration 366, loss = 26955.82258263 Iteration 367, loss = 26942.95314348 Iteration 368, loss = 26929.01295124 Iteration 369, loss = 26920.69256523 Iteration 370, loss = 26894.50786370 Iteration 371, loss = 26873.25487920 Iteration 372, loss = 26863.65562887 Iteration 373, loss = 26852.85895305 Iteration 374, loss = 26824.15774957 Iteration 375, loss = 26823.73577282 Iteration 376, loss = 26794.48862393 Iteration 377, loss = 26771.75857844 Iteration 378, loss = 26751.41502844 Iteration 379, loss = 26738.00076690 Iteration 380, loss = 26729.97778458 Iteration 381, loss = 26695.14345160 Iteration 382, loss = 26676.91703866 Iteration 383, loss = 26657.80468923 Iteration 384, loss = 26628.87591687 Iteration 385, loss = 26627.95157960 Iteration 386, loss = 26593.72373282 Iteration 387, loss = 26573.09590591 Iteration 388, loss = 26567.42610868 Iteration 389, loss = 26557.44376995 Iteration 390, loss = 26530.69743487 Iteration 391, loss = 26499.47446217 Iteration 392, loss = 26493.47724349 Iteration 393, loss = 26461.27777197 Iteration 394, loss = 26469.66731068 Iteration 395, loss = 26440.94418324 Iteration 396, loss = 26416.27172838 Iteration 397, loss = 26414.17375414 Iteration 398, loss = 26375.68115115 Iteration 399, loss = 26373.72989775 Iteration 400, loss = 26347.47559160 Iteration 401, loss = 26325.62295744 Iteration 402, loss = 26321.30016351 Iteration 403, loss = 26308.15668385 Iteration 404, loss = 26279.97535592 Iteration 405, loss = 26275.14757659 Iteration 406, loss = 26254.18020489 Iteration 407, loss = 26228.79205931 Iteration 408, loss = 26190.92715254 Iteration 409, loss = 26154.97288546 Iteration 410, loss = 26133.75838172 Iteration 411, loss = 26101.63286111 Iteration 412, loss = 26101.21606662 Iteration 413, loss = 26076.59752128 Iteration 414, loss = 26055.19130269 Iteration 415, loss = 26033.83650084 Iteration 416, loss = 26032.40788128 Iteration 417, loss = 26011.72207681 Iteration 418, loss = 26003.02701578 Iteration 419, loss = 26009.83351247 Iteration 420, loss = 25989.81391990 Iteration 421, loss = 25970.78420588 Iteration 422, loss = 25960.20011779 Iteration 423, loss = 25945.71302981 Iteration 424, loss = 25948.04617523 Iteration 425, loss = 25949.24089369 Iteration 426, loss = 25935.22189807 Iteration 427, loss = 25937.81558238 Iteration 428, loss = 25925.55132853 Iteration 429, loss = 25925.90829218 Iteration 430, loss = 25916.13261024 Iteration 431, loss = 25903.39944175 Iteration 432, loss = 25879.24429122 Iteration 433, loss = 25889.28099986 Iteration 434, loss = 25877.23155796 Iteration 435, loss = 25882.68455777 Iteration 436, loss = 25862.04723612 Iteration 437, loss = 25868.90459661 Iteration 438, loss = 25845.82570131 Iteration 439, loss = 25857.72186330 Iteration 440, loss = 25846.54040497 Iteration 441, loss = 25835.36786719 Iteration 442, loss = 25822.44476199 Iteration 443, loss = 25824.87078707 Iteration 444, loss = 25824.08177790 Iteration 445, loss = 25802.22165400 Iteration 446, loss = 25804.55322724 Iteration 447, loss = 25795.99870852 Iteration 448, loss = 25786.80443072 Iteration 449, loss = 25778.08997011 Iteration 450, loss = 25795.83934567 Iteration 451, loss = 25780.21012879 Iteration 452, loss = 25759.95200255 Iteration 453, loss = 25748.25249322 Iteration 454, loss = 25748.93116353 Iteration 455, loss = 25750.34183202 Iteration 456, loss = 25747.17106457 Iteration 457, loss = 25754.47003538 Iteration 458, loss = 25743.19943102 Iteration 459, loss = 25721.43196087 Iteration 460, loss = 25721.45271120 Iteration 461, loss = 25712.01810271 Iteration 462, loss = 25693.42159229 Iteration 463, loss = 25698.04298405 Iteration 464, loss = 25693.84961379 Iteration 465, loss = 25682.03018258 Iteration 466, loss = 25679.83394642 Iteration 467, loss = 25685.16969146 Iteration 468, loss = 25674.20679049 Iteration 469, loss = 25668.67408771 Iteration 470, loss = 25660.94999459 Iteration 471, loss = 25655.75941282 Iteration 472, loss = 25649.50734625 Iteration 473, loss = 25630.12871355 Iteration 474, loss = 25649.93510555 Iteration 475, loss = 25633.20881012 Iteration 476, loss = 25623.02295952 Iteration 477, loss = 25612.49189560 Iteration 478, loss = 25616.76106225 Iteration 479, loss = 25601.83991999 Iteration 480, loss = 25597.90332214 Iteration 481, loss = 25596.87748192 Iteration 482, loss = 25594.64636317 Iteration 483, loss = 25591.00886941 Iteration 484, loss = 25594.96398318 Iteration 485, loss = 25574.35805817 Iteration 486, loss = 25570.42723982 Iteration 487, loss = 25564.94056423 Iteration 488, loss = 25567.71797542 Iteration 489, loss = 25556.84305181 Iteration 490, loss = 25552.77234747 Iteration 491, loss = 25553.55688146 Iteration 492, loss = 25530.85308514 Iteration 493, loss = 25533.57686848 Iteration 494, loss = 25534.88325880 Iteration 495, loss = 25528.37573945 Iteration 496, loss = 25525.43413133 Iteration 497, loss = 25510.98490038 Iteration 498, loss = 25519.01721996 Iteration 499, loss = 25516.45141698 Iteration 500, loss = 25499.69907327 Iteration 501, loss = 25512.48559400 Iteration 502, loss = 25489.35385476 Iteration 503, loss = 25486.40709811 Iteration 504, loss = 25477.80291586 Iteration 505, loss = 25497.71678994 Iteration 506, loss = 25478.76926970 Iteration 507, loss = 25470.50151931 Iteration 508, loss = 25467.67964768 Iteration 509, loss = 25462.08127982 Iteration 510, loss = 25465.09565516 Iteration 511, loss = 25454.66521103 Iteration 512, loss = 25452.47999306 Iteration 513, loss = 25446.63091412 Iteration 514, loss = 25470.93078236 Iteration 515, loss = 25439.94520522 Iteration 516, loss = 25427.69041543 Iteration 517, loss = 25435.95845839 Iteration 518, loss = 25418.84995961 Iteration 519, loss = 25425.59410262 Iteration 520, loss = 25423.46788387 Iteration 521, loss = 25437.90245157 Iteration 522, loss = 25420.33543905 Iteration 523, loss = 25412.26852669 Iteration 524, loss = 25417.72832428 Iteration 525, loss = 25417.13351334 Iteration 526, loss = 25399.22633979 Iteration 527, loss = 25387.78294624 Iteration 528, loss = 25394.21022136 Iteration 529, loss = 25384.41798181 Iteration 530, loss = 25375.19790394 Iteration 531, loss = 25390.33043670 Iteration 532, loss = 25377.29727747 Iteration 533, loss = 25381.64428521 Iteration 534, loss = 25376.02586368 Iteration 535, loss = 25365.78698248 Iteration 536, loss = 25367.20207037 Iteration 537, loss = 25354.59830082 Iteration 538, loss = 25361.89393413 Iteration 539, loss = 25357.72636975 Iteration 540, loss = 25343.96043983 Iteration 541, loss = 25343.70453873 Iteration 542, loss = 25333.24940135 Iteration 543, loss = 25333.69272809 Iteration 544, loss = 25331.60969745 Iteration 545, loss = 25333.95408750 Iteration 546, loss = 25335.72874293 Iteration 547, loss = 25320.91002667 Iteration 548, loss = 25356.05038290 Iteration 549, loss = 25310.83576156 Iteration 550, loss = 25338.47930501 Iteration 551, loss = 25324.82324756 Iteration 552, loss = 25300.51228252 Iteration 553, loss = 25300.12542236 Iteration 554, loss = 25300.78973632 Iteration 555, loss = 25295.69511300 Iteration 556, loss = 25284.93007328 Iteration 557, loss = 25310.35650979 Iteration 558, loss = 25276.64072082 Iteration 559, loss = 25286.80131237 Iteration 560, loss = 25295.53533824 Iteration 561, loss = 25267.64826531 Iteration 562, loss = 25271.78036647 Iteration 563, loss = 25283.15516671 Iteration 564, loss = 25271.68005664 Iteration 565, loss = 25270.85736219 Iteration 566, loss = 25249.40849138 Iteration 567, loss = 25253.63581260 Iteration 568, loss = 25248.32201360 Iteration 569, loss = 25253.21292655 Iteration 570, loss = 25247.06795339 Iteration 571, loss = 25252.03564667 Iteration 572, loss = 25232.27038339 Iteration 573, loss = 25232.44176150 Iteration 574, loss = 25226.67342395 Iteration 575, loss = 25229.59768956 Iteration 576, loss = 25224.05333763 Iteration 577, loss = 25227.75036526 Iteration 578, loss = 25231.78440538 Iteration 579, loss = 25211.28520887 Iteration 580, loss = 25205.10580325 Iteration 581, loss = 25215.40401959 Iteration 582, loss = 25211.38839268 Iteration 583, loss = 25188.29736186 Iteration 584, loss = 25195.03899373 Iteration 585, loss = 25191.15729741 Iteration 586, loss = 25186.27383746 Iteration 587, loss = 25179.10905185 Iteration 588, loss = 25170.85120924 Iteration 589, loss = 25164.69337636 Iteration 590, loss = 25180.59540433 Iteration 591, loss = 25194.04430543 Iteration 592, loss = 25171.21050653 Iteration 593, loss = 25158.33956274 Iteration 594, loss = 25151.50745024 Iteration 595, loss = 25150.63516103 Iteration 596, loss = 25143.81203908 Iteration 597, loss = 25130.64210277 Iteration 598, loss = 25122.41995324 Iteration 599, loss = 25121.00254209 Iteration 600, loss = 25106.56699721 Iteration 601, loss = 25101.09902339 Iteration 602, loss = 25074.65532051 Iteration 603, loss = 25074.72131852 Iteration 604, loss = 25058.03154942 Iteration 605, loss = 25039.19543870 Iteration 606, loss = 25035.94231743 Iteration 607, loss = 25027.47246606 Iteration 608, loss = 25023.11608132 Iteration 609, loss = 25007.81005790 Iteration 610, loss = 25006.88826638 Iteration 611, loss = 25007.85039475 Iteration 612, loss = 24982.21640540 Iteration 613, loss = 25026.09737132 Iteration 614, loss = 24982.99261859 Iteration 615, loss = 24976.93260611 Iteration 616, loss = 24967.60286184 Iteration 617, loss = 24963.40203498 Iteration 618, loss = 24959.34976490 Iteration 619, loss = 24953.76730372 Iteration 620, loss = 24946.15402195 Iteration 621, loss = 24939.73876740 Iteration 622, loss = 24934.62561823 Iteration 623, loss = 24945.89068511 Iteration 624, loss = 24942.78655984 Iteration 625, loss = 24935.66310369 Iteration 626, loss = 24947.16464298 Iteration 627, loss = 24939.13025923 Iteration 628, loss = 24938.25437489 Iteration 629, loss = 24921.93860029 Iteration 630, loss = 24916.82940477 Iteration 631, loss = 24946.57363940 Iteration 632, loss = 24931.29559087 Iteration 633, loss = 24910.44166699 Iteration 634, loss = 24900.90120149 Iteration 635, loss = 24914.84226350 Iteration 636, loss = 24893.87018312 Iteration 637, loss = 24902.83817069 Iteration 638, loss = 24892.56710134 Iteration 639, loss = 24884.59037451 Iteration 640, loss = 24898.64306634 Iteration 641, loss = 24890.12604855 Iteration 642, loss = 24878.60725128 Iteration 643, loss = 24879.43405136 Iteration 644, loss = 24874.57308693 Iteration 645, loss = 24876.54971319 Iteration 646, loss = 24864.52964926 Iteration 647, loss = 24864.36472442 Iteration 648, loss = 24860.77064935 Iteration 649, loss = 24857.22326575 Iteration 650, loss = 24849.74945469 Iteration 651, loss = 24856.53207852 Iteration 652, loss = 24849.52398525 Iteration 653, loss = 24856.52275963 Iteration 654, loss = 24861.85469654 Iteration 655, loss = 24832.95863970 Iteration 656, loss = 24872.31274127 Iteration 657, loss = 24838.35241404 Iteration 658, loss = 24843.38026924 Iteration 659, loss = 24827.60618561 Iteration 660, loss = 24821.04918884 Iteration 661, loss = 24832.68054509 Iteration 662, loss = 24820.78800333 Iteration 663, loss = 24809.76256712 Iteration 664, loss = 24811.41424367 Iteration 665, loss = 24843.72223924 Iteration 666, loss = 24805.59866598 Iteration 667, loss = 24801.87534955 Iteration 668, loss = 24802.36881239 Iteration 669, loss = 24799.48783064 Iteration 670, loss = 24785.01934688 Iteration 671, loss = 24794.64810837 Iteration 672, loss = 24783.22004826 Iteration 673, loss = 24780.01118040 Iteration 674, loss = 24779.17887389 Iteration 675, loss = 24784.43734321 Iteration 676, loss = 24763.56170699 Iteration 677, loss = 24767.07034965 Iteration 678, loss = 24770.05939601 Iteration 679, loss = 24759.60374069 Iteration 680, loss = 24752.36345441 Iteration 681, loss = 24766.59627028 Iteration 682, loss = 24753.86080962 Iteration 683, loss = 24743.99611079 Iteration 684, loss = 24756.49967224 Iteration 685, loss = 24744.78803430 Iteration 686, loss = 24747.46184188 Iteration 687, loss = 24734.53935576 Iteration 688, loss = 24752.47522208 Iteration 689, loss = 24725.36912904 Iteration 690, loss = 24730.63230895 Iteration 691, loss = 24718.85235949 Iteration 692, loss = 24713.09033601 Iteration 693, loss = 24713.78764826 Iteration 694, loss = 24710.04079808 Iteration 695, loss = 24700.80500848 Iteration 696, loss = 24700.39055791 Iteration 697, loss = 24693.53425229 Iteration 698, loss = 24700.25663297 Iteration 699, loss = 24689.79890116 Iteration 700, loss = 24693.84980253 Iteration 701, loss = 24685.58483639 Iteration 702, loss = 24699.72121420 Iteration 703, loss = 24687.72010551 Iteration 704, loss = 24685.94756821 Iteration 705, loss = 24673.18341824 Iteration 706, loss = 24683.86228783 Iteration 707, loss = 24661.47103577 Iteration 708, loss = 24671.32680167 Iteration 709, loss = 24662.00752716 Iteration 710, loss = 24654.39817893 Iteration 711, loss = 24655.41092000 Iteration 712, loss = 24657.30504395 Iteration 713, loss = 24645.61828494 Iteration 714, loss = 24644.41109442 Iteration 715, loss = 24645.38200975 Iteration 716, loss = 24635.80409678 Iteration 717, loss = 24633.10974023 Iteration 718, loss = 24629.91165381 Iteration 719, loss = 24627.08077932 Iteration 720, loss = 24619.94136265 Iteration 721, loss = 24618.92495669 Iteration 722, loss = 24613.55219659 Iteration 723, loss = 24607.08735792 Iteration 724, loss = 24608.01429120 Iteration 725, loss = 24600.65505611 Iteration 726, loss = 24595.76805094 Iteration 727, loss = 24591.68829623 Iteration 728, loss = 24595.69117411 Iteration 729, loss = 24591.64841373 Iteration 730, loss = 24599.24061184 Iteration 731, loss = 24606.08810677 Iteration 732, loss = 24577.06244550 Iteration 733, loss = 24580.77290431 Iteration 734, loss = 24570.61422863 Iteration 735, loss = 24574.59390296 Iteration 736, loss = 24565.11561322 Iteration 737, loss = 24567.60151009 Iteration 738, loss = 24562.57728990 Iteration 739, loss = 24564.79167117 Iteration 740, loss = 24554.43637210 Iteration 741, loss = 24551.19165610 Iteration 742, loss = 24557.28882273 Iteration 743, loss = 24543.55795921 Iteration 744, loss = 24542.79950178 Iteration 745, loss = 24547.13221380 Iteration 746, loss = 24544.47278965 Iteration 747, loss = 24536.69916030 Iteration 748, loss = 24532.46453016 Iteration 749, loss = 24542.74007634 Iteration 750, loss = 24522.38873249 Iteration 751, loss = 24516.09917204 Iteration 752, loss = 24523.06061738 Iteration 753, loss = 24507.32512904 Iteration 754, loss = 24504.33550737 Iteration 755, loss = 24528.37315133 Iteration 756, loss = 24503.87403256 Iteration 757, loss = 24497.56923384 Iteration 758, loss = 24498.03041028 Iteration 759, loss = 24502.82457198 Iteration 760, loss = 24493.86802315 Iteration 761, loss = 24491.04608615 Iteration 762, loss = 24522.09588040 Iteration 763, loss = 24500.43010830 Iteration 764, loss = 24488.24080269 Iteration 765, loss = 24482.01530103 Iteration 766, loss = 24484.02603401 Iteration 767, loss = 24486.00519190 Iteration 768, loss = 24482.90339598 Iteration 769, loss = 24475.11918917 Iteration 770, loss = 24473.28825828 Iteration 771, loss = 24466.63533731 Iteration 772, loss = 24456.57073881 Iteration 773, loss = 24461.20467373 Iteration 774, loss = 24467.90616649 Iteration 775, loss = 24458.14548508 Iteration 776, loss = 24473.82048549 Iteration 777, loss = 24468.75329883 Iteration 778, loss = 24453.59278857 Iteration 779, loss = 24453.50201867 Iteration 780, loss = 24443.37060552 Iteration 781, loss = 24448.93703801 Iteration 782, loss = 24432.92188085 Iteration 783, loss = 24446.36469153 Iteration 784, loss = 24463.19158233 Iteration 785, loss = 24445.01984375 Iteration 786, loss = 24431.01604595 Iteration 787, loss = 24425.21789531 Iteration 788, loss = 24427.64337467 Iteration 789, loss = 24417.21006445 Iteration 790, loss = 24416.87656311 Iteration 791, loss = 24411.01579967 Iteration 792, loss = 24414.12594192 Iteration 793, loss = 24443.50993289 Iteration 794, loss = 24402.25232548 Iteration 795, loss = 24406.74760505 Iteration 796, loss = 24393.53387282 Iteration 797, loss = 24387.08421531 Iteration 798, loss = 24390.81885053 Iteration 799, loss = 24387.47174338 Iteration 800, loss = 24378.98005249 Iteration 801, loss = 24396.90733691 Iteration 802, loss = 24390.01111380 Iteration 803, loss = 24378.11644272 Iteration 804, loss = 24373.92602575 Iteration 805, loss = 24377.39801742 Iteration 806, loss = 24373.05400394 Iteration 807, loss = 24361.87015008 Iteration 808, loss = 24362.79250422 Iteration 809, loss = 24357.96082104 Iteration 810, loss = 24356.44574675 Iteration 811, loss = 24355.95940729 Iteration 812, loss = 24343.90793424 Iteration 813, loss = 24352.66178479 Iteration 814, loss = 24345.48936424 Iteration 815, loss = 24349.93229960 Iteration 816, loss = 24333.59004443 Iteration 817, loss = 24348.68715409 Iteration 818, loss = 24335.32142407 Iteration 819, loss = 24346.11046871 Iteration 820, loss = 24326.55427263 Iteration 821, loss = 24337.60168659 Iteration 822, loss = 24319.59419306 Iteration 823, loss = 24327.79362353 Iteration 824, loss = 24316.43860333 Iteration 825, loss = 24322.19006752 Iteration 826, loss = 24319.01263414 Iteration 827, loss = 24309.94089924 Iteration 828, loss = 24314.24017932 Iteration 829, loss = 24312.75592200 Iteration 830, loss = 24326.26390741 Iteration 831, loss = 24308.56817973 Iteration 832, loss = 24294.91910074 Iteration 833, loss = 24297.12222467 Iteration 834, loss = 24305.33156496 Iteration 835, loss = 24296.95935906 Iteration 836, loss = 24283.55679538 Iteration 837, loss = 24285.89731039 Iteration 838, loss = 24296.86882051 Iteration 839, loss = 24287.71120540 Iteration 840, loss = 24283.63627423 Iteration 841, loss = 24280.81729564 Iteration 842, loss = 24276.22795184 Iteration 843, loss = 24285.89416241 Iteration 844, loss = 24291.16074638 Iteration 845, loss = 24270.47977327 Iteration 846, loss = 24278.11967021 Iteration 847, loss = 24270.97976999 Iteration 848, loss = 24269.95239441 Iteration 849, loss = 24268.46623175 Iteration 850, loss = 24260.01155858 Iteration 851, loss = 24271.35473411 Iteration 852, loss = 24266.55202515 Iteration 853, loss = 24275.60133218 Iteration 854, loss = 24256.77434556 Iteration 855, loss = 24255.54742516 Iteration 856, loss = 24258.84243046 Iteration 857, loss = 24255.21465686 Iteration 858, loss = 24261.58296071 Iteration 859, loss = 24260.02597687 Iteration 860, loss = 24256.25502170 Iteration 861, loss = 24252.28565636 Iteration 862, loss = 24246.90990294 Iteration 863, loss = 24259.07099746 Iteration 864, loss = 24252.48650511 Iteration 865, loss = 24249.94861326 Iteration 866, loss = 24228.69117551 Iteration 867, loss = 24250.80594988 Iteration 868, loss = 24250.75870206 Iteration 869, loss = 24256.27649540 Iteration 870, loss = 24246.46973747 Iteration 871, loss = 24230.67458755 Iteration 872, loss = 24233.08209150 Iteration 873, loss = 24231.16535921 Iteration 874, loss = 24234.77151509 Iteration 875, loss = 24234.47458957 Iteration 876, loss = 24229.70800448 Iteration 877, loss = 24235.05766618 Training loss did not improve more than tol=0.000010 for 10 consecutive epochs. Stopping. train R2 score= 0.8749091825129898 test R2 score= 0.8632855258750649 y_train_predict= [2974.637 3009.637 2907.553 ... 3606.43 3663.414 3663.414] y_test_predict= [3149.635 3205.635 3240.635 ... 3736.532 3800.895 3757.986] 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. 0. ] [2.7 1. 0. ... 0. 0. 0. ] [2.7 1. 0. ... 0. 0. 0. ]]] 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. 0. ] [4.2 0. 0. ... 0. 0. 0. ] [4.6 0. 0. ... 0. 0. 0. ]]]
さて,結果を見てみましょう.ニューラルネットワークモデルの予測精度(決定係数$R^2$)は,学習用データは0.88を,検証用データは0.87となりました.気温,曜日,時間を説明変数とすることで,精度良く電力消費量を予測できると言えるでしょう.
さて,2.4.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.4.5節のセルまで戻る必要がありますので気をつけて下さい.
次に,時系列図を作成して調べてみましょう.予測データと正解データを適切な形にグラフ化することで,その特徴がみえてきます.数値だけでは分かりにくい例外的な特徴も,グラフにすることによって抽出することできる.ここでは,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ヶ月分のデータの範囲が指定されています.
それでは以下のセルを実行してみましょう.
# 予測モデルの結果(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
から散布図を作成しています.
それでは以下のセルを実行してみましょう.
# 予測モデルの結果(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に近く,ニューラルネットワークの高いパフォーマンスが確認できました.しかし,実用面を考えるとより精度が高くなることが望まれます.
是非,皆さんも,説明変数を変更したり,ニューラルネットワークのハイパーパラメータを変更したりすることで,より精度の高い予測モデルの開発を目指して下さい.
● 予測モデルの開発のためには,データ分析により説明変数と目的変数の設計を行います.
● データを適切に要約・整理します.
● データを適切に標準化します.
● 線形重回帰モデルやニューラルネットワークモデルなどの統計的モデルを利用します.
● 統計的モデルにより重み係数やバイアスの最適化を行います.
● 決定係数$R^2$により予測モデルの精度を検証します.
● 時系列図や散布図を確認することで予測モデルの特性を考察し,更なる高精度化に向けて検討します.
高精度な気象予測データと皆さん自身で開発した予測モデルを組み合わせることによって,これまで誰もできなかった新しいビジネスが実現できるようになるかもしれませんね.