非線形モデル① 非線形ガウス基底モデル

非線形モデル① 非線形ガウス基底モデル

非線形モデル

 個々の基底関数が入力データ $x$ 以外の内部パラメータ $p,\ q,\ r,\ ...$ を含んでいる場合、これを
 
\[\phi(x,\ p,\ q,\ r,\ ...)\]
と表すことにします。たとえば、ガウス基底の場合は $\mu$ と $\sigma$ が内部パラメータにあたります。線形モデルにおいては基底関数同士の線形結合
 
\[f(x)=\sum_{k}a_k\phi_k(x,\ p_k,\ q_k,\ r_k,\ ...)\]
を作り、内部パラメータを固定して $a_k$ についてのみ最適化しました。一方で $a_k$ と内部パラメータ $\ p_k,\ q_k,\ r_k,\ ...$ をすべて最適化する場合は 非線形モデル になります。

 非線形モデルは内部パラメータを自由に調整することで、線型モデルよりも柔軟な関数近似を行なうことができますが、最適化手順を簡潔な数式で表すことが難しい(線型代数では扱えない)ので、一般的には(この講座で最初に扱った)勾配降下法 のアルゴリズムを用いて最適化を行ないます(後述するように、この記事では別の方法を選択します)。

 前回まで扱っていたランダムガウス基底モデルでは、乱数を使って内部パラメータを調整していましたが、最適化はあくまで係数 $a_k$ について行なっていたので線形モデルです。乱数を使うことによって、非線形モデルの柔軟性を線形モデルに組み込んでいたのです。しかし、安定した結果を得るためには試行回数を十分に大きくとる必要があるので、大規模なデータを扱う場合には計算に時間がかかるという短所があります。

Powell's method

 今回は乱数を使わずに、$\mu$ と $\sigma$ もパラメータとして動かす非線形ガウス基底モデル
 
\[f(x)=a_0+\sum_{k=1}^{3}a_k\exp\left\{-\frac{(x-\mu_k)^2}{2\sigma_{k}^{2}}\right\}\]
を考えます。前回までのモデルとは異なり、$\mu$ と $\sigma$ の取り得る幅に制限はありません。それぞれの基底関数の中心位置も幅も自在に変えながら重ね合わせて、与えられたデータに最も適合する合計 10 個のパラメータを決定することになります。
 
 これまではライブラリに頼らずにゼロからモデルを構築してきましたが、今回は少し楽をすることにして、Scipy を活用して最適化計算を行なうことにします。

 scipy.optimize.minimize() に関数と入力データを渡すと、最適化されたパラメータを簡単に取得することができます。まずは計算に必要な関数を定義しておきましょう。

# https://python.atelierkobato.com/powell/

# 非線形モデル リストM4-A-1

import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import minimize

# ガウス関数を定義
def gauss(x, mu = 0, sigma = 1):
    return np.exp(-(x - mu)**2 / (2 * sigma**2))

# ガウス基底の線形結合
def gauss_base(x, p):
    f = p[0]
    for i in range(1, 9, 3):
        f += p[i] * gauss(x, p[i + 1], p[i + 2])
    return f

# 平均2乗誤差関数(コスト関数)
def mse_gauss_base(p, x, y):
    f = gauss_base(x, p)
    mse = np.mean((f - y)**2)
    return mse

# 最適化されたパラメータを取得する関数
# 最適化手法として powell method を指定
def fit_gauss(p_init, x, y):
    z = minimize(mse_gauss_base, p_init, args = (x, y), method = "powell")
    return z.x

 minimize()関数の第1引数には最小化したい平均2乗誤差関数を渡しています。第2引数にはパラメータベクトルの初期値、args には入力データベクトルと目標データを指定しています。method には "powell" を指定していますが、これは最適化手法として Powell's method (パウエル法) を選択したことを意味します。Powell's method は勾配降下法を用いずに関数の極値を見つけ出すことのできる手法です。
 
 リスト M4-A-1 で定義した関数に実際にデータを入力してパラメータを最適化し、近似曲線を描いてみます。

# https://python.atelierkobato.com/powell/

# 非線形モデル リストM4-A-2
# 動作条件:リストM4-A-1の実装

