【tensorflowjs x kerasの使い方】正規分布の評価〜非線形回帰解析②

前回

線型回帰の手順を解説しました。それでは今回はいよいよ待望(?)の解析的な手法による非線形回帰の実装のやり方を説明していこうと思います。

【tensorflowjs & kerasの使い方】燃費を予測するサンプルを理解する〜非線形回帰解析①

tensorflowjsとkerasで、データを非線形回帰解析する手順を考えていきます。

【効果的学習法レポート】シェルスクリプトをこれから学びたい人のためのオススメ書籍＆教材特集

ガウス分布ノイズの実験データの生成

今回、tensorflowとkerasで、ガウス分布ノイズの度数分布の正規化済み生データから以下の解析式で非線形回帰を試みます。
既に関連記事として、

シェルコマンド x 機械学習 | awkで機械学習で使える高速データ処理〜ガウス分布ノイズの生成方法

の内容の中でも解説しておりますので、実験データセットとして今回利用します。
どうやってデータセットを生成しているか興味がお有りの方はそちらの内容をご参照ください。

【シェルスクリプトで機械学習】Awkで機械学習で使える高速データ処理〜ガウス分布ノイズの生成方法

ガウス分布ノイズのデータセット生成の方法を例にとり、Awkで学習用の訓練データを膨大かつ高速で下処理するテクニックを考えます。

それでは正規分布$$N(0,1)$$に従う正規乱数を100万回のサンプリングで度数分布データを得るawkスクリプトを以下に再掲しておきます。

開く/閉じる

            $ awk -v seed="${RANDOM}" -v iter=1000000 -v inp_mu=0.0 -v inp_sigma=1.0 '
    function box_muller(mu, sigma, arr_len_) {
        PI = 3.14159265359;
        count = 0;
        while(count < arr_len_) {
            x1 = rand();
            x2 = rand();
            tmp_rand_val = sqrt(-2 * log(x1)) * cos(2 * PI * x2);
            gauss_arr[count] = mu + sigma * tmp_rand_val;
            count++;
        }
    }
    BEGIN{
        srand(seed); #シード値をリセット
        box_muller(inp_mu, inp_sigma, iter);
        arr_len = length(gauss_arr);
        for (i = 0; i < arr_len; i++) {
            print gauss_arr[i];
        }
    }
' | awk -v binsize=0.2 '
    BEGIN {
        OFS=","
    }
    {
        if(binsize <= 0) exit;
        if($1 < 0) {
            frequency[int($1 / binsize) - 1] ++;
        } else {
            frequency[int($1 / binsize)] ++;
        }
    }
    END {
        for(i in frequency) {
            print (i + 0.5) * binsize, frequency[i] | "sort -n";
        }
    }
' | awk -F "," '
    BEGIN {
        OFS=",";
        count = 0;
    }
    {
        xarr[count] = $1;
        yarr[count] = $2;
        if(ymax<$2) ymax=$2;
        count++;
    }
    END {
        xarr_len = length(xarr);
        for (i = 0; i < xarr_len; i++) {
            print i, xarr[i], yarr[i] / ymax
        }
    }
' > norm_dist_exp.csv

データをプロットすると以下のような図になります。

なお既にこのデータセットは正規化していますので、本来の正規分布の振幅部$$A = {1 \over \sqrt{2\pi\sigma^2}}$$の情報は失わて$$A = 1$$としてしまったため、正規化されたガウス分布の確率密度関数で以下のような解析式を使います。

\displaystyle{ f(x) = \exp \biggl( -\frac{(x - \mu)^2}{2 \sigma^2} \biggr) }

Eq. (1)

我々としてはこのダミーデータセットは平均値$$\mu = 0$$と標準偏差$$\sigma = 1$$が正解であることをソースコードの中身から既に知っていますが、解析してくれるtensorflow/kerasモデルはそれを知りません。
非線形回帰解析によって見事正解を導き出すことができるのかを以降で検証してみます。

機械学習の基本と関連する技術をバランス良く学ぶ図解即戦力機械学習&ディープラーニングのしくみと技術がこれ1冊でしっかりわかる教科書

【効果的学習法レポート】シェルスクリプトをこれから学びたい人のためのオススメ書籍＆教材特集

