[chainer+optuna]ニューラルネットワークのハイパーパラメータを最適化して回帰モデルを実装

# １．はじめに
　初めまして。ijykと申します。
G検定2018#2に合格しましたが、実践的なことは何もできていないので、勉強もかねてpythonで色々と遊んでいます。

　ニューラルネットワーク(NN)を使った回帰(予測)した事例はいくつも見かけますが、
・chainerとoptunaを使った回帰モデルとハイパーパラメータのチューニング
・ハイパーパラメータ最適化後の回帰モデル実装、可視化までの流れ
はあまり見かけなかったので、今回、勉強も兼ねて実装したものをまとめました。
すみませんが、用語の説明などは省かせて頂きます。

# ２．やりたいこと
　1)ニューラルネットワークの回帰（予測）モデルを作りたい
　2)ハイパーパラメータを最適化したい
　3)最適化したハイパーパラメータで予測したい
　4)予測結果（精度とか誤差とか）を可視化したい

・・・てんこ盛り！

# ３．さあやってみよう
## 1)回帰（予測）モデルの作成
　今回は、シンプルにボストンの住宅価格のデータセットを使用します。
他のデータセットで実装したい場合は、X,Yのデータを変えれば使えます。

```python
from sklearn.datasets import load_boston
from sklearn.preprocessing import StandardScaler　#標準化に使います
import pandas as pd
import numpy as np

#boston住宅価格を読み込む
boston=load_boston()
X=boston['data']
Y=boston['target'].reshape(len(X),1)
x_col=boston['feature_names']

#データを標準化
x_scalar=StandardScaler()
X=x_scalar.fit_transform(X)
y_scalar=StandardScaler()
Y=y_scalar.fit_transform(Y)

#pandasに格納
X=pd.DataFrame(X)
Y=pd.DataFrame(Y)
Y.columns=['target']
X.columns=x_col
data_std=pd.concat([X,Y],axis=1)

#chainer用にfloat32型に変換
x=np.array(data_std[x_col]).astype('float32')
y=np.array(data_std['target']).astype('float32')
y_n=y.ndim

#目的変数の次元が1つだとreshapeが必要
y=y.reshape(len(y),y_n)
```

次に、chainerに放り込むデータセットをいじります。


```python
import chainer
import chainer.functions as F
import chainer.links as L
from chainer import report 
from chainer.datasets import split_dataset_random

#chainerに用いるデータセットをタプルに格納
dataset=chainer.datasets.TupleDataset(x,y)

#教師データとテストデータに分割
n_train=int(len(dataset)*0.7)
n_train_examples = n_train
n_test_examples = len(dataset) - n_train

#データセットの作成
rng = np.random.RandomState(0)
train, test = split_dataset_random(dataset, n_train, seed=0)
train = chainer.datasets.SubDataset(
    train, 0, n_train_examples, order=rng.permutation(len(train)))
test = chainer.datasets.SubDataset(
    test, 0, n_test_examples, order=rng.permutation(len(test)))

#データセット内の説明変数と目的変数を分ける
X_train, Y_train = map(list, zip(*train))
X_test, Y_test = map(list, zip(*test))
```
chainerを使ってモデルを作っていきますが、scikit-learnのMLPRegressorとは違い、少しコードが複雑なので素人の私には理解するのに時間がかかりました。defは未だに分からん...
今回は、optunaを使ってハイパーパラメータを最適化するので、最適化したいパラメータは変数化します。
〇〇=trial.suggest_△△（◇◇）　とすればOK。
最適化するパラメータは
　・隠れ層の数
　・各層のニューロン数
　・バッチ数
　・活性化関数
　・エポック数
　・最適化関数
　・重み減衰
です。少なくとも上3つは寄与が大きいのでしたほうが良いそうです。とりあえず色々と変数化してみました。


