線型基底関数モデル① 多項式近似

線型基底関数モデル① 多項式近似

線型基底関数モデル

新しい年齢と体重のデータ

 前回までは 25歳以下に限定した年齢と体重のデータを使っていました。今回は年齢制限をなくして、幅広い年齢層から 20人を無作為抽出したデータを使用します。

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

# 線形基底関数モデル リストM3-A-1

# モジュールをインポート
import numpy as np
import matplotlib.pyplot as plt

# 20人の年齢と体重のデータセット
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])

 
 まずは与えられたデータをプロットしてみます。

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

# 線形基底関数モデル リストM3-A-2
# 動作条件:リストM3-A-1の実装

# Figureを作成
fig = plt.figure(figsize = (8, 6))

# FigureにAxesを1つ追加
ax = fig.add_subplot(111)

# Axesのタイトルを設定
ax.set_title("Age vs. Weight", 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)

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

plt.show()

Python 年齢と体重のデータプロット

 このグラフを直線で近似するのは明らかに無理があります。人間の身体的成長には限界があるので、ある年齢で体重の増加が頭打ちになることは常識的観点からもわかります。そこで、上のグラフを近似できるような別の関数を探すことにします。
 

基底関数モデル

 とはいえ、「どのような関数をもってくればデータにフィットさせることができるのか」を判断することは難しいので、いくつかの関数を組合わせて近似曲線をつくることを考えます。すなわち、基底関数 とよばれる $M$ 個の関数のセット
 
\[\phi_0(x),\quad \phi_1(x),\quad ... \quad \phi_{M-1}(x)\]
を用意して、各関数に 重み をつけて足し合わせ、
 
\[f(x)=a_0\phi_0(x)+a_1\phi_1(x)+a_2\phi_2(x)+\cdots +a_{M-1}\phi_{M-1}(x)+a_M\]
という関数を近似関数とします。基底関数の選択に関しては、扱う問題ごとにプログラマーに委ねられますが、判断しにくい場合は多項式やガウス関数、三角関数などの汎用性の高い関数を与えておきます。係数 $a_k$ は基底関数の重みとよばれるパラメータで、こちらは最小二乗法のアルゴリズムにしたがって計算機のほうで調整してくれます。うまくフィットしない場合は、基底関数の種類や数を変えて再計算させるなど、ある程度の試行錯誤は必要です。

 パラメータ $a_k$ の求め方は前回と同じです。
 次のような $N$ 個の入力データ $\boldsymbol{x}$ と目標データ $\boldsymbol{y}$ が与えられたとします。
 
\[\boldsymbol{x}=\begin{bmatrix}
x_0\\x_1\\\vdots\\x_{N-1}\end{bmatrix}
,\quad\boldsymbol{y}=\begin{bmatrix}
y_0\\y_1\\\vdots\\y_{N-1}\end{bmatrix}\]
 ここでは記述を簡単にするために、$2$ つの基底関数を用いて
 
\[f(x)=a_0\phi_0(x)+a_1\phi_1(x)+a_2\]
という関数を用意します。常に $1$ をとるような $3$ つめの基底関数 $\phi_2(x)$ を定義すると、
 
\[f(x)=a_0\phi_0(x)+a_1\phi_1(x)+a_2\phi_2(x)\]
と表すことができます。入力データベクトルに対しては、
 
\[\boldsymbol{f}(\boldsymbol{x})=
a_0\begin{bmatrix}
\phi_0(x_0)\\\phi_0(x_1)\\\vdots\\\phi_0(x_{N-1})
\end{bmatrix}
+a_0\begin{bmatrix}
\phi_1(x_0)\\\phi_1(x_1)\\\vdots\\\phi_1(x_{N-1})
\end{bmatrix}
+a_2\begin{bmatrix}
\phi_2(x_0)\\\phi_2(x_1)\\\vdots\\\phi_2(x_{N-1})
\end{bmatrix}\]
となります。ここで、
 
\[\boldsymbol{\phi_0}=\begin{bmatrix}
\phi_0(x_0)\\\phi_0(x_1)\\\vdots\\\phi_0(x_{N-1})\end{bmatrix}
,\quad\boldsymbol{\phi_1}=\begin{bmatrix}
\phi_1(x_0)\\\phi_1(x_1)\\\vdots\\\phi_1(x_{N-1})\end{bmatrix}
,\quad\boldsymbol{\phi_2}=\begin{bmatrix}
\phi_2(x_0)\\\phi_2(x_1)\\\vdots\\\phi_2(x_{N-1})\end{bmatrix}
=\begin{bmatrix}1\\1\\\vdots\\1\end{bmatrix}\]
というベクトルを定義すると、
 
