2. アメダス気象データ分析チャレンジ!(Python版) 予測的分析

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

2.1 ニューラルネットワークとは?

今回利用する ニューラルネットワーク は,教師あり学習 の一種です.教師あり学習では,教師データ(正解データ) による学習(訓練)が必要となります.例えば,スパムメールフィルタでは,事前にそのメールであるかどうかの正解を教えてあげる必要があります.教師あり学習には 「分類(classification)問題」「回帰(regression)問題」 があります.「分類問題」とは,例えば「気温(24℃)」という入力データから「アイスが売れる(1)」か「アイスが売れない(0)」といった2クラスや「今日売れそうな商品(ビール,コーラ,コーヒー,お茶・・・)」といった多クラスを分類するような問題を指します.また,「回帰問題」とは,例えば「気温(24℃)」という入力データから「アイスが何個売れる(100個)」という出力値を予測するような問題を指します.

前回の診断的分析の中でも気温と電力消費量との間の散布図を作成し,線形回帰モデル(ライブラリscikit-learn(sklearn))により回帰直線を引きました.これも回帰問題です.線形回帰モデルでは説明変数と目的変数との間に線形の関係があることを仮定されていますが,非線形な関係であってもその関係を表現できることができるのがニューラルネットワーク(非線形回帰モデル) です.

そもそもモデルとは,入力値(説明変数)$x$(例えば,気温)を与えることで,出力値(目的変数)$y=f(x)$(例えば,電力消費量)を得るような関数$f$ を指します.

FIG12.png 図:モデルとは?

ニューラルネットワークは,今流行の人工知能(AI)の基本的な技術となっていますので,まずはニューラルネットワークの基本的な仕組みから理解していきましょう.

ニューラルネットワークの歴史を簡単に振り返ると,まず,1943年に 「ニューロン」 という人間の脳を模した数理モデルが提唱され,1957年にデータによる学習が可能な 「パーセプトロン」 というニューラルネットワークが開発されます.パーセプトロンは1960年代に 第1次ブーム を引き起こしますが,XOR問題という弱点が指摘され 冬の時代 を迎えます.その後,1986年に誤差逆伝播法が提案されると,再び 第2次ブーム を迎えますが,コンピュータの能力が十分でなかったり,学習に必要なデータがなかったり,勾配消失などの問題があったりと,再び 冬の時代 となります.その後,2006年になり,トロント大学のヒントン教授によるランプ関数ReLU(Rectified Linear Unit)と呼ばれる活性化関数とドロップアウト正則化などの技術革新に加えて,コンピュータ性能の爆発的な向上により大規模なニューラルネットワークと大量のデータを処理できるようになり現在の 第3次ブーム を迎えることになります.現在も続く第3次ブームで,多層構造をもち「畳み込みネットワーク」を有するニューラルネットワークが広く普及するようになり,ディープラーニング(深層学習) と呼ばれるようになりました.ディープラーニングは音声認識や画像診断の分野で広く利用されるようになり産業応用分野でも急速に利用が広がっています.  

FIG01.png 図:ニューラルネットワークの歴史

ニューラルネットワークでは,人間の脳内にある 神経細胞(ニューロン) とそのつながり,つまり神経回路網を人工ニューロンという数式的なモデルで表現しています.ニューロンが他のニューロンと繋がっている状態を シナプス結合 といい,電気信号が伝わることによって情報的なつながり ネットワーク を形成しています.シナプスとシナプスの間の情報伝達効率は,前後のニューロンの活動状態によって決まり保持されます.このようなシナプスの情報伝達効率が生物の記憶や学習において重要な役割を果たしています.人間の脳では,1個のニューロンにおよそ1000個のニューロンが結合し,全部で1000億個のニューロンが存在していると言われており,人間の複雑な知性を司っていると考えられています.

FIG02.png 図:ニューロンの構造

神経細胞(ニューロン)をコンピュータ上でモデル化したものを人工ニューロンやパーセプトロンともいいます.まず,中央のニューロン(破線の円)に着目すると,まず,入力側($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)かを判断する関数です.

FIG03.png 図:パーセプトロンの構造

このような単純な計算をしているニューロンが複雑に組み合わさりネットワークを組むことで,次々と電気信号が下流側へと伝達され,高度な「分類」や「回帰」が可能となります.これが ニューラルネットワーク です.ニューラルネットワークを理解する上で重要なのは,「入力層」「中間層(隠れ層)」「出力層」 です.そして,ニューロン同士をつなぐ 「重み」 です.ここで,「重み」 は,情報と情報のつながりの強さを表すパラメータです.この「重み」の大きさは,「入力層」と「出力層」に教師データを入力して 学習 することによって設定されます.情報(説明変数) が入力されるニューロンは 入力層 といいます.複数の説明変数があれば層は複数ニューロンで構成されます.入力層から情報を受け取った 「中間層(隠れ層)」 のニューロンは受け取った情報を自分なりに評価して,伝言ゲームのようにして,その結果を次の層のニューロンに伝達します.それを繰り返した末に,最終的に,情報(目的変数)「出力層」 に伝わります.このような,入力層から出力層へと一方向的に情報が伝わるニューラルネットワークのことを,「フィードフォワード(順伝播型)ニューラルネットワーク」 という.

下のフィードフォワードニューラルネットワークの例では,傘の画像の1つ1つのピクセルの輝度が配置された入力層から入っていきます.出力層には,「傘」「ゴルフクラブ」「杖」のニューロンが配置されています.入力層から入った情報は,中間層を伝って,最終的には,出力層の中の「傘」のニューロンに大きな電気信号(情報)が伝わり,画像が80%の確率で「傘」であると認識できたわけです(分類問題の場合).

FIG04.png 図:フィードフォワードニューラルネットワークの構造