```python
#ニューラルネットワークの構造
def MLP(trial):
#活性化関数と隠れ層の数をoptuna用に変数化
    activation=trial.suggest_categorical('activation', ['relu', 'tanh','sigmoid'])
    n_fc_layers = trial.suggest_int('n_fc_layers', 1, 5)
    if activation=='relu':
        act=F.relu
    elif activation=='tanh':
        act=F.tanh
    else:
        act=F.sigmoid
#設定した隠れ層の数に基づきネットワークを作成
    layers = []
    for i in range(n_fc_layers):
#ニューロン数も変数化
        n_units = int(trial.suggest_int('n_units_l{}'.format(i), 4, 128))
        layers.append(L.Linear(None, n_units))
        layers.append(act)
    layers.append(L.Linear(None, y_n))#y_n･･･目的変数の次元
    return chainer.Sequential(*layers)
```
分類はL.Classifierですが、回帰は別の関数がありますのでそちらを使います。
さらに、最適化関数も選べるように設定します。

```python
class MyRegressor(chainer.Chain):

    def __init__(self, predictor):
        super(MyRegressor, self).__init__(predictor=predictor)

    def __call__(self, x, y):
#予測するときは〇〇.predictor(x)使用
#平均絶対誤差(MAE)と、平均二乗誤差(MSE)を返す
        pred = self.predictor(x)
        abs_error = F.sum(abs(pred - y)) / len(x.data)
        loss = F.mean_squared_error(pred, y)

#平均絶対誤差(MAE)と平均二乗誤差(MSE)を記録
        report({'abs_error': abs_error, 'squared_error': loss}, self)
        return loss

#最適化関数を決める関数
def create_optimizer(trial, model):
# 最適化関数の選択
    optimizer_name = trial.suggest_categorical('optimizer', ['Adam', 'MomentumSGD'])
    if optimizer_name == 'Adam':
        adam_alpha = trial.suggest_loguniform('adam_alpha', 1e-5, 1e-1)
        optimizer = chainer.optimizers.Adam(alpha=adam_alpha)
    else:
        momentum_sgd_lr = trial.suggest_loguniform('momentum_sgd_lr', 1e-5, 1e-1)
        optimizer = chainer.optimizers.MomentumSGD(lr=momentum_sgd_lr)

    weight_decay = trial.suggest_loguniform('weight_decay', 1e-10, 1e-3)
    optimizer.setup(model)
    optimizer.add_hook(chainer.optimizer.WeightDecay(weight_decay))
    return optimizer
```
ハイパーパラメータを最適化するときはobjectiveという関数を作ります。基本的にはここに最適化したいパラメータや、出力する誤差、学習実行のコードを組み込んでおきます。

```python
def objective(trial):
# モデルのインスタンス化
    model = MyRegressor(MLP(trial))
    optimizer = create_optimizer(trial, model) # モデルとoptimizerを紐付ける
    batchsize=trial.suggest_int('batchsize',10,len(y))
    epoch=trial.suggest_int('epoch',10,50)
# Iteratorの設定

    train_iter = chainer.iterators.SerialIterator(train, batchsize)
    test_iter = chainer.iterators.SerialIterator(test, batchsize, repeat=False, shuffle=False)
# Trainerの設定
    updater = chainer.training.StandardUpdater(train_iter, optimizer)
    trainer = chainer.training.Trainer(updater, (epoch, 'epoch'))
    trainer.extend(chainer.training.extensions.Evaluator(test_iter, model))
    log_report_extension = chainer.training.extensions.LogReport(log_name=None)
    trainer.extend(chainer.training.extensions.PrintReport(
        ['epoch', 'main/squared_error', 'validation/main/squared_error',
         'main/abs_error', 'validation/main/abs_error', 'elapsed_time']))
    trainer.extend(log_report_extension)

# 学習の実行
    trainer.run()

# 学習結果の保存
    log_last = log_report_extension.log[-1]
    for key, value in log_last.items():
        trial.set_user_attr(key, value)

# 最終的なバリデーションの値を返す
    val_err = log_report_extension.log[-1]['validation/main/squared_error']
    return val_err
```