KerasモデルのLayerをカスタマイズする

通常、tensorflowjsのSquentialモデルは、Layarクラスのインスタンスを一つのブロックに見立てて、複数のブロックを入力層から順に繋げていって出力層までちゃんと繋がるように組立しなくてはいけません。
デフォルトで利用できる

Layersクラス

とも呼ばれて、現在は以下の種類が対応しています。

開く/閉じる

            linear(デフォルト)
elu
hardSigmoid
relu
relu6
selu
sigmoid
softmax
softplus
softsign
tanh

前回の線型回帰では活性化関数を

linear

にしていました。
これはバイアス項ありの場合、内部に$$y = ax + b$$という解析式と等価になりますので、入力値$$x$$に対し一次関数として変換された出力値$$y$$が得られる、というお話でした。
今回は正規分布関数を活性化関数として扱いたいのですが、でもこれってtensorflowjsの活性化関数のリストに無いやつじゃん、と思われるかもしれません。
...そうです。
残念ながら無いタイプの活性化関数は自分でカスタマイズは避けて通れません。
ということで少しだけプログラミングする手間がかかりますが、以降の内容から低レベルなAPIである

tf.Layersクラス

を拡張をやっていきます。

カスタムLayerの作り方

まずはおらが専用モデルのソースコードファイルを追加して、以下の内容で編集してみます。
とりあえずファイル名は

normDistLayer.js

にしてみます。

開く/閉じる

            import * as tf from '@tensorflow/tfjs';

export class NormDistLayer extends tf.layers.Layer {
    constructor(config) {
        super(config);
        this.mu0 = config.mu0;
        this.sigma0 = config.sigma0;
    }

    /**
    * build関数では最適化パラメータの初期化処理等を行う
    */
    build(inputShape) {
        // 最適化パラメーターはスカラなのでconst shape = [];としてもOK
        const shape = [inputShape[inputShape.length - 1]];
        this.mu = this.addWeight('mu', shape, 'float32', tf.initializers.constant({value: this.mu0}));
        this.sigma = this.addWeight('sigma', shape, 'float32', tf.initializers.constant({value: this.sigma0}));
    }

    /**
    * call関数で活性化関数を内部で定義する。
    * 基本はTensor型の入力値を入れて、Tensor型の出力を出す。
    * メモリーリークを避けるためにtidyは必ず使うこと
    */
    call(input) {
        return tf.tidy(() => {
            const p1 = tf.pow(input[0].sub(this.mu.read()), 2);
            return tf.div(p1, this.sigma.read().pow(2).mul(-2)).exp();
        });
    }

    /**
    * getConfig関数はLayerインスタンスの読み込みと書き出しの時に
    * 必須となるプロパティをJSON形式でシリアライズ・デシリアライズする
    */
    getConfig() {
        const config = super.getConfig();
        Object.assign(config, {
            mu0: this.mu0,
            sigma0: this.sigma0
        });
        return config;
    }

    /**
    * Layersのライブラリに登録するクラス名
    */
    static get className() {
        return 'NormDistLayer';
    }
}

基本的に

tf.layers.Layerクラス

を継承し、

build

call

getConfig

などのメンバ関数をオーバーライトしていけばカスタムレイヤーが作成できます。
このカスタムレイヤーを使う側の設定は以下のように書けます。
なお

前回の線型回帰解析で利用したコード

を

nonlinear_regression.js

として再利用しています。

開く/閉じる

            import * as tf from '@tensorflow/tfjs';
//👇先程作成したカスタムレイヤークラスをインポート
import { NormDistLayer } from './normDistLayer';

//👇カスタムクラスをtensorflowで呼び出せるように登録する(重要)
tf.serialization.registerClass(NormDistLayer);