\[\boldsymbol{f}(\boldsymbol{x})=a_0\boldsymbol{\phi_0}+a_1\boldsymbol{\phi_1}+a_2\boldsymbol{\phi_2}\]
と表せます。$\boldsymbol{f}(\boldsymbol{x})$ と実際の目標データ $\boldsymbol{y}$ の差の平方、
 
\[|\boldsymbol{y}-\boldsymbol{f}(\boldsymbol{x})|^2=|\boldsymbol{y}-(a_0\boldsymbol{\phi_0}+a_1\boldsymbol{\phi_1}+a_2\boldsymbol{\phi_2})|^2\]
を最小にするようにパラメータ $a_k$ を決めます。$a_0\boldsymbol{\phi_0}+a_1\boldsymbol{\phi_1}+a_2\boldsymbol{\phi_2}$ は $N$次元空間内で
 
\[\boldsymbol{\phi_0},\quad \boldsymbol{\phi_1},\quad \boldsymbol{\phi_2}\]
の 3本のベクトルが張る空間内 $V$ を動くベクトルです。この点を $\mathrm{P}$ とおきます。$\boldsymbol{y}$ は $N$次元空間内の固定された点です。この点を $\mathrm{Q}$ とおくと、
 
\[\vec{\mathrm{PQ}}=\boldsymbol{y}-(a_0\boldsymbol{\phi_0}+a_1\boldsymbol{\phi_1}+a_2\boldsymbol{\phi_2})\]
となります。$|\vec{\mathrm{PQ}}|$ が最小となるような点 $\mathrm{P}$ は、$\mathrm{Q}$ から 空間 $V$ に下ろした垂線の足 $\mathrm{H}$ と一致します。空間に下ろす垂線をイメージするのは難しいですが、要するに $\vec{\mathrm{PQ}}$ は $\boldsymbol{\phi_0},\ \boldsymbol{\phi_1},\ \boldsymbol{\phi_2}$ の 3本のベクトルと内積をとると、いずれも $0$ になるということです。
 
\[\begin{align*}(\boldsymbol{y}-a_0\boldsymbol{\phi_0}-a_1\boldsymbol{\phi_1}-a_0\boldsymbol{\phi_2})\cdot\boldsymbol{\phi_0}=0\\[6pt]
(\boldsymbol{y}-a_0\boldsymbol{\phi_0}-a_1\boldsymbol{\phi_1}-a_0\boldsymbol{\phi_2})\cdot\boldsymbol{\phi_1}=0\\[6pt]
(\boldsymbol{y}-a_0\boldsymbol{\phi_0}-a_1\boldsymbol{\phi_1}-a_0\boldsymbol{\phi_2})\cdot\boldsymbol{\phi_2}=0
\end{align*}\]
 これを行列形式で書き直すと次のようになります。
 
\[\begin{bmatrix}\boldsymbol{\phi_0}\cdot\boldsymbol{\phi_0}
& \boldsymbol{\phi_1}\cdot\boldsymbol{\phi_0}
& \boldsymbol{\phi_2}\cdot\boldsymbol{\phi_0}\\
\boldsymbol{\phi_0}\cdot\boldsymbol{\phi_1}
& \boldsymbol{\phi_1}\cdot\boldsymbol{\phi_1}
& \boldsymbol{\phi_2}\cdot\boldsymbol{\phi_1}\\
\boldsymbol{\phi_0}\cdot\boldsymbol{\phi_2}
& \boldsymbol{\phi_1}\cdot\boldsymbol{\phi_2}
& \boldsymbol{\phi_2}\cdot\boldsymbol{\phi_2}\end{bmatrix}
\begin{bmatrix}a_0\\a_1\\a_2\end{bmatrix}
=\begin{bmatrix}\boldsymbol{y}\cdot\boldsymbol{\phi_0}
\\\boldsymbol{y}\cdot\boldsymbol{\phi_1}
\\\boldsymbol{y}\cdot\boldsymbol{\phi_2}\end{bmatrix}\tag{1}\]
 記述を簡単にするために $N=3$ とします(一般の $N$ でも以下の変形手順は同じです)。
 基底関数行列 $\Phi$ とパラメータベクトル $\boldsymbol{p}$ を
 
\[\Phi=\begin{bmatrix}\boldsymbol{\phi_0}(x_0)&\boldsymbol{\phi_1}(x_0)&\boldsymbol{\phi_2}(x_0)\\
\boldsymbol{\phi_0}(x_1)&\boldsymbol{\phi_1}(x_1)&\boldsymbol{\phi_2}(x_1)\\
\boldsymbol{\phi_0}(x_2)&\boldsymbol{\phi_1}(x_2)&\boldsymbol{\phi_2}(x_2)\end{bmatrix}
,\quad \boldsymbol{p}=\begin{bmatrix}a_0\\a_1\\a_2\end{bmatrix}\]
のように定義すると、式 (1) は
 