## 2)ハイパーパラメータ最適化
ハイパーパラメータの最適化手法はグリッドサーチやランダムサーチなどがありますが、今回はoptunaというベイズ最適化ライブラリを使って最適化します。ベイズ最適化をすごくざっくり要約すると、誤差が小さくなりそうなパラメータを人のカンコツのように「こっちのほうが良くなりそうだ」といった挙動で探してくれる優れものです。厳密にいうと、確率論的にパラメータを選択していくわけですが、ここでは割愛します(詳しく知りたい方は最下部参考⑩などを読んでみてください)。
　ハイパーパラメータ最適化は学習を都度繰り返すので、途中で精度がサチュレートしたときに打ち切りたいときがあります。optunaではそれが可能なので必要に応じて使用します。

```python
import optuna

#最適化
study = optuna.create_study()

#枝刈りを行う場合は以下のコメントアウトを外して使用
#study = optuna.create_study(pruner=optuna.pruners.MedianPruner(n_warmup_steps=5))

study.optimize(objective, n_trials=50)　　　#n_trials･･･探索回数

#最適化されたパラメータの出力
print('Number of finished trials: ', len(study.trials))

print('Best trial: ', )
trial = study.best_trial

#ハイパーパラメータの出力と同時に呼び出せるように格納
print('Params: ')
param_n=[]
param_v=[]
for key, value in trial.params.items():
    print('{}:{}'.format(key, value))
    param_n.append(key)
    param_v.append(value)
    
print('User attrs: ')
attrs_n=[]
attrs_v=[]
for key, value in trial.user_attrs.items():
    print('{}:{}'.format(key, value))
    attrs_n.append(key)
    attrs_v.append(value)
```
変数化したものはすべて出力されます。下は計算結果です。毎回結果が変わり、同じハイパーパラメータは得られないのであしからず。

| Params |値　|
|:--:|:-|
|activation|relu|
|n_fc_layers|3|
|n_units_l0|57|
|n_units_l1|127|
|n_units_l2|54|
|optimizer|Adam|
|adam_alpha|1.634e-2|
|weight_decay|9.576e-8|
|batchsize|52|
|epoch|39|

| User attrs |値　|
|:--:|:-|
|main/abs_error|1.086e-1|
|main/squared_error|2.251e-2|
|validation/main/abs_error|1.907e-1|
|validation/main/squared_error|6.801e-2|
|epoch|39|
|iteration|266|
|elapsed_time|3.085|

## 3)最適化したハイパーパラメータで予測
最適化したハイパーパラメータが得られたので、NNで予測してみます。sklearnとchainerで実装しました。
### (1)sklearn(MLPRegressor)での実装
```python
from sklearn.neural_network import MLPRegressor
from sklearn.model_selection import train_test_split

#NN回帰モデル
model=MLPRegressor()
#sklearnは隠れ層が複数の場合タプルにしないといけないので
#タプルに出力させる
n_units=()
for i in range(param_v[1]):
    n_units=n_units+(param_v[i+2],)#「,」がないとエラーが出る
#最適化関数別にパラメータを設定
if param_v[param_v[1]+2]=='Adam':
    model.set_params(activation=param_v[0],alpha=param_v[param_v[1]+3],
                     batch_size=int(param_v[param_v[1]+5]),hidden_layer_sizes=n_units,
                     max_iter=int(attrs_v[5]),solver='adam')
else:
    model.set_params(activation=param_v[0],learning_rate_init=param_v[param_v[1]+3],
                     batch_size=int(param_v[param_v[1]+5]),hidden_layer_sizes=n_units,
                     max_iter=int(attrs_v[5]),solver='sgd')
#モデルの学習
model.fit(X_train,Y_train)

#予測データ、元データの標準化を戻す。絶対誤差も計算。
Y_test_in=y_scalar.inverse_transform(Y_test).transpose()
X_test_pre=y_scalar.inverse_transform(model.predict(X_test))
error_test=Y_test_in-X_test_pre
Y_train_in=y_scalar.inverse_transform(Y_train).transpose()
X_train_pre=y_scalar.inverse_transform(model.predict(X_train))
error_train=Y_train_in-X_train_pre


#予測結果
print('予測精度(テスト):{}'.format(model.score(X_test,Y_test)))
print('予測精度(教師):{}'.format(model.score(X_train,Y_train)))
```
予測精度(テスト):0.864
予測精度(教師):0.977

