Python Keras(TensorFlow)によるニューラルネットワーク回帰分析

　TensorFlow2.0が'19/10/01に正式リリースされてKeras（ケラス）が統合された。これにより、Sequential（シーケンシャル）にNeural Networkを構築する方法が主流になりそうに思う。それは順番に層を積み重ねて記述する方法であり、define by runでデバックが容易とされる。以前のバージョンの主流はテンソルフローグラフのセッション形式であり、define and runと呼び、一旦一通りモデルを作成する必要があった。この背景にはPyTorchの躍進がある模様→PyTorch vs. TensorFlow、ディープラーニングフレームワークはどっちを使うべきか問題 (1/2)：気になるニュース＆ネット記事 - ＠IT

　本記事では、シーケンシャルに記述する方法に加えて、Kerasの6つの最適化アルゴリズム（オプティマイザ）で実行して各々を比較できる。また、学習したモデルをオブジェクトファイル「.h5」で保存するまでを行うコードを記載した。

●本プログラムでは、下図のように「zzz_create_NNmodel.py」適当な名で保存して実行するだけで試せます。分析データはsklearnに含まれているボストンデータセット(Boston Housing dataset)を利用しています。

●下図は実行後の結果です。学習履歴などをcsvファイルや画像グラフに保存するだけでなく、学習した結果の学習器モデル（拡張子「.h5」）も保存する仕様です。

▼世代が進行するにつれて、損失関数が収束する様子をグラフ化

　下図の縦軸は損失関数として平均二乗誤差 (Mean Squared Error, MSE)を選択した場合である。値が小さい程、精度が高いことを意味する。横軸は繰り返し数(iterations, epoch)などと呼ばれる。kerasのEarlyStoppingを使用して、損失関数が収束してサチッた場合に途中でストップする枝刈り機能を利用している。

▼オプティマイザ間の損失関数や評価関数、解析時間の比較結果をcsvファイル出力

詳細は本コードを見て頂ければわかるが、tensorflowのmodel.compileにおいて、「loss」に設定する評価指標を収束判定として使用する。これを損失関数と呼ぶ。一方の「metrics」に設定する評価指標は観察するだけであり、評価関数と呼ぶ。

▼オプティマイザ別に損失関数や評価関数の変化の過程をcsvファイル出力

▼オプティマイザ別に予測結果をcsvファイル出力

▼モデル構成を可視化

多層パーセプトロン（Multilayer perceptron）である。上から順に、入力層（特徴量が13種類）→隠れ層1（ニューロン数1000）→隠れ層2（ニューロン数800）：ドロップアウト判定→隠れ層3（ニューロン数100）→出力層（1つ。本例題はボストン住宅価格のため）

▼実測値と予測の比較　その1

青色が実測値で赤色が予測結果である。オプティマイザ別に適合度合いを視覚的に確認できる。但し、各々のハイパーパラメータを合わせた結果ではない。

▼実測値と予測の比較　その2

横軸が実測で縦軸が予測結果である。比例y=x線上に近ければ傾向の再現性が比較的高いと言える。

▼途中の世代の隠れ層の重み等のモデル状態のデータを保存する

下図は、世代進行30ステップ毎にオプティマイザ別に出力した結果である。読み出し方法は、別途記事を書く予定。

■本プログラム
　実行にあたっては、Anaconda環境であれば本コードをコピペするだけで可能なはずとは思います。
以下に、特記事項を記載します。

1. 学習過程の枝刈り
callbackのEarlyStoppingを使用する。monitor に設定した損失関数の値が、 patienceの値の回数連続でmin_deltaの値以上に改善しない場合に学習を止める。

2. 評価関数metricsについて
model.compileの引数である評価関数metricsは、損失関数lossと異なります。評価結果の値が訓練に直接使われることはありません。

3. plotグラフのアスペクト比を正方形に保つ
plt.gca().set_aspect('equal', adjustable='box')で維持できます。

4. 散布図に式を記述する方法
seabornとscipyを使用しています。

#!/usr/bin/env python3
# -*- coding: utf-8 -*-
import os, sys, glob
import pandas as pd
import numpy as np
import scipy as sp
from scipy import optimize
import math
import tensorflow as tf
from tensorflow.keras.models import Sequential
from tensorflow.keras.layers import Activation, Dropout, Flatten, Dense
from tensorflow.keras.callbacks import Callback
from tensorflow.keras.callbacks import EarlyStopping
from tensorflow.keras.callbacks import ModelCheckpoint
from tensorflow.keras.callbacks import CSVLogger
from tensorflow.keras.utils import plot_model