\[\Phi^T\Phi\boldsymbol{p}=\Phi^T\boldsymbol{y}\]
と表すことができます。両辺に左側から $\Phi^T\Phi$ の逆行列をかけると
 
\[\boldsymbol{p}=(\Phi^T\Phi)^{-1}\Phi^T\boldsymbol{y}\tag{2}\]
となります。$(\Phi^T\Phi)^{-1}\Phi^T$ は前回記事でも現れた疑似逆行列(一般化逆行列)です。
 

多項式近似

 曲線を $x$ の $n$ 次式で近似することを多項式近似とよびます。前回までは $x$ の $1$ 次式で近似していたので、線形回帰も多項式近似に含まれます。今回は曲線を $2$ 次式で近似することを試みます。すなわち基底関数として
 
\[\phi_0(x)=x^2,\quad \phi_1(x)=x\]
を選んで、データを
 
\[ax^2+bx+c\]
という2次関数で近似するようなパラメータ $a,\ b,\ c$ を決定します。

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

# 線形基底関数モデル リストM3-A-3
# 動作条件:リストM3-A-1の実装

# Φ0とΦ1の設定
def phi_0(x):
    return x ** 2

def phi_1(x):
    return x

# Φの転置行列
# 3行目の要素はすべて1
m1 = np.array([phi_0(x), phi_1(x), np.ones(len(x))])

# m1と[m1の転置行列]の積
m2 = m1.dot(m1.T)

# m2の逆行列
m3 = np.linalg.inv(m2)

# m3とm1の積
m4 = m3.dot(m1)

# パラメータベクトルの決定
p = m4.dot(y.T)

# 近似曲線の係数を表示
print("[a, b, c] = {}".format(np.round(p, 4)))

# 平均2乗誤差mseを計算
fx = p[0]*x**2 + p[1]*x + p[2]
mse = np.mean((fx-y)**2)

# 標準偏差SDを表示
print("SD = {0:.2f} kg".format(np.sqrt(mse)))

# Figureを作成
fig = plt.figure(figsize = (8, 6))

# FigureにAxesを1つ追加
ax = fig.add_subplot(111)

# Axesのタイトルを設定
ax.set_title("Age vs. Weight", 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)

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

# 近似曲線の描画
x2 = np.linspace(0, 70, 100)
y2 = p[0]*x2**2 + p[1]*x2 + p[2]
ax.plot(x2, y2, color = "red")

plt.show()
[a, b, c] = [-0.0343  3.2275 10.0832]
SD = 5.68 kg

線形基底関数モデル 年齢と体重の多項式近似

 標準偏差 SD はフィッティング精度を示す値で、この値が大きいほど実際のデータと近似曲線のずれが大きいことを意味します。多項式近似からは少し離れてしまいますが、3 つめの基底関数として対数関数 $\log x$ を加えてみると、もう少しだけ精度が上がります。

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

# 線形基底関数モデル リストM3-A-4
# 動作条件:リストM3-A-1の実装

# Φ0とΦ1の設定
def phi_0(x):
    return x ** 2

def phi_1(x):
    return x

def phi_2(x):
    return np.log(x)

# Φの転置行列
# 3行目の要素はすべて1
m1 = np.array([phi_0(x), phi_1(x), phi_2(x), np.ones(len(x))])

# m1と[m1の転置行列]の積
m2 = m1.dot(m1.T)

# m2の逆行列
m3 = np.linalg.inv(m2)

# m3とm1の積
m4 = m3.dot(m1)

# パラメータベクトルを決定
p = m4.dot(y.T)

# 係数を表示
print("[a, b, c, d] = {}".format(np.round(p, 4)))

# 平均2乗誤差mseを計算
fx = p[0] * phi_0(x) + p[1] * phi_1(x) + p[2] * phi_2(x) + p[3]
mse = np.mean((fx-y)**2)

# 標準偏差SDを表示
print("SD = {0:.2f} kg".format(np.sqrt(mse)))

# Figureを作成
fig = plt.figure(figsize = (8, 6))

# FigureにAxesを1つ追加
ax = fig.add_subplot(111)

# Axesのタイトルを設定
ax.set_title("Age vs. Weight", 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)

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

# 近似曲線を描画
x2 = np.linspace(1, 70, 100)
y2 = p[0] * phi_0(x2) + p[1] * phi_1(x2) + p[2] * phi_2(x2) + p[3]
ax.plot(x2, y2, color = "red")

plt.show()
[a, b, c, d] = [-0.0243  2.0948  9.5511  2.0261]
SD = 4.90 kg

Python機械学習 多項式近似

 基底関数の組合わせ方は無数にあるので、最適な近似関数を見つけることは事実上不可能ですが、今後の記事でオーバーフィッテングの問題と絡めて、より良いモデルを選択する方法(ベストではなくベターなモデルに改良していく方法)について解説していく予定です。