まずまずといったところでしょうか。
過学習気味なので、エポック数などを調整するなどの工夫が必要そうです。

### (2)chainerでの予測

```python
#ハイパーパラメータを呼び出す
activation=param_v[0]
batchsize=int(param_v[param_v[1]+5])
epoch=int(attrs_v[4])
optimizer_name=param_v[param_v[1]+2]
weight_decay=param_v[param_v[1]+4]
n_fc_layers=param_v[1]

#最適化関数別にパラメータを設定
if param_v[param_v[1]+2]=='Adam':
    adam_alpha=param_v[param_v[1]+3]
else:
    momentum_sgd_lr=param_v[param_v[1]+3]
```
ここから先は最適化時のコードをいじっただけです。
詳細は省きます。

```python
def MLP():
    if activation=='relu':
        act=F.relu
    elif activation=='tanh':
        act=F.tanh
    layers = []
    for i in range(n_fc_layers):
        n_units = param_v[i+2]
       #dropout_ratio= trial.suggest_loguniform('dropout_r{}'.format(i), 0.5, 1)
        layers.append(L.Linear(None, n_units))
        layers.append(act)
        #layers.append(F.dropout(act,ratio=dropout_ratio))
    layers.append(L.Linear(None, y_n))
    return chainer.Sequential(*layers)

class MyRegressor(chainer.Chain):

    def __init__(self, predictor):
        super(MyRegressor, self).__init__(predictor=predictor)

    def __call__(self, x, y):
        pred = self.predictor(x)
        abs_error = F.sum(abs(pred - y)) / len(x.data)
        loss = F.mean_squared_error(pred, y)
        report({'abs_error': abs_error, 'squared_error': loss}, self)
        return loss

def create_optimizer( model):
    if optimizer_name == 'Adam':
        optimizer = chainer.optimizers.Adam(alpha=adam_alpha)
    else:
        optimizer = chainer.optimizers.MomentumSGD(lr=momentum_sgd_lr)
    optimizer.setup(model)
    optimizer.add_hook(chainer.optimizer.WeightDecay(weight_decay))
    return optimizer

model = MyRegressor(MLP())
optimizer = create_optimizer(model) 
train_iter = chainer.iterators.SerialIterator(train, batchsize)
test_iter = chainer.iterators.SerialIterator(test, batchsize, repeat=False, shuffle=False)

updater = chainer.training.StandardUpdater(train_iter, optimizer)
trainer = chainer.training.Trainer(updater, (epoch, 'epoch'))
trainer.extend(chainer.training.extensions.Evaluator(test_iter, model))
log_report_extension = chainer.training.extensions.LogReport(log_name=None)
trainer.extend(chainer.training.extensions.PrintReport(
    ['epoch', 'main/squared_error', 'validation/main/squared_error',
    'main/abs_error', 'validation/main/abs_error', 'elapsed_time']))
trainer.extend(log_report_extension)

#学習の実行
trainer.run()

#予測データ、元データの標準化を戻します。絶対誤差も計算します。
#chainerの予測データはVariable形で出力されるので、「.array」でnumpyに変換しています
test_pre=y_scalar.inverse_transform(model.predictor(np.array(X_test).astype('float32')).array)
test_Y=y_scalar.inverse_transform(np.array(Y_test).astype('float32'))
train_pre=y_scalar.inverse_transform(model.predictor(np.array(X_train).astype('float32')).array)
train_Y=y_scalar.inverse_transform(np.array(Y_train).astype('float32'))
error_te=test_pre-test_Y
error_tr=train_pre-train_Y

#予測精度
test_sc=F.r2_score(test_pre,test_Y)
train_sc=F.r2_score(train_pre,train_Y)
print('予測精度(テスト):{}'.format(train_sc.array))
print('予測精度(教師):{}'.format(test_sc.array))
```
予測精度(テスト):0.955
予測精度(教師):0.955

