最小二乗法の適用例

  • $x$と$y$は本質的には線形関係を持っている（$y=ax+b$）．しかし，$x$に対する$y$の値を計測する時に必ず誤差が生じる．
  • $(x,y)=(1,2),(2,4),(3,5),(4,7)$の観測が得られた際に，最小二乗法にもとづいて $a, b$ を求めよ．

解法１

$y = ax + b = w^T x$と表す。ここで$w = \begin{pmatrix}a \\ b \end{pmatrix}$ である（つまり、$w^T = (a \ b)$。

問題設定から$\boldsymbol{x}= \begin{pmatrix}1 \\ 2 \\ 3 \\ 4\end{pmatrix}$ 　$\boldsymbol{y}= \begin{pmatrix}2 \\4 \\5 \\7\end{pmatrix}$

教科書p.153の(11.11)$\frac{\partial E}{\partial w_j} = -2 \sum_i x_{i,j}(y_i - {w^*}^T x_i)$ を解いてみる。

$E(w) = \sum_i(y_i - w^Tx_i)^T(y^i - w^Tx_i)$に問題の設定を代入すると $E(w) = (2-(a+b)*1)^2 + (4-(2a+b))^2 + (5-(3a+b))^2 + (7-(4a+b))^2$ であり、 これを最小とする $w^T=(a, b)$を求める

$\frac{\partial E}{\partial a} =-2(2-a-b) - 4(4-2a-b) -6(5-3a-b) - 8(7-4a-b) = 60a + 20b - 106$
極値では $= 0$ となるので $30a + 10b = 53$

$\frac{\partial E}{\partial b} =-2(2-a-b) - 2(4-2a-b) -2(5-3a-b) - 2(7-4a-b) = 20a + 8b - 36 = 0$
書き換えれば $5a + 2b = 9$

この２つの式を解いて$a,b$の値を求める： $a = 1.6, b = 0.5$

解法２

教科書 p.153 の(11.15)を用いる。問題設定は解法1と同じである。 ただこのままだと値がスカラーなので、ダミーの定数をいれてベクトルとする。 つまり $X = \begin{pmatrix}1 & 1 & 1 & 1\\ 1 & 2 & 3 & 4\end{pmatrix}, y=\begin{pmatrix} 1 & 1 & 1 & 1\\ 2 & 4 & 5 & 7\end{pmatrix}$と表す。

In [1]:
import numpy as np
X=np.array([[1,1],[1,2],[1,3],[1,4]]).T
y=np.array([[1,2],[1,4],[1,5],[1,7]]).T
In [2]:
print(X)
print(y)

[[1 1 1 1]
 [1 2 3 4]]
[[1 1 1 1]
 [2 4 5 7]]
In [3]:
XX= X.dot(X.T)
XX
Out[3]:
array([[ 4, 10],
       [10, 30]])
In [4]:
XX_inv=np.linalg.inv(XX)
XX_inv
Out[4]:
array([[ 1.5, -0.5],
       [-0.5,  0.2]])
In [5]:
Xy = X.dot(y.T)
Xy
Out[5]:
array([[ 4, 18],
       [10, 53]])
In [6]:
w_ast=XX_inv.dot(Xy)
In [7]:
w_ast
Out[7]:
array([[1. , 0.5],
       [0. , 1.6]])
In [8]:
w_ast.T.dot(X)
Out[8]:
array([[1. , 1. , 1. , 1. ],
       [2.1, 3.7, 5.3, 6.9]])

解法3: scikit-learn

In [2]:
# scikit-learnを用いる
import numpy as np
from sklearn.linear_model import LinearRegression
X=np.array([[x+1] for x in range(4)])
y=np.array([2,4,5,7])
model=LinearRegression()
model.fit(X,y)
print("slope="  ,model.coef_[0], "intercept= ", model.intercept_)

slope= 1.6000000000000003 intercept=  0.4999999999999991
In [7]:
%matplotlib inline
import matplotlib.pyplot as plt
a=model.coef_[0]
b=model.intercept_
plt.plot(X,y,"go")
plt.plot(X,[a*x + b for x in X], "b-")
Out[7]:
[<matplotlib.lines.Line2D at 0x2937d3f9cf8>]