線型回帰
次のように $N$ 人の年齢 $x_k$ と体重 $y_k$ のデータセットが用意されたとします。
\[\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}\]
ベクトル $\boldsymbol{x}$, $\boldsymbol{y}$ をそれぞれ入力データ、目標データとよびます。
ベクトルの各成分 $x_k$, $y_k$ をそれぞれ 入力変数、目標変数 とよびます。
データに存在しない年齢 $x$ が与えられたときに、体重 $y$ を予測できるような関係式を見つけることが 線型回帰 の目的です。
平均二乗誤差の定義
20人の年齢と体重のサンプルデータを用意しました。
ただし、25歳以下という年齢制限を設けてあります。
Matplotlib を使って、横軸に年齢、縦軸に体重をプロットしてみます。
# Mean_Squared_Error
# In[1]
# モジュールのインポート
import numpy as np
import matplotlib.pyplot as plt
# 20人の年齢と体重のデータセット
# 年齢は25歳以下に限定
x = np.array([9, 21, 6, 18, 5,
6, 7, 12, 10, 20,
17, 3, 24, 9, 18,
23, 12, 13, 16, 13])
y = np.array([29.01, 78.1, 18.78, 57.94, 14.81,
22.34, 26.31, 41.69, 35.25, 60.92,
69.08, 16, 66.39, 27.06, 68.03,
61.45, 42.07, 40.64, 59.62, 32.22])
# Figureを作成
fig = plt.figure(figsize = (8, 6))
# FigureにAxesを１つ追加
ax = fig.add_subplot(111)
# Axesのタイトルを設定
ax.set_title("Age vs. Weight", fontsize = 16)
# 目盛線を表示
ax.grid()
# 軸範囲を設定
ax.set_xlim([0, 25])
ax.set_ylim([0, 90])
# 軸ラベルを設定
ax.set_xlabel("Age", fontsize = 14)
ax.set_ylabel("Weight [kg]", fontsize = 14)
# 身長と体重データの散布図
ax.scatter(x, y, color = "blue")
plt.show() 
痩せすぎの人がいる一方で肥満体質の人もいるのでバラつきはありますが、全体として年齢が大きいほど体重が重くなるという傾向が見てとれます。そこで「年齢と体重は比例関係にある」と考えて、「このデータを直線で近似しよう」というのが線型回帰モデルの考え方です。上のグラフを見れば点が直線上に乗らないことは一目瞭然なので、かなり強引な手法に思えますが、ともかくも年齢 $x$ と体重 $w(x)$ の間には
\[w(x)=ax+b\]
の関係があると仮定して、実際の年齢データ $x_k$ をこの仮想的な関係式の変数 $x$ に代入します。
\[w(x_k)=ax_k+b\]
$w(x_k)$ と実際の体重データ $y_k$ の差は
\[w(x_k)-y_k=ax_k+b-y_k\]
となります。この差分の2乗を人数分足し合わせて $N$ で割った値
\[Q(a,b)=\frac{1}{N}\sum_{k=0}^{N-1}\{w(x_k)-y_k\}^2=\frac{1}{N}\sum_{k=0}^{N-1}(ax_k+b-y_k)^2\]
を定義します。$Q$ は 平均二乗誤差 (MSE;Mean Square Error) とよばれる量で、この $Q$ を最小とするようなパラメータ $a,\ b$ を 最急降下法 によって見つけます。
平均二乗誤差の可視化
最急降下法 は極小値を探すアルゴリズムなので、複数の極小値をもつような関数については、それなりに工夫が必要となります。そこでまず、$Q$ がどのような形の曲面であるかを Matplotlib で事前に調べておきます。
# In[2]
# モジュールのインポート
import itertools
from mpl_toolkits.mplot3d import axes3d
# 平均2乗誤差MSE
def mse(x, y, a, b):
return np.mean((a * x + b - y)**2)
# Figureを追加
fig = plt.figure(figsize = (10, 6))
# 3DAxeSを追加
ax = fig.add_subplot(111, projection="3d")
# Axesのタイトルを設定
ax.set_title("Mean Square Error", fontsize = 16,
position = (0.5, 1.1))
# 軸ラベルを設定
ax.set_xlabel("a", size = 14, color = "blue")
ax.set_ylabel("b", size = 14, color = "blue")
ax.set_zlabel("Q", size = 14, color = "blue")
# 解像度(resolution)を設定
res = 100
# 格子点の作成
a = np.linspace(-25, 25, res)
b = np.linspace(0, 100, res)
A, B = np.meshgrid(a, b)
# すべての要素が0のres×resサイズの配列
Q = np.zeros((len(a), len(b)))
# 配列Qに格子点の値を格納
for j, k in itertools.product(range(res), range(res)):
Q[j, k] = mse(x, y, a[j], b[k])
# 視点の設定
ax.view_init(45, 45)
# 曲面を描画
ax.plot_surface(A, B, Q, color = "red")
plt.show() 
表示されたグラフを見ると、$Q(a,b)$ は谷のような形をしていて、明らかに極小値は１つしかないことがわかります。
平均二乗誤差を最小にするパラメータ
平均二乗誤差 $Q(a,b)$ の定義式を再掲します。
\[Q(a,b)=\frac{1}{N}\sum_{k=0}^{N-1}\{w(x_k)-y_k\}^2=\frac{1}{N}\sum_{k=0}^{N-1}(ax_k+b-y_k)^2\]
目的は最急降下法を使って $Q$ を最小にするような $a,\ b$ を決定することです。
あるステップ数 $t$ におけるパラメータ $a,\ b$ を、それぞれ $a(t),\ b(t)$ と表すことにして、パラメータベクトルを
\[\boldsymbol{p}(t)=\binom{a(t)}{b(t)}\]
と定義すると、ステップ数 $t+1$ におけるパラメータベクトルは
\[\boldsymbol{p}(t+1)=\boldsymbol{p}(t)-\alpha\nabla Q|_{\boldsymbol{p}(t)}\]
と表すことができます。$\nabla Q$ は
\[\nabla Q=\begin{bmatrix}\cfrac{\partial Q}{\partial a}\\
\cfrac{\partial Q}{\partial b}\\\end{bmatrix}\]
で表される $Q$ の勾配 (gradient) です。ステップ $t+1$ におけるパラメータベクトルを成分で表すと、
\[\binom{a(t+1)}{b(t+1)}=\begin{bmatrix}a(t)-\alpha\left.\cfrac{\partial Q}{\partial a}\right|_{a(t),b(t)}\\
b(t)-\alpha\left.\cfrac{\partial Q}{\partial b}\right|_{a(t),b(t)}\\\end{bmatrix}\]
となります。$Q$ を $a$ と $b$ で偏微分すると
\[\begin{align*}
&\frac{\partial Q}{\partial a}=\frac{2}{N}\sum_{k=0}^{N-1}(ax_k+b-y_k)x_k\\[6pt]
&\frac{\partial Q}{\partial b}=\frac{2}{N}\sum_{k=0}^{N-1}(ax_k+b-y_k)\end{align*}\]
となるので、
\[\begin{bmatrix}a(t+1)\\b(t+1)\end{bmatrix}=\begin{bmatrix}
a(t)-\cfrac{2}{N}\displaystyle\sum_{k=0}^{N-1}\{w(x_k)-y_k\}x_k\\
b(t)-\cfrac{2}{N}\displaystyle\sum_{k=0}^{N-1}\{w(x_k)-y_k\}\end{bmatrix}\]
と表されることになります。
# In[3]
# MSEの勾配
def dmse(x, y, a, b):
d_a = 2 * np.mean((a * x + b - y) * x)
d_b = 2 * np.mean(a * x + b - y)
return np.array([d_a, d_b])
# 初期値の設定
a_init, b_init = (-20, 0)
# x,yに初期値を代入
a, b = (a_init, b_init)
# 学習率α
alpha = 0.001
# 収束判定条件
eps = 0.01
# 繰返し最大数
k_max = 100000
for k in range(1, k_max):
a -= alpha * dmse(x, y, a, b)[0]
b -= alpha * dmse(x, y, a, b)[1]
# 勾配の大きさの計算
d = np.sqrt(dmse(x, y, a, b)[0]**2 + dmse(x, y, a, b)[1]**2)
# 勾配の大きさがeps以下になったらループ終了
if d < eps:
break
# 平均2乗誤差を計算
m = mse(x, y, a, b)
# 数値を丸める
a = np.round(a, 5)
b = np.round(b, 5)
d = np.round(d, 6)
m = np.round(m, 3)
print("パラメータベクトル (a, b) = ({0}, {1})".format(a, b))
print("勾配の大きさ d = {}".format(d))
print("平均2乗誤差 mse = {}".format(m))
fig = plt.figure(figsize = (8, 6))
ax = fig.add_subplot(111)
# タイトルを設定
ax.set_title("Age vs. Weight", fontsize = 16)
# 目盛線を表示
ax.grid()
# 軸範囲を設定
ax.set_xlim([0, 25])
ax.set_ylim([0, 90])
# 軸ラベルを設定
ax.set_xlabel("Age", fontsize = 14)
ax.set_ylabel("Weight [kg]", fontsize = 14)
# 身長と体重データの散布図
ax.scatter(x, y, color = "blue")
# 近似直線をプロット
x2 = np.linspace(0, 25, 100)
y2 = a * x2 + b
ax.plot(x2, y2, color = "red")
plt.show()パラメータベクトル (a, b) = (3.01098, 3.93668) 勾配の大きさ d = 0.009997 平均2乗誤差 mse = 45.194
学習率を小さくしているので、環境によっては計算に少し時間がかかるかもしれません。eps 判定で終了していることを確認するために最終地点における勾配の大きさ d を表示させています。繰返し最大数は適当な時間で処理を終了させるための安全弁に過ぎません。勾配の大きさがある程度小さくなっていないと、窪みの底に落ちていないので注意してください。
赤いラインが近似直線です。平均二乗誤差は 45 なので、近似直線と各点との距離の平均は $\sqrt{45}=6.7$ ぐらいになっています。