from sklearn.model_selection import train_test_split
from sklearn import datasets
import time

import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt
import seaborn as sns

######################################## csvファイルの任意のデータセットをロード
######################################## 対応書式は1行目に列名、それ以下はデータとなってるもの
def my_load_csv(path, check_column, check_duplicate_columns,
                start_row, input_columns, out_start_column,
                target_Yname, target_value2ln):
    df = pd.read_csv(path, skiprows=start_row-1, encoding="utf-8")
    print("df -> " + str(df))
    print(df.dtypes) # pandasの列毎の型確認。その列に文字列があるかで決まる。但し、object型は要素毎に異なる可能性がある。
    df.replace('none', '0', inplace=True) # 置換。 inplace=Trueでdfが置換される(df=と新規に代入しなくてよい）
    df.replace('FALSE', '0', inplace=True)
    df_buf1 = df.iloc[:, :out_start_column]
    df_buf2 = df.iloc[:, out_start_column:] # ilocとreplaceは同時に使用できなかったため、一旦分割して処理する
    
    df_buf2.replace('0', np.nan, inplace=True) #文字列0をnumpy形式のNaNにする
    df_buf2.replace(0, np.nan, inplace=True) #数値0をnumpy形式のNaNにする    
    
    columns_list = df_buf2.columns
    #df_buf2[columns_list] = pd.to_numeric(df_buf2[columns_list], errors='coerce') #文字列をNaNへ置換する
    
    df = pd.concat([df_buf1, df_buf2], axis=1) #データフレームの列を連結して、ひとつのdfへ戻す
    df.to_csv("data_nan.csv", sep=",", header=True, index=False)
    print('check NaN -> ' + str(df.isnull().any())) # ひとつでもNaNがある列はTrue、そうでなければFalseと表示される。
    print(df.dtypes)
    
    if check_duplicate_columns: #指定した列名があった場合、重複行を削除する
        df2 = df.drop_duplicates(check_duplicate_columns, keep='first') #最初の行が重複とみなされずに残す
    else:
        df2 = df
    print("df2 -> " + str(df2))
    
    if check_column[0]:
        df3 = df2.dropna(subset=[check_column[1]]) # NaNやNoneなどの欠損値を除外する場合
    else:
        df3 = df2[df2[check_column[1]].str.contains(check_column[2], na=False)] # 指定文字列の行を残す場合
    print("df3 -> " + str(df3))
    
    if target_value2ln: # 目的変数Yを対数に変換する場合
        df4 = df3.assign(tmp = np.log(df3[target_Yname]))
        target_Yname = target_Yname + '_ln'
        tmp_columns = df4.columns.values
        print("tmp_columns -> " + str(tmp_columns))
        tmp_columns[-1] = target_Yname
        df4.columns = tmp_columns
        print("df4 -> " + str(df4))
        df4.to_csv("data_df4.csv", sep=",", header=True, index=False)
        DF = df4
    else:
        DF = df3
        
    x_DF = DF.iloc[:, input_columns[0]:input_columns[1]]
    print("x_DF -> " + str(x_DF))
    y = DF.loc[:,target_Yname] 
    DF_columns = x_DF.columns
    print("DF_columns -> " + str(DF_columns))
    DF.to_csv("data_DF.csv", sep=",", header=True, index=False)
    print('check2 NaN -> ' + str(DF.isnull().any()))
    
    return x_DF, y, DF_columns, target_Yname

######################################## 説明変数Xを標準化
######################################## 各因子のデータを平均0,標準偏差を1へ変換
def my_scaler(X_df, my_feature_names):
    scaler = StandardScaler()
    x_sc = scaler.fit_transform(X_df)  
    # sklearnの上記のようなfit関連の実行後、データの書式がnumpy形式となりカラム(列)名が無くなる。
    # そのため、再度、pandas形式に変換する(次の2行)。これはグラフの軸名で使用するため。
    x_sc = pd.DataFrame(x_sc) # pandasデータフレーム形式へ変換
    x_sc.columns = my_feature_names # カラム名を追記
    x_sc.to_csv('data_scalerX.csv', sep=',', index=False, encoding='utf-8', header=True)
    x_sc.describe().to_csv('data_describe.csv', index=True, encoding='utf-8', mode='w', header=True)

    return x_sc

