最小二乗法をフルスクラッチで実装する

ココで言うフルスクラッチは行列演算以外である.
行列演算は本質ではないので, 行列演算ライブラリは使うことにする.

前回の記事で線形モデルの最小二乗法について, パラメータの解が

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

  1. a, b, σを決定する
  2. [-10, 10]の範囲でデータxをランダムに生成する
  3. 決定したa, bからt=ax+bを計算する
  4. 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()

f:id:umashika5555:20181221092149p:plain
data00
上図のように実際に目で観ても綺麗にフィットできていることがわかる.


次回はこの線形モデルを更に発展させ, 基底関数, 正則化あたりをやっていきたい.