# 入力データベクトル
x = np.array([35, 16, 22, 43, 5,
              66, 20, 13, 52, 1,
              39, 62, 45, 33, 8,
              28, 71, 24, 18, 3])

# 目標データベクトル
y = np.array([85.19, 58.93, 64.27, 68.91, 21.27,
              68.88, 60.07, 55.18, 88.08, 8.89,
              82.31, 81.18, 80.76, 78.98, 37.55,
              75.9, 72.39, 69.51, 62.04, 12.47])

# ★★★★★★★★★★

# パラメータベクトルの初期値を設定
p_init = np.full(10, 50)

# パラメータベクトルを最適化
p = fit_gauss(p_init, x, y)

# 平均2乗誤差を計算
mse = mse_gauss_base(p, x, y)

# 標準偏差を計算
sd = np.round(np.sqrt(mse), 3)

# ★★★★★★★★★★

# FigureとAxesを作成
fig = plt.figure(figsize = (8, 6))
ax = fig.add_subplot(111)

# Axesのタイトルを設定
ax.set_title("Gaussian base (Powell method)", fontsize = 16)

# 目盛線を表示
ax.grid()

# テキストボックス書式辞書を作成
boxdic = {"facecolor" : "white", "edgecolor" : "gray",}

# テキストボックスを表示
ax.text(4, 88, "SD = {}".format(sd),
        size = 14, linespacing = 1.5, bbox = boxdic)

# 軸範囲と軸ラベルを設定
ax.set_xlim([0, 70])
ax.set_ylim([0, 100])
ax.set_xlabel("Age", fontsize = 14)
ax.set_ylabel("Weight [kg]", fontsize = 14)

# 身長と体重データの散布図
ax.scatter(x, y, color = "blue")

# 近似曲線を描画
x2 = np.linspace(0, 70, 100)
y2 = gauss_base(x2, p)
ax.plot(x2, y2, color = "red")

plt.show()

Powell's method
 

4分割交差検証

 一般に非線形モデルにおける平均2乗誤差関数(コスト関数)は複数の極値をもつことが多いので、初期値の選び方によって結果が変わる可能性もあります(このモデルの平均2乗誤差関数もデコボコしてます)が、最小値でなくても、交差検証の結果と処理時間を他のモデルと比較して、精度が良ければ採用するという方向で進めます。以下に 4分割交差検証のコードも載せておきます。

# https://python.atelierkobato.com/powell/

# 非線形モデル リストM4-A-3
# 動作条件:リストM4-A-1, リストM4-A-2の実装

# データの数
N = x.shape[0]

# 分割数を設定
K = 4

# データをK分割
xd = np.array_split(x, K)
yd = np.array_split(y, K)

# ★★★★★★★★★★

# パラメータベクトルの初期値を設定
p_init = np.full(10, 50)
p_train = np.zeros(10)

# 標準偏差の初期値を設定
sd_test = 0

for j in range(K - 1):
    x_test = xd[j]
    x_train = np.delete(x, slice(j * int(N/K), (j + 1) * int(N/K)))
    y_test = yd[j]
    y_train = np.delete(y, slice(j * int(N/K), (j + 1) * int(N/K)))
    
    # パラメータベクトルを最適化
    p_train = fit_gauss(p_init, x_train, y_train)

    # 現在選択しているテストデータの平均2乗誤差を計算
    mse_test = mse_gauss_base(p_train, x_test, y_test)

    # 現在選択しているテストデータの標準偏差を加算
    sd_test += np.sqrt(mse_test)

# 標準偏差の平均値を計算
sd_mean = np.round(sd_test/K, 3)

print("SD_test = {}".format(sd_mean))
SD_mean = 4.269

 各テストデータの標準偏差の平均値は 4.269 です。
 ランダムガウス基底モデルのときのように「パラメータを調整しながら標準偏差平均値の小さいほうを採用する」という処理を行なっていないので、4.269 は十分に評価できる精度といえます。また、こちらのモデルのほうが圧倒的に処理が速いので、大規模なデータの処理に向いているといえます。とはいえ、基底関数の選び方にはまだ工夫の余地がありそうなので、次回記事では別の非線形モデルを紹介します。

 ≫ 非線形モデル② 複数モデルの比較検討