######################################## 説明変数Xを正規化
######################################## 各因子のデータ範囲を0～1へ変換
def my_mscaler(X_df, my_feature_names):
    mscaler = MinMaxScaler()
    x_msc = mscaler.fit_transform(X_df)  
    # sklearnの上記のようなfit関連の実行後、データの書式がnumpy形式となりカラム(列)名が無くなる。
    # そのため、再度、pandas形式に変換する(次の2行)。これはグラフの軸名で使用するため。
    x_msc = pd.DataFrame(x_msc) # pandasデータフレーム形式へ変換
    x_msc.columns = my_feature_names # カラム名を追記
    x_msc.to_csv('data_scalerX.csv', sep=',', index=False, encoding='utf-8', header=True)
    x_msc.describe().to_csv('data_describe.csv', index=True, encoding='utf-8', mode='w', header=True)

    return x_msc

######################################## ニューラルネットワークのモデル
def build_model(input_factors, myLoss, observation_metrics, myOptimizer):
    model = Sequential()
    model.add(Dense(units=1000, input_dim=input_factors, activation='relu'))
    model.add(Dense(units=800, activation='relu'))
    model.add(Dropout(rate=0.5))
    model.add(Dense(units=100, activation='relu'))
    model.add(Dense(1))

    model.compile(loss=myLoss, optimizer=myOptimizer, metrics=observation_metrics)
    # lossに指定した評価指標は、損失関数と呼ぶ。
    # 一方のmetricsに指定した評価指標は、観察してるだけで評価関数と訳されてる模様。
    model.summary()

    return model

######################################## 実行
def run(X_df, Y, model, name, fit_parameter, observation_metrics, i=1):
    myBatch_size, myEpochs = fit_parameter

    start = time.time()
    plot_model(model, to_file = str(i) + "_" + name + "_model.png", show_shapes=True)

    ### check poit
    checkpoint_path = os.path.join("ckpt", str(i) + "_" + name + "_cp-{epoch:03d}.ckpt")
    checkpoint_dir = os.path.dirname(checkpoint_path)
    callbacks = []
    callbacks.append(EarlyStopping(monitor='loss',
                                   min_delta=0,
                                   patience=20,
                                   verbose=0,
                                   mode='auto'))
    callbacks.append(ModelCheckpoint(filepath = checkpoint_path,
                                     monitor = 'loss',
                                     verbose = 0,
                                     save_best_only = True,
                                     save_weights_only = False,
                                     mode = 'min',
                                     period = 30))
    callbacks.append(CSVLogger(str(i) + "_" + name + "_history.csv"))

    ### fit
    history = model.fit(X_df, Y,
                        batch_size = myBatch_size,
                        epochs = myEpochs,
                        verbose=0,
                        callbacks = callbacks)
    print(history.history)
    y_list = history.history["loss"]
    ax1.plot(y_list, label=str(i) + ". " + name) #, color="r")
    ax1.legend(loc = "upper right")
    #plt.show
    fig.savefig('00_history.png')

    ### evaluate
    out_loss_and_metrics = model.evaluate(X_df, Y, batch_size = myBatch_size, verbose=1)
    out_loss = out_loss_and_metrics[0]
    out_metrics = out_loss_and_metrics[1:]

    out_time = time.time() - start
    print('optimizer:', name)
    print('loss:{0:.3f}'.format(out_loss))
    print('metrics:', out_metrics)
    print('time:{0:.3f}sec'.format(out_time))
    out_metrics_str = [str(s) for s in out_metrics]
    out_metrics_str_with_csv = ','.join(out_metrics_str)
    buf = [name, str(out_loss), str(out_metrics_str_with_csv), str(out_time)]
    out_data = ','.join(buf)
    with open("00_out_data.csv", 'a') as f:
        if i == 1:
            observation_metrics_with_csv = ','.join(observation_metrics)
            f.write("name,loss," + observation_metrics_with_csv + ",time" + '\n')
        f.write(out_data + '\n')

    ### predict
    pred_Y_np = model.predict(X_df) #.flatten()
    pred_Y_np_1d = np.ravel(pred_Y_np)
    pred_Y = pd.Series(pred_Y_np_1d)

    ### plot graph
    plot_scatter(name, X_df, Y, target_Yname, X_for_plot_graph, pred_Y, i)
    plot_reg(name, Y, target_Yname, X_for_plot_graph, pred_Y, i)

    ### save model
    model.save(str(i) + "_" + name + "_model.h5")