export async function evalNonlinaer(rawData) {
    const model = tf.sequential();

    //👇非線形回帰の場合、カスタムレイヤーが入力層かつ出力層
    model.add(new NormDistLayer({
        inputShape: [1],
        mu0: 0.1,
        sigma0: 1.1
    }));

    console.log('Analysis started.');

    const originalData = [];
    const tmpLabels = [];
    for (const row in rawData) {
        tmpLabels.push(rawData[row][0]);
        originalData.push({x: parseFloat(rawData[row][1]), y: parseFloat(rawData[row][2])});
    };

    // 入力用データセット配列をTensor型に変換
    const tensorData = convertToTensor(originalData);

    // 非線形回帰解析の実行
    await trainModel(model, tensorData.inputs, tensorData.labels);

    // 解析済みモデルから結果を得る
    const predictedData = testModel(model);

    console.log('Done.');

    // Chart.js描画用にラベル・元データセット・解析済みデータセットを返す
    return { tmpLabels, originalData, predictedData };
}

/// モデル・入力（インプット）テンソル・出力（ラベル）テンソルを指定してモデルの学習を行う関数
async function trainModel(model, inputs, labels) {
    model.compile({
        // 最適化法はAdamに指定
        optimizer: tf.train.adam(),
        // 損失関数をMSEに指定
        loss: tf.losses.meanSquaredError,
        // 学習とテストに用いる指標の名前を定義
        metrics: ['mse'],
    });

    // バッチサイズ
    const batchSize = 48;
    // エポック数
    const epochs = 500;

    // エポック毎に解析中のパラメーターを確認出来るようにcallbacksを指定
    return await model.fit(inputs, labels, {
        batchSize,
        epochs,
        shuffle: true,
        callbacks: {
            onEpochEnd: async(epoch, logs) => {
                // 繰り返し回数 - 損失 - 誤差の順にを出力
                console.log(`Epoch#${epoch} : Loss(${logs.loss}) : mse(${logs.mse})`);

                // 現在平均値μと標準偏差σの値を出力
                const muval = model.layers[0].getWeights()[0].dataSync()[0];
                const sgmval = model.layers[0].getWeights()[1].dataSync()[0];
                console.log(`mu=${muval} : sgm=${sgmval}`);
            }
        }
    });
}

/// 学習済みモデルから解析後のプロット用データを生成
function testModel(model) {
    const [xs, preds] = tf.tidy(() => {
        const xs = tf.linspace(-5, 5, 100);
        const preds = model.predict(xs.reshape([100, 1]));
        return [xs.dataSync(), preds.dataSync()];
    });

    return Array.from(xs).map((val, i) => {
        return {
            x: val,
            y: preds[i]
        }
    });
}

/// 学習データをテンソルに変換する関数
function convertToTensor(data) {
    return tf.tidy(() => {
        // データのシャッフル
        tf.util.shuffle(data);

        // データを配列に格納してから展開してテンソルに変換
        const inputs = data.map(d => d.x)
        const labels = data.map(d => d.y);
        const inputTensor = tf.tensor2d(inputs, [inputs.length, 1]);
        const labelTensor = tf.tensor2d(labels, [labels.length, 1]);

        return { inputs: inputTensor, labels: labelTensor }
    });
}

後はメインコードでこのnonlinear_regression.jsを読み込んで実行すると非線形回帰解析が走ります。

機械学習の基本と関連する技術をバランス良く学ぶ図解即戦力機械学習&ディープラーニングのしくみと技術がこれ1冊でしっかりわかる教科書

【効果的学習法レポート】シェルスクリプトをこれから学びたい人のためのオススメ書籍＆教材特集

正規化済みの正規分布の解析

先程のソースコードを使ってプログラムをビルド・実行してみます。
上記で作成していた

norm_dist_exp.csv

のデータセットを用いた解析を走らせると、以下のような結果を得ました。

開く/閉じる

            ...
Epoch#499 : Loss(0.000004175022240815451) : mse(0.000004175022695562802)
mu=0.0016643998678773642 : sgm=0.9995725154876709

ほとんど理論値に近い値を精度良く推測できていると思います。

まとめ

今回はtensorflowjs/kerasでも、非線形解析のソルバーがきちんと実装できることを実例を持って示してみました。
解析関数が理論的に厳密な形で決めてしまえるのであれば、今回ご紹介したカスタムレイヤーを使わない手はないです。
tensorflowjsには非常に優れた最適化ソルバーも充実しているため、上手くモデルにハマれば通常のデープラーニングで重い演算よりも圧倒的な速度で厳密解に収束してくれるのも利点と言えます。