最小二乗法を実装する - メモ的な何か

予測モデルをつくるときに最も簡単な方法としては「最小二乗法による線形モデル」と「k近傍モデル」である.
この2つは統計学や機械学習を学んでいない人でも割と思いつきそうな方法である.
線形モデルとはパラメータに対して線形なモデルという意味である.
例えば, y = ax + b という1次関数を中学で学んだ.
今, データxが観測されたときに, yを予測するというタスクを行う.
このときa(傾き), b(切片)の値によって予想するべき値yが異なってくる.
すなわちa, bが決定されれば全てのxの入力に対してyという出力ができるようになる.
このようなa, b をパラメータという.
これを多次元に拡張すると
$y = w_0 + w_1x_1 + w_2x_2 + \ldots w_nx_n$
これは
$y = w_01 + w_1x_1 + w_2x_2 + \ldots w_px_p$
と見て
$y = \boldsymbol{x}^T \boldsymbol{w}$
とベクトル表記で書くことができる.
これは, p+1 次元の一つのデータに対する予測の式なので, n個のデータを同時に扱う場合には
$\boldsymbol{y} = \boldsymbol{X}^T \boldsymbol{w}$
とする. $\boldsymbol{X}$ は(p+1, n)行列であり転置すると(n, p+1)行列, パラメータベクトル $\boldsymbol{w}$ は(p+1)ベクトルである. ((p+1, 1)行列という表現もできる)
この線形モデルを, 最小二乗法によりフィットする. (予測と実際の誤差を最小にする)
$min: RSS(\boldsymbol{w}) = (\boldsymbol{y} - \boldsymbol{X}\boldsymbol{w})^T(\boldsymbol{y}-\boldsymbol{X}\boldsymbol{w})$
$\begin{align*} RSS(\boldsymbol{w}) &= (\boldsymbol{y}^T-(\boldsymbol{X}\boldsymbol{w})^T)(\boldsymbol{y}-\boldsymbol{X}\boldsymbol{w})\\ &= \boldsymbol{y}^T\boldsymbol{y}-\boldsymbol{y}^T\boldsymbol{X}\boldsymbol{w}-(\boldsymbol{X}\boldsymbol{w})^T\boldsymbol{y}+(\boldsymbol{X}\boldsymbol{w})^T(\boldsymbol{X}\boldsymbol{w})\\ &= \boldsymbol{y}^T\boldsymbol{y}-(\boldsymbol{y}^T\boldsymbol{X})\boldsymbol{w}-\boldsymbol{w}^T(\boldsymbol{X}^T\boldsymbol{y})+\boldsymbol{w}^T(\boldsymbol{X}^T\boldsymbol{X})\boldsymbol{w} \end{align*}$
この式はパラメータ $\boldsymbol{w}$ に関して二次式であるので最小値が存在する.
パラメータ $\boldsymbol{w}$ で微分して $\boldsymbol{0}$ になるパラメータ $\boldsymbol{w}$ を求める.
$\begin{align*} &\boldsymbol{0} - \frac{\partial (\boldsymbol{y}^T\boldsymbol{X})\boldsymbol{w}}{\partial \boldsymbol{w}} - \frac{\boldsymbol{w}^T(\boldsymbol{X}^T\boldsymbol{y})}{\partial \boldsymbol{w}} + \frac{\partial \boldsymbol{w}^T(\boldsymbol{X}^T\boldsymbol{X})\boldsymbol{w}}{\partial \boldsymbol{w}} = \boldsymbol{0}\\ &-(\boldsymbol{y}^T\boldsymbol{X})^T-(\boldsymbol{X}^T\boldsymbol{y}) + (\boldsymbol{X}^T\boldsymbol{X} + (\boldsymbol{X}^T\boldsymbol{X})^T)\boldsymbol{w} = \boldsymbol{0}\\ &-\boldsymbol{X}^T\boldsymbol{w}-\boldsymbol{X}^T\boldsymbol{y} + (\boldsymbol{X}^T\boldsymbol{X} + \boldsymbol{X}^T\boldsymbol{X})\boldsymbol{w} = \boldsymbol{0}\\ &-2\boldsymbol{X}^T\boldsymbol{y} + 2\boldsymbol{X}^T\boldsymbol{X}\boldsymbol{w} = \boldsymbol{0}\\ &2\boldsymbol{X}^T\boldsymbol{X}\boldsymbol{w} = 2\boldsymbol{X}^T\boldsymbol{y}\\ &\boldsymbol{w} = (\boldsymbol{X}^T\boldsymbol{X})^{-1}\boldsymbol{X}^T\boldsymbol{y} \end{align*}$