# 相関性の確認1 オリジナル生データを青点で、回帰式によるデータを赤点でプロット
def plot_scatter(name, X_df, Y, target_Yname, X_for_plot_graph, pred_Y, cnt=1):
    #fig = plt.figure(num=2, figsize=(8, 6), dpi=100)
    fig = plt.figure(dpi=120)
    ax = fig.add_subplot(1,1,1)
    ax.set_title(str(cnt) + ' ' + name, fontsize=14)
    ax.set_xlabel(X_for_plot_graph, fontsize=16)
    ax.set_ylabel(target_Yname + '(target)', fontsize=16)
    ax.scatter(X_df.loc[:,X_for_plot_graph], Y, facecolors='none', edgecolors='b', label='Raw data')
    ax.scatter(X_df.loc[:,X_for_plot_graph], pred_Y, facecolors='none', edgecolors='r', label='RandomForest')
    ax.legend(loc='lower right', fontsize=14)
    ax.tick_params(labelsize=16)
    plt.tight_layout()
    #plt.show()
    plt.savefig(str(cnt) + "_" + name + '_scatter.png')
    plt.close() # Close a figure window

# 相関性の確認2 横軸Y：縦軸Pred_Y
def plot_reg(name, Y, target_Yname, X_for_plot_graph, pred_Y, cnt=1):
    save_fname = target_Yname + "_{:02g}_".format(cnt) + name + '_real_vs_pred.png'

    # 欠損値を省いたことにより、連番でなくなっている可能性のあるインデックスをリセット
    Y = Y.reset_index(drop=True) # drop=Trueにより、元を削除して、新規インデックスを0から連番で作成
    DF = pd.concat([Y.rename(target_Yname), pred_Y.rename(target_Yname + "_pred")], axis=1)
    DF.to_csv('data_PredY_vs_Y_' + "{:02g}".format(cnt) + ".csv", sep=',', index=False, header=True, encoding='utf-8')
    
    popt, pcov = optimize.curve_fit(approximation_expression, Y, pred_Y)
    
    ax = plt.figure(num=0, dpi=120).gca() 
    #plt.rcParams["axes.labelsize"] = 15
    ax.set_title("{:02g}_".format(cnt) + name, fontsize=14)
    ax.set_xlabel(target_Yname, fontsize=16)
    ax.set_ylabel("Pred\n" + target_Yname, fontsize=16)
    rp = ax.scatter(x = target_Yname,
                    y = target_Yname + "_pred",
                    data = DF,
                    facecolors="none",
                    edgecolors='black') # purple

    # 生データの最小値と最大値
    y_min = DF[target_Yname].min()
    y_max = DF[target_Yname].max()
    # 予測データの最小値と最大値
    y_pred_min = DF[target_Yname + '_pred'].min()
    y_pred_max = DF[target_Yname + '_pred'].max()
    # プロットするためのxデータの範囲を決める
    x_min = min(y_min, y_pred_min)
    x_max = max(y_max, y_pred_max)
    print(x_min)
    print(x_max)
    # グラフのレンジを決める
    x_range = x_max - x_min # データの範囲が100以下か以上かで分ける
    print("x_range = x_max - x_min = " + str(x_range))
    if x_max > 1:
        min_lim = 0
        if x_range <= 10:
            max_lim = math.floor(x_max + 1)
        else:
            max_lim = math.floor(x_max + 10)
        rp.axes.set_xlim(min_lim, max_lim)
        rp.axes.set_ylim(min_lim, max_lim)
    else:
        max_lim = 0
        if x_range <= 10:
            max_lim = math.floor(x_max - 1)
        else:
            max_lim = math.floor(x_max - 10)
        rp.axes.set_xlim(min_lim, max_lim)
        rp.axes.set_ylim(min_lim, max_lim)

    # 近似式プロットのためのデータを作成
    x_approximation = np.linspace(min_lim, max_lim, 10) # numpyでxデータを作成
    y_approximation = popt[0] * x_approximation # 近似式に代入してyデータを作成
    line_approximation = ax.plot(x_approximation, y_approximation,
                                 linestyle = 'dashed', linewidth = 3, color='r') 
    
    rp.axes.set_aspect('equal', adjustable='box')
    plt.grid(True)
    # 凡例を指定して記入する場合は、リストで指定する
    ax.legend([line_approximation[0]], ["y = {:.2f}x".format(popt[0])],
               loc='upper left',
               numpoints=1,
               fontsize=15)
    plt.tick_params(labelsize=15)
    plt.tight_layout() 
    plt.savefig(save_fname) 
    plt.close()

### 近似式 y = ax 用の関数 
def approximation_expression(x, a):
    return a * x

def create_instance_history_fig(loss, epochs):
    global fig, ax1

    fig = plt.figure()
    fig.suptitle("Optimizer comparison", fontsize = 14)
    ax1 = fig.add_subplot(1,1,1)
    ax1.set_xlabel("epoch", size = 14)
    ax1.set_ylabel(loss, size = 14)
    ax1.set_xlim(0, epochs)
    ax1.set_yscale('log')
    ax1.grid(True)
    legend = ax1.legend()
    legend.set_title('Optimizer')
    fig.set_dpi(120)