さて,精度の高いニューラルネットワークを作る上で,ニューロン同士をつなぐ多数の重み$w$の 学習 が重要となってきます(正確にはバイアス$b$も学習の対象となりますが,重み$w$の一種として考えることとしましょう).重み$w$の学習とは,様々な画像に対応する出力値と正解値との間の 誤差が最小 となるに重み$w$をチューニングしていく過程のことを指します.順伝播によって得られた出力値と予め用意された正解値との間で評価された誤差$E$(下の例では2乗和誤差)から,層を逆方向に遡るようにして重み$w$の値を少しずつ 更新 していきます.伝言ゲームの答え合わせをするようなイメージです.たくさんの教師データによりこのような更新の手続きを繰り返すことで次第にニューラルネットワークの誤差が最小になるように 最適化 されていきます.

FIG05.png 図:ニューラルネットワークの学習

さて,このような誤差$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法など様々なものがあります.

FIG06.png 図:勾配降下法の仕組み.

「勾配降下法」によって重み$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とドロップアウト正則化)で解決され,今日のディープラーニングのブームのきっかけとなっています.

FIG07.png 図:誤差逆伝播法の仕組み.

ニューラルネットワークにおける活性化関数(Activation function)は,ニューロンを興奮させるための関数であり,ニューロンの興奮状態を発火(1)か非発火(0)を表す信号に変換します.この活性化関数は,非線形変換 である必要があり,これによりニューラルネットワークの複雑な表現力を可能にします(線形変換は何回重ねてもても線形の変化しか表現できない).

使われる活性化関数は時代とともに変化しています.初期には,「ステップ関数」 という活性化関数が用いられていましたが,「誤差逆伝播法」が登場してからは 「シグモイド関数(ロジスティック関数)」「tanh関数」 が使われるようになりました.また,最近のディープニューラルネットワークでは勾配消失の問題を解決させた 「ReLU」 がよく使われるようになりました.

また,出力層で使われる活性化関数は,二クラスの分類問題 の場合は 「シグモイド関数(ロジスティック関数)」多クラスの分類問題 の場合は 「ソフトマックス関数」回帰問題の場合は 「恒等関数」が挙げられます.

FIG08.png 図:様々な活性化関数.

以下,ニューラルネットワークによく使われるシグモイド関数,tanh関数,ReLU関数,恒等関数のグラフをPythonのプログラミングによって示しています.恒等関数を除くいずれも$z=0$付近で 発火 している様子が見て取れます.

2.2 sklearn(scikit-learn)のニューラルネットワーク(neural_networkクラス)

ここでは機械学習ライブラリ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は,乱数発生器の番号をしてします.重み付けやバイアスの初期値やバッチサンプリング(ソルバーsgdadamのみ)で乱数を発生させる場合に,指定した番号の乱数を発生させます.これを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倍にするモデルができているでしょうか??

予測された結果を見ると,説明変数「4」「5」「6」 に対して,いずれの 予測値も「7.983」「9.972」「11.962」 となっています.今回は,正解値「8」「10」「12」とほぼ一致する結果が得られています.平均二乗誤差MSEも小さく適切にモデル化できていると言えるでしょう.みなさんも,入力する訓練データや検証データの数値を変えて試してみてください.

さて,実際の気象データと電力消費データからsklearn(scikit-learn)のニューラルネットワークにより統計モデルを作成してみましょう.

2.3 使用するデータ

2.3.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版)診断的分析」の講義内容をご覧下さい.

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

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

2.3.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版)診断的分析」の講義内容をご覧下さい.

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

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

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

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

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

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

2.4.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(インラインで画像を表示)も忘れずにつけましょう.

2.4.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(

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とすれば行の方向で削除されます.

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

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

2.4.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(

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行目).

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

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

2.4.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~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()メソッドによりdfhourを結合しています.

次に,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()メソッドによりdfdayを結合しています.

また,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()メソッドによりdfweekを結合しています.

最後に,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()メソッドによりdfmonthを結合しています.

また,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つ下にずらすことで上端で欠測が生じるため).

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

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

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

次に,作成されたデータフレーム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行目では,データフレーム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行~最終行までが検証用データとなります.

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

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

2.4.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に対してそれぞれ標準化を行います.

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

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

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

これで入力データの準備が完了しました.ここでは機械学習ライブラリ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は,乱数発生器の番号をしてします.重み付けやバイアスの初期値やバッチサンプリング(ソルバー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])

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

さて,結果を見てみましょう.ニューラルネットワークモデルの予測精度(決定係数$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節のセルまで戻る必要がありますので気をつけて下さい.

2.4.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ヶ月分のデータの範囲が指定されています.

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

ニューラルネットワークにより予測された電力消費量の時系列図の結果を見てみると,正解データと予測データは良い一致を示していることが分かります.「気温」,「曜日」,「時刻」というたった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から散布図を作成しています.

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

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

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

2.5 まとめ

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

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

2.6 著作権について

Copyright (c) 気象データ×IT勉強会,吉野 純,岐阜大学工学部附属応用気象研究センター 2022 All rights reserved.
<利用条件>
本書は、本書に記載した要件・技術・方式に関する内容が変更されないこと、および出典を明示いただくことを前提に、無償でその全部または一部を複製、翻案、翻訳、転記、引用、公衆送信等して利用できます。なお、全体または一部を複製、翻案、翻訳された場合は、本書にある著作権表示および利用条件を明示してください。
<免責事項>
本書の著作権者は、本書の記載内容に関して、その正確性、商品性、利用目的への適合性等に関して保証するものではなく、特許権、著作権、その他の権利を侵害していないことを保証するものでもありません。本書の利用により生じた損害について、本書の著作権者は、法律上のいかなる責任も負いません。