ベクトルに対する微分の以下の式は覚えておいたほうが良い.
$\begin{align*} \frac{\partial A^T \boldsymbol{x}}{\partial \boldsymbol{x}} = \frac{\partial \boldsymbol{x}^T A}{\partial \boldsymbol{x}} &= A\\ \frac{\partial \boldsymbol{x}^TA\boldsymbol{x}}{\partial \boldsymbol{x}} &= (A+A^T)\boldsymbol{x}\\ if\ A\ is\ symmetric\ matrix\\ &= 2A\boldsymbol{x} \end{align*}$

線形モデルの最小二乗法について, パラメータの解が
$\begin{align*} \boldsymbol{w} = (\boldsymbol{X}^T\boldsymbol{X})^{-1}\boldsymbol{X}^T\boldsymbol{y} \end{align*}$
で得られることが分かった.
これをPythonによって実装してみる. 今回は一次関数のパラメータa,b を求め直線フィッティングを行う.
まずデータを生成する. データ生成のアルゴリズムは以下の通りに行った.

a, b, σを決定する
[-10, 10]の範囲でデータxをランダムに生成する
決定したa, bからt=ax+bを計算する
t=ax+b+εを計算する. ただしε~N(0, σ^2)

### 単純な線形モデルy = ax + b についてのパラメータa,b を求める問題におけるデータを生成する
### 
### a, b, σの値を決定する
### [-x_min, x_max]のデータxをランダムにN個用意 
### 各データに対して, N(ax+b, σ^2)に従うガウス分布からtを生成
### 

import numpy as np
from numpy.random import seed
import matplotlib.pyplot as plt

seed(0)

a, b, sigma = 2, 3, 3
N = 100
x_min, x_max = -10, 10

# データ生成
X = (x_max-x_min) * np.random.rand(N) + x_min
Epsiron = np.random.normal(0, sigma**2, N)# N(0, sigma^2)のN個の正規乱数
T = (a*X + b) + Epsiron

# データをファイルに出力
f = open("data00.csv", "w")
for (x, t) in zip(X, T):
    f.write("{}, {}\n".format(x,t))
f.close()

# データをプロット
plt.scatter(X, T)
plt.show()

以下の図のようになった.

f:id:umashika5555:20181221091859p:plain — データの様子

このデータから,
$\begin{align*} \boldsymbol{w} = (\boldsymbol{X}^T\boldsymbol{X})^{-1}\boldsymbol{X}^T\boldsymbol{y} \end{align*}$
に従ってパラメータを決定する.

import numpy as np
import matplotlib.pyplot as plt

# データの読み込み
name_data = "data00"
path_data = name_data + ".csv"
data = np.loadtxt(path_data,delimiter=",")
x, t = data[:,0], data[:,1]
p = x.shape[0]
ones = np.ones(p)
x_ = np.c_[ones, x]

# パラメータベクトルの計算
A = np.dot(x_.T, x_)
w = np.dot(np.dot(np.linalg.inv(A), x_.T), t)
print(w)

# データのプロット
plt.scatter(x, t)

# 回帰曲線の描画
b, a = w
xx = np.linspace(np.min(x)-2, np.max(x)+2, 100)
yy = a * xx + b
plt.plot(xx, yy, color="red")

plt.xlabel("x")
plt.ylabel("y")
plt.title("a=2, b=3, σ=2")
plt.savefig(name_data.format(a,b))
plt.show()