######################################## 指定ファイルを複数一括削除
def my_del_files():
    for f1 in glob.glob('*.png'):
        os.remove(f1)
    for f2 in glob.glob('data_*'):
        os.remove(f2)

######################################## main関数
def main():
    ### 初期化(前回実行時に生成されたファイルを削除)
    my_del_files()
    
    ### データセットからデータ抽出
    data = my_load_csv(path, check_column, check_duplicates_columns,
                       start_row, input_columns, out_start_column,
                       target_Yname, target_value2ln)
    
    X = data[0]
    Y = data[1]
    feature_names = data[2]
    print('### extraction data')
    print(X)
    print(Y)
    print(feature_names)

    ### 説明変数Xを標準化,正規化する？
    if scalerjudge == 1: # 標準化する場合
        X = my_scaler(X, feature_names)
    elif scalerjudge == 2: # 正規化する場合
        X = my_mscaler(X, feature_names)
    else:
        pass
    print('### scalerX data')
    print(X)

    create_instance_history_fig(myLoss, fit_parameter[1])
    i = 1
    for name, myOptimizer in all_optimizer.items():
        print(name)
        print(myOptimizer)
        model = build_model(input_factors, myLoss, observation_metrics, myOptimizer)
        run(X, Y, model, name, fit_parameter, observation_metrics, i)
        i += 1

if __name__ == "__main__":
    ### parameter
    path = "boston_XYdata.csv" # 読み込む生データcsvファイル    
    input_factors = 13
    target_Yname = 'PRICE' # 目的変数Y
    X_for_plot_graph = 'RM'

    ### 生データcsvファイルの読み込み設定
    start_row = 1 # 読み込み開始行。1行目なら1, 2行目なら2, …
    target_value2ln = False # Yを対数変換するか？
    check_flag = True # True:NaNやNoneなどの欠損値を除外する場合。False:指定文字列を残す場合
    check_column = (check_flag, target_Yname, 'ok') # 要素3は、Falseの場合に機能する。指定文字列を残す
    check_duplicates_columns=[] # ダブりを省く列名を指定。複数可。省かない場合は[]と空リストにする
    
    input_columns = (0, input_factors) # 設計因子の列番号(開始, 終了)
    out_start_column = 14 # この列以降のデータについて、欠損値の処理を行う
    scalerjudge = 0 # データXを標準化する場合は1, 正規化する場合は2, しない場合はその他任意

    ### データ分析の設定
    myLoss = 'mean_squared_error'
    observation_metrics = ['mean_squared_error', 'mean_absolute_error',
                  'mean_absolute_percentage_error', 'mean_squared_logarithmic_error',
                  'squared_hinge', 'hinge', 'categorical_hinge', 'logcosh']
    fit_parameter = (30, 250) #(batch_size, epochs)

    ### オプティマイザの設定
    all_optimizer = {
        #'SGD': tf.keras.optimizers.SGD(lr=0.01, momentum=0.0, decay=0.0, nesterov=False),
        'RMSprop': tf.keras.optimizers.RMSprop(lr=0.001, rho=0.9, epsilon=None, decay=0.0),
        'Adagrad': tf.keras.optimizers.Adagrad(lr=0.01, epsilon=None, decay=0.0),
        'Adadelta': tf.keras.optimizers.Adadelta(lr=1.0, rho=0.95, epsilon=None, decay=0.0),
        'Adam': tf.keras.optimizers.Adam(lr=0.001, beta_1=0.9, beta_2=0.999, epsilon=None, decay=0.0, amsgrad=False),
        'Adamax': tf.keras.optimizers.Adamax(lr=0.002, beta_1=0.9, beta_2=0.999, epsilon=None, decay=0.0),
        'Nadam': tf.keras.optimizers.Nadam(lr=0.002, beta_1=0.9, beta_2=0.999, epsilon=None, schedule_decay=0.004),
    }

    ### call function
    main()

    print('finished')

●参考資料

下記はkerasの本家で簡潔に書かれていてわかりやすい。

keras.io

下記二つはテンソルフロー2.0の簡潔な説明である。

初心者のための TensorFlow 2.0 入門  |  TensorFlow Core

エキスパートのための TensorFlow 2.0 入門  |  TensorFlow Core

そして、下記はテンソルフローのトップページ 

www.tensorflow.org

指標は下記です。

keras.io

●関連リンク

hk29.hatenablog.jp

以上

<広告>

リンク