ニューラルネットワークによる回帰分析をしてみよう
ここでは重回帰分析に失敗するようなデータについて，ニューラルネットワーク（ディープラーニング／深層学習）による回帰分析を行って推定や予測を行います．このデータの詳細はここを参照してください．なお，多くの参考書やWebの記事ではニューラルネットワーク（ディープラーニング／深層学習）による分類（文字認識や画像認識など）が説明されています．これらと同じような考え方で回帰分析も可能ですが多くの記事ではあまり回帰モデルに触れていないようです．
データの読み込み
まずは必要なモジュールをインポートします．なお，ニューラルネットワーク（ディープラーニング／深層学習）を行うためのモジュールは様々ありますが，ここでは Keras を使うことにします．したがって，このあたりを参考に，pip install keras コマンドで Keras を（同時に TensorFlow も pip install tensorflow コマンドで）インストールしてください．
モジュールのインポート
import pandas as pd
import numpy as np
import seaborn as sns
import matplotlib.pyplot as plt
from keras.models import Sequential
from keras.layers import Dense
# 高解像度ディスプレイ用
from IPython.display import set_matplotlib_formats
# from matplotlib_inline.backend_inline import set_matplotlib_formats # バージョンによってはこちらを有効に
set_matplotlib_formats('retina')
次に CSV ファイルを読み込んで表示してみます．このデータは，\(x_1\), \(x_2\), \(x_3\) の値が決まれば何らかの法則で \(y\) の値がおおよそ決まるような500件のデータです．（実際にはここに示すような恣意的なデータです．）
ファイルを読み込んで表示する
url = "https://github.com/rinsaka/sample-data-sets/blob/master/mra-02.csv?raw=true"
df = pd.read_csv(url)
print(df)
no x1 x2 x3 y 0 0 0.412 3.968 3.083 -0.764 1 1 2.562 4.609 0.683 -0.648 2 2 0.874 4.577 2.704 1.005 3 3 0.869 3.275 4.463 0.200 4 4 1.326 4.764 4.729 0.212 .. ... ... ... ... ... 495 495 3.552 4.549 1.786 -0.806 496 496 4.754 1.685 1.406 1.019 497 497 0.826 4.907 3.840 0.641 498 498 1.948 3.455 1.806 1.088 499 499 2.671 2.491 4.405 -0.398 [500 rows x 5 columns]
相関係数と散布図
ニューラルネットワークによる回帰分析を行う前に，変数間の相関係数を計算し，散布図も描いて確認してみよう．まずは，ここと同じ方法で相関係数を計算します．これには Pandas のデータフレームを NumPy 配列に変換し，これを転置して np.corrcoef 関数に与えると良いでしょう．
相関係数
# x1, x2, x3, y 列だけとりだして NumPy 配列に変換
xy = df.loc[:,['x1','x2','x3','y']].values
# NumPy配列を転置して相関係数を求める
print(np.corrcoef(xy.T))
[[ 1. 0.02061493 -0.0424868 0.03623328] [ 0.02061493 1. 0.00665035 0.00484668] [-0.0424868 0.00665035 1. 0.06805707] [ 0.03623328 0.00484668 0.06805707 1. ]]
上の結果から，どの変数間にも相関関係がなさそうであることがわかります．
次に，変数ごとの散布図を描いて変数間の関連を確認します．複数の散布図をまとめて描くには seaborn パッケージを使うと良いでしょう．このとき Pandas のデータフレームには no 列が含まれているので，この列を削除してから散布図を描きます．なお，df.drop でデータフレームの行や列を削除できるが，axis=1とすると列を削除し，axis=0や省略すると行を削除します．
データフレームから散布図を描く
# no列を削除したデータフレームについて， x1, x2, x3, y それぞれの組み合わせで散布図を描く
sns.pairplot(df.drop('no', axis=1))
plt.show()
上の散布図を確認すると，やはり変数間の相関関係はなさそうです．
ニューラルネットワークによる回帰分析
データの準備
まず，500件のデータをトレーニングデータ (80%) と テストデータ (20%) に分割します．ここで説明したような方法でシャッフルすることも可能ですが，同じデータが重複して各データに含まれる可能性があるので，別の方法で行います．なお，以下のコードはよりスッキリと記述できるはずです．
まず Pandas のデータフレームを NumPy 配列（行列）に変換します．なお，no の列は本来不要（さらに読み込むとしても整数にすべき）ですが，あとで行列をシャッフルしたときに正しくシャッフルできているかを確認するためだけにここでは追加しています．
xy_data = df.loc[:,['no', 'x1', 'x2', 'x3', 'y']].values
print(xy_data)
[[ 0.000e+00 4.120e-01 3.968e+00 3.083e+00 -7.640e-01] [ 1.000e+00 2.562e+00 4.609e+00 6.830e-01 -6.480e-01] [ 2.000e+00 8.740e-01 4.577e+00 2.704e+00 1.005e+00] ... [ 4.970e+02 8.260e-01 4.907e+00 3.840e+00 6.410e-01] [ 4.980e+02 1.948e+00 3.455e+00 1.806e+00 1.088e+00] [ 4.990e+02 2.671e+00 2.491e+00 4.405e+00 -3.980e-01]]
行列をシャッフルする前に乱数生成器を初期化しておきます．
rng = np.random.default_rng(1)
# rng = np.random.default_rng()
行列をランダムに並べ替え（シャッフル）します．念の為に，no 列（先頭列）が重複していないことを確認すると良いでしょう．
rng.shuffle(xy_data)
print(xy_data)
[[ 2.750e+02 1.318e+00 4.880e+00 4.960e+00 2.160e-01] [ 2.810e+02 2.219e+00 4.486e+00 7.590e-01 -8.920e-01] [ 1.700e+02 2.712e+00 3.654e+00 2.913e+00 1.169e+00] ... [ 1.870e+02 6.730e-01 2.275e+00 4.701e+00 6.690e-01] [ 1.020e+02 2.700e-01 1.922e+00 3.633e+00 6.850e-01] [ 2.550e+02 4.817e+00 3.011e+00 4.381e+00 1.062e+00]]
しきい値 (threshold) をリストの長さの 80% に設定し，トレーニングデータに格納するデータのインデックスと，テストデータに格納するデータのインデックスを生成します．このとき，np.round() 関数で丸め込んでもなお実数型 (float64) なので，スライスするとエラーになります．
threshold = np.round(500*.8)
print(threshold)
print(threshold.dtype)
400.0 float64
エラーなく分割するには次のようにしきい値を整数化しておく必要があります．
threshold = np.round(500*.8).astype(np.int64)
print(threshold)
print(threshold.dtype)
400 int64
しきい値 threshold を使ってデータを2つに分割します．
xy_train = xy_data[0:threshold]
xy_test = xy_data[threshold:]
print(xy_train.shape)
print(xy_test.shape)
(400, 5) (100, 5)
さらに分割したデータを x のデータ行列と y のデータ配列に分割します．このとき，no 列は含めないように注意します．
x_train = xy_train[:, 1:4]
y_train = xy_train[:, 4]
x_test = xy_test[:, 1:4]
y_test = xy_test[:,4]
トレーニングデータを表示し，そのサイズを確認します．
print(x_train)
print(y_train)
print(x_train.shape)
print(y_train.shape)
[[1.318 4.88 4.96 ] [2.219 4.486 0.759] [2.712 3.654 2.913] ... [1.676 1.491 2.421] [3.358 2.624 1.849] [0.054 4.307 2.63 ]] [ 0.216 -0.892 1.169 -0.101 -0.11 1.126 1.02 -0.144 -0.791 0.978 ．．．（中略）．．． 0.015 -0.458 -0.649 -0.598 0.888 -0.822 0.774 0.734 1.015 -0.139] (400, 3) (400,)
同様にテストデータも表示し，そのサイズとともに確認します．
print(x_test)
print(y_test)
print(x_test.shape)
print(y_test.shape)
[[4.05 4.195 0.209] [2.037 2.504 3.917] ．．．（中略） [4.817 3.011 4.381]] [ 0.648 -0.222 1.067 1.144 -0.373 -0.712 -0.807 0.798 -0.922 0.938 ．．．（中略） 0.995 -0.391 0.046 -0.938 0.389 -0.405 -0.756 0.669 0.685 1.062] (100, 3) (100,)
モデルの作成
データの準備ができたので，ここでは多層ニューラルネットワークのモデルを Keras で作成します．今回作成するモデルは次の構造です．
- 入力層
- ニューロン数 = 3
- 中間層1
- ニューロン数 = 64
- 中間層2
- ニューロン数 = 32
- 出力層
- ニューロン数 = 1
ここで，今回のデータは x1, x2, x3 からなる3次元のデータであるので，入力層のニューロン数が3です．また文字認識のような分類モデルでは分類したいクラス分だけ出力層にニューロンを設置するが，ここでの回帰モデルでは y の値を予測したいので，ニューロン数が1です．さらに，とりあえず2つの中間層を置き，中間層1，中間層2のニューロン数をそれぞれ 64，32 とします．さらに中間層の活性化関数には ReLU 関数を，出力層の活性化関数には恒等関数 (linear) を指定します．（中間層の数やニューロン数などにはまだ検討の余地があります．）
では，このモデルを Keras によって記述します．
nn1 = 64
nn2 = 32
model = Sequential()
model.add(Dense(nn1, activation='relu', input_dim=3))
model.add(Dense(nn2, activation='relu'))
model.add(Dense(1, activation='linear'))
次に，モデルをコンパイルします．なお，損失関数には平均2乗誤差 (mean_squared_error) を使い，最適化アルゴリズムには RMSProp オプティマイザを指定します．（ここも検討の余地は十分あります）．
model.compile(optimizer='rmsprop',
loss='mean_squared_error',
metrics=['accuracy'])
モデルのサマリーを確認します．2369個のパラメータを最適化することになります．このパラメータ数が2369個になる理由はここを参照してください．
model.summary()
Model: "sequential" _________________________________________________________________ Layer (type) Output Shape Param # ================================================================= dense (Dense) (None, 64) 256 _________________________________________________________________ dense_1 (Dense) (None, 32) 2080 _________________________________________________________________ dense_2 (Dense) (None, 1) 33 ================================================================= Total params: 2,369 Trainable params: 2,369 Non-trainable params: 0 _________________________________________________________________
モデルの学習
実際にトレーニングデータを用いてモデルを学習させる．ここでミニバッチのサイズを50，エポックを3000にしているが，このあたりも検討が必要です．実行したときの途中経過は train_history に格納されます．学習の結果，損失関数の値 (loss) が徐々に小さくなり，収束していく様子が確認できます．
train_history = model.fit(x_train, y_train,
batch_size=50,
epochs=3000,
verbose=1)
Epoch 1/3000 8/8 [==============================] - 0s 984us/step - loss: 0.5675 - accuracy: 0.0000e+00 Epoch 2/3000 8/8 [==============================] - 0s 1ms/step - loss: 0.5208 - accuracy: 0.0000e+00 Epoch 3/3000 8/8 [==============================] - 0s 995us/step - loss: 0.5023 - accuracy: 0.0000e+00 Epoch 4/3000 8/8 [==============================] - 0s 1ms/step - loss: 0.4928 - accuracy: 0.0000e+00 ．．．（中略）．．． Epoch 2995/3000 8/8 [==============================] - 0s 743us/step - loss: 0.0017 - accuracy: 0.0000e+00 Epoch 2996/3000 8/8 [==============================] - 0s 641us/step - loss: 0.0053 - accuracy: 0.0000e+00 Epoch 2997/3000 8/8 [==============================] - 0s 760us/step - loss: 0.0032 - accuracy: 0.0000e+00 Epoch 2998/3000 8/8 [==============================] - 0s 733us/step - loss: 0.0032 - accuracy: 0.0000e+00 Epoch 2999/3000 8/8 [==============================] - 0s 703us/step - loss: 0.0036 - accuracy: 0.0000e+00 Epoch 3000/3000 8/8 [==============================] - 0s 716us/step - loss: 0.0033 - accuracy: 0.0000e+00
学習経過は train_history に格納されています．次のようにすると，train_history.history として辞書形式で保存され，'loss' と 'accuracy' のキーがあることがわかります．
train_history.history.keys()
dict_keys(['loss', 'accuracy'])
学習経過を表示してみます．
print(train_history.history)
{'loss':
[0.5674909353256226, 0.5208286643028259, 0.5022525191307068,
．．．（中略）．．
0.003228976158425212, 0.003617800772190094, 0.003272315254434943],
'accuracy':
[0.0, 0.0, 0.0, ．．．（中略）．．． 0.0, 0.0, 0.0]}
学習は3000エポックまで行ったので，結果は3000あります．
len(train_history.history['loss'])
3000
上の学習過程でエポックごとの損失関数値をプロットしてみましょう．3000エポックまで学習しましたが，1500エポック程度で収束していたようです．
fig, ax = plt.subplots(1, 1, figsize=(6, 4))
ax.plot(train_history.history['loss'])
ax.set_xlabel('Epoch')
ax.set_ylabel('loss')
ax.set_xlim(0, 3000)
ax.set_ylim(-0.01, 0.5)
# plt.savefig('nn-regression-train-1.png', dpi=300, facecolor='white')
plt.show()
200エポックまでの拡大図です．この段階ではまだ収束はしていません．
fig, ax = plt.subplots(1, 1, figsize=(6, 4))
ax.plot(train_history.history['loss'][:200])
ax.set_xlabel('Epoch')
ax.set_ylabel('loss')
ax.set_xlim(0, 200)
ax.set_ylim(-0.01, 0.5)
# plt.savefig('nn-regression-train-2.png', dpi=300, facecolor='white')
plt.show()
2000エポックまでの拡大図です．概ね収束しています．
fig, ax = plt.subplots(1, 1, figsize=(6, 4))
ax.plot(train_history.history['loss'][:2000])
ax.set_xlabel('Epoch')
ax.set_ylabel('loss')
ax.set_xlim(0, 2000)
ax.set_ylim(-0.01, 0.2)
# plt.savefig('nn-regression-train-3.png', dpi=300, facecolor='white')
plt.show()
なお，ここで説明したような方法を使う場合は，学習データに過剰に適合していないか（つまり過学習していないか）を検証する必要があります．
予測とその評価
3000エポックまで学習したモデルに学習に利用しなかった100件のテストデータを与えて予測を行います．
y_pred = model.predict(x_test)
print(y_pred)
[[ 0.29134643] [-0.30296147] [ 1.022902 ] ．．．（中略）．．． [ 0.6234139 ] [ 0.61915994] [ 0.9628694 ]]
次に，ある2組のテストデータに対する正解値と，予測値を並べて表示してみます．
print(y_test[0], y_pred[0][0])
print(y_test[1], y_pred[1][0])
0.648 0.29134643 -0.222 -0.30296147
正解値と予測値の散布図を作成します．この図から，重回帰分析では全く歯が立たなかったデータに対して，かなり正確に予測できていることがわかりました．
fig, ax = plt.subplots(1, 1, figsize=(6, 4))
ax.scatter(y_test, y_pred, alpha=0.5)
ax.set_xlabel('y_data')
ax.set_ylabel('y_pred')
ax.set_xlim(-1.1,1.3)
ax.set_ylim(-1.1,1.3)
# plt.savefig('nn-regression-predict.png', dpi=300, facecolor='white')
plt.show()
