非線形モデル

非線形モデル

非線形モデル

 個々の基底関数が入力データ $x$ 以外の内部パラメータ $p_1,\ p_2,\ p_3,\ ...$ を含んでいる場合、これを
 
\[\phi(x,\ p_1,\ p_2,\ p_3,\ ...)\tag{1}\]
と表すことにします。たとえば、ガウス基底の場合は $\mu$ と $\sigma$ が内部パラメータにあたります。線形モデルにおいては基底関数同士の線形結合
 
\[f(x)=\sum_{k}a_k\phi_k\tag{2}\]
を作り、内部パラメータを固定して $a_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 を活用して最適化計算を行なうことにします。

 最初に必要なモジュールをインポートしておきます。

# Machine_Learning_06

# In[1]

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

 ガウス関数とガウス関数の線形結合、平均二乗誤差関数を定義します。

# In[2]

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

# ガウス関数の線形結合
# pはリストやタプルで渡す
# p[0],p[1],[2],p[3],p[4],...はa,a1,μ1,σ1,...に対応する10個のパラメータ
def gaussian_basis(x, p):
    func = p[0]
    for i in range(1, 9, 3):
        func += p[i] * gauss(x, p[i+1], p[i+2])
    return func

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

 scipy.optimize.minimize() を使って、最適化されたパラメータを返す関数を定義します。

# In[3]

# 平均二乗誤差を最小にするパラメータpを返す関数
# 最適化手法としてpowell's methodを指定
def fit_gaussian(p, x, y):
    z = minimize(mse_gaussian_basis, p, args=(x, y), method="powell")
    return z.x

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

 fit_gaussian() に年齢 x と体重 y のデータを入れて、パラメータを最適化してみます。

# In[4]

# 小数点以下3桁まで表示
np.set_printoptions(precision=3, floatmode="fixed")

# 入力データベクトル
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])

# パラメータの初期値
# 10個のパラメータをすべて50に設定しておく
p = np.full(10, 50)

# 平均二乗誤差関数を最小にするようにpを最適化
p = fit_gaussian(p, x, y)

# 最適化されたパラメータベクトルpを表示
print(p)
[-1.429e+02  1.029e+02  1.779e+01  2.061e+01  1.443e+02  6.406e+01
  2.689e+01  6.826e+01  1.076e+02  8.117e+08]

 回帰曲線とデータの間の平均二乗誤差と標準偏差を求めてみましょう (曲線がどれぐらいデータにフィットしているかの目安となります)。

# In[5]

# 平均二乗誤差
mse = mse_gaussian_basis(p, x, y)

# 標準偏差
sd = np.sqrt(mse)

print("平均二乗誤差: {:.3f}".format(mse))
print("標準偏差: {:.3f}".format(sd))
平均二乗誤差: 20.446
標準偏差: 4.522

 回帰曲線をプロットしてフィッティングの様子を確認します。

# In[6]

# FigureとAxes
fig = plt.figure(figsize = (8, 6))
ax = fig.add_subplot(111)
ax.set_title("Fitting curve", fontsize = 16)
ax.grid()
ax.set_xlim([0, 70])
ax.set_ylim([0, 100])
ax.set_xlabel("Age", fontsize = 14)
ax.set_ylabel("Weight [kg]", fontsize = 14)

boxdic = {"facecolor" : "white", "edgecolor" : "gray",}
ax.text(4, 88, "SD = {:.3f}".format(sd),
        size = 14, linespacing = 1.5, bbox = boxdic)

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

# 回帰曲線をプロット
x2 = np.linspace(0, 70, 100)
y2 = gaussian_basis(x2, p)
ax.plot(x2, y2, color = "red")

plt.show()

Powell's method
 

たった1秒の最強スキル パソコン仕事が10倍速くなる80の方法

新品価格
¥1,382から
(2019/8/6 17:57時点)

4分割交差検証

 以前の講座 で定義した学習用データとテストデータに分ける関数を再掲します。

# In[7]

# 学習用データとテストデータに分ける関数
def train_test(x, y, k, r):
    idx = np.arange(0, x.shape[0])
    idx_test = idx[np.fmod(idx, k) == r]
    idx_train = idx[np.fmod(idx, k) != r]
    x_test = x[idx_test]
    x_train = x[idx_train]
    y_test = y[idx_test]
    y_train = y[idx_train]
    return x_train, x_test, y_train, y_test

 以下は 4分割交差検証のコードです。

# In[8]

# 分割数
k = 4

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

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

for r in range(k):
    x_train, x_test, y_train, y_test = train_test(x, y, k, r)

    # パラメータを最適化
    p_train = fit_gaussian(p_init, x_train, y_train)

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

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

# 標準偏差の平均値
sd_mean = sd_test / k

print("標準偏差の平均値: {:.3f}".format(sd_mean))
標準偏差の平均値: 6.070