chainerはかなり良い結果になりました。エポック数や重み減衰のパラメータが設定できる分精度がよくなっているのでしょうか。
## 4)予測結果の可視化

　matplotlibを使って可視化します。可視化するとどこが外れ値とか、本当に精度よく予測できてるかよくわかります。
### (1)sklearn
```python
#予測結果の可視化
import matplotlib.pyplot as plt
plt.scatter(Y_train_in,X_train_pre)
plt.scatter(Y_test_in,X_test_pre)
plt.xlabel("target")
plt.ylabel("predict")
plt.show()

#誤差の可視化
plt.scatter(Y_train_in,error_train)
plt.scatter(Y_test_in,error_test)
plt.ylim(-20,20)
plt.xlabel("target")
plt.ylabel("error")
plt.show()
```
青色が教師データ、オレンジ色がテストデータです。
#### 元データvs予測データ
![あ.png](https://qiita-image-store.s3.amazonaws.com/0/350979/80386f53-287c-bee8-f652-9826e0bac0ba.png)
#### 元データvs絶対誤差
![い.png](https://qiita-image-store.s3.amazonaws.com/0/350979/4f950e73-b3ae-5a3a-d140-b070b69a0457.png)

いくつか外れ値が見られます。右端の値については大きく外れていますね。やや過学習気味にも見えます。。。
### (2)chainer
```python
plt.scatter(train_Y,train_pre)
plt.scatter(test_Y,test_pre)
plt.xlabel("target")
plt.ylabel("predict")
plt.show()

plt.scatter(train_Y,error_tr)
plt.scatter(test_Y,error_te)
plt.ylim(-20,20)
plt.xlabel("target")
plt.ylabel("error")
plt.show()
```
同じく青色が教師データ、オレンジ色がテストデータです。
#### 元データvs予測データ
![image.png](https://qiita-image-store.s3.amazonaws.com/0/350979/daec2972-5cbe-1e13-281f-d0687f7da601.png)
#### 元データvs絶対誤差
![う.png](https://qiita-image-store.s3.amazonaws.com/0/350979/bd24caf7-a88c-5be5-5243-ed406ad9bfc9.png)
かなりリニアになりましたが右端のデータがばらついています。それでもまずまずの精度はでてそうです。

# ４．まとめ
　ディープラーニングを行う上では、一番基本となる教師あり学習について記載しました。ちなみに、データの前処理に標準化を行っていますが、標準化しなかった場合予測精度はR2=0.5以下とかなり悪かったです。今回は、ハイパーパラメータの最適化に注力しましたが、データセットの前処理（クレンジング）も必要不可欠ということですね。今後はそのあたりも含めて検討したいです。
　大真面目に色々書きましたが、データセットなどユニークなものでも検討してみたいですね。

以下は参考にさせて頂いたページとライブラリのURLです。

①chainerのTrainerやらUpdaterやらの仕組みを理解したかった
　https://qiita.com/sumsum88/items/a62b6950533e8dbb7e02
②ChainerとOptunaでハイパーパラメータ探索
　https://qiita.com/yoshiyoshi0505/items/2e39ffdab305dba95164
③chainer
　https://chainer.org/
④chainerのMyRegressorについて
　https://docs.chainer.org/en/stable/reference/util/generated/chainer.report.html
⑤optuna
　https://optuna.org/
⑥numpy
　http://www.numpy.org/
⑦pandas
　https://pandas.pydata.org/
⑧scikit-learn
　https://scikit-learn.org/stable/
⑨matplotlib
　https://matplotlib.org/
⑩ベイズ最適化入門
　https://qiita.com/masasora/items/cc2f10cb79f8c0a6bbaa

投稿 '19/2/5
改訂 '22/4/12