多次元での線形回帰
前回は2次元での線形回帰について考えました。
今回は、特徴量ベクトルが一般次元の時について考えて行きたいと思います。
数学的な意味
一般次元の場合、線形回帰モデルは次のような式で表すことができる。
y = ω_0 + ω_1x_1 + ω_2x_2 + ・・・ω_dx_d + ε
ここで、
(x_0, x_1, ・・・x_d)^T
は入力変数、
ω_0, ω_1, ・・・ω_d
はパラメータ、yはターゲットになります。またεはノイズであり、必ず発生してしまう誤差ですので、深くは考えません。
ここで、ベクトル
x = (x_0, x_1, ・・・x_d)^T
に対して、要素1を付加させたベクトルをXと考え、ベクトル
ω = (ω_0, ω_1, ・・・ω_d)^T
と定義すると、yは次のように表すことができる。
y = ω^T X
ここで、n×d行の時の行列を考え、線形和を考えます。そうすると、
y =
\begin{bmatrix}
x_{11} & x_{12} & ・・・ x_{1d}\\
x_{21} & x_{22} & ・・・ x_{2d}\\
・・ & ・・ & ・・・・ \\
x_{n1} & x_{n2} & ・・・ x_{nd}\\
\end{bmatrix}
\begin{bmatrix}
ω_1 \\
ω_2 \\
...\\
ω_n \\
\end{bmatrix}
= Xω
と表記することができます。ここで、実際のデータとの誤差を、二乗して和を求め、その和を最小にすることを考えます。最小二乗法ですね。
つまり、以下の式が最小になるようなwを求めます。
E = || y - Xω ||^2
= (y-Xω)^T(y-Xω)
= (y^T - X^Tω^T)(y - Xω)
ここからこの上式を整理すると、
E = y^Ty -2y^TXω + X^Tω^TXω
と変形することができます。
Eを最小にするということなので、wに関する偏微分を考えます。wの偏微分を行いやすくするために、w以外を定数項とみなします。
E = y^Ty -2y^TXω + X^Tω^TXω
= y^Ty + (-2y^TX)ω + (X^TX)ωω^T
= a + bω + ωCω^T
ここで、
a= y^Ty
b= -2y^TX
C= X^TX
と置換をしました。
ここから、偏微分を行っていきます。
\frac{δL}{δω} = \frac{δ}{δω} (a + bω + ωCω^T) \\
\frac{δ}{δω}a + \frac{δ}{δω}b+\frac{δ}{δω} (ωCω^T) \\
= 0 + b^T + ω^T(A+A^T)
となった。ここで、置換をしていたa, b, Cを代入します。
ΔE = (-2y^TX)^T + ω^T(X^TX+(X^TX)^T)
この式をもう少し綺麗にしてみます。式を綺麗にする際に転置の公式を用いていきます。
ΔE = -2yX^T + 2ωX^TX
かなり綺麗になりました。
さて、二次関数を考えた時、最小になる場所について思い出してみます。
二次関数の最小値を考える際、最小となるところは微分をし、その傾きが0になったところが最小でした。
その考え方を今回も同じように使います。そうすると、ΔE=0の時が、最小になるωであるということがわかります。
ωについて解いていきたいと思います。
ΔE = -2yX^T + 2ωX^TX = 0
\\
2yX^T = 2ωX^TX\\
ω = (X^TX)^{-1} yX^T
これで、多次元の時の最小にするωを求めることができました。
実装してみる
これをPythonを使って実装してみます。
import numpy as np
from scipy import linalg
class LinearRegression:
def __init__(self):
self.w = None
def fit(self, X, t):
Xtil = np.c_[np.ones(X.shape[0]), X]
A = np.dot(Xtil.T, Xtil)
b = np.dot(Xtil.T, t)
self.w = linalg.solve(A, b)
def predict(self, X):
if X.ndim == 1:
X= X.reshape(1, -1)
Xtil = np.c_[np.ones(X.shape[0]), X]
return np.dot(Xtil, self.w)
このアルゴリズムではfitメソッドでは訓練データによって学習を行います。
また、predictメソッドを、入力値に対して予測をしたいときに使います。
次に呼び出して予測をしてみます。
import linearreg
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import axes3d
n = 100
scale = 10
np.random.seed(0)
X = np.random.random((n, 2)) * scale
w0 = 4
w1 = 2
w2 = 3
y = w0 + w1 * X[:, 0] + w2 * X[:, 1] + np.random.randn(n)
print(y)
model = linearreg.LinearRegression()
model.fit(X, y)
print("係数:" , model.w )
aaa = model.predict(np.array([1 , 1]))
print("(1, 1)に対する予測値: ",aaa )
xmesh, ymesh = np.meshgrid(np.linspace(0, scale, 20), np.linspace(0, scale, 20))
zmesh = (model.w[0] +model.w[1] * xmesh.ravel() + model.w[2] * ymesh.ravel()).reshape(xmesh.shape)
fig = plt.figure()
ax = fig.add_subplot(111, projection="3d")
ax.scatter(X[:, 0], X[:, 1], y, color = "k")
ax.plot_wireframe(xmesh, ymesh, zmesh, color ="k")
plt.show()
これを実行し、表示させてみます。
これを実行する際、まず要素が乱数であるサイズ100 × 2の行列を生成している。
また、線形和
w_0 + w_1x_1 + w_2x_2
に関して、
w_0 = 4, w_1 =2, w_2 =3
と適当に数字を決め、その線形和に乱数を加え、ノイズを再現しています。
点群が平面に近似されていることがわかりました。
まとめ
線形回帰について理論を理解することができました。
次回は、データセットを用いて、実際に線形回帰を実装してみたいと思います。
この記事を書くに当たって、以下の文献を参考にしました。
参考文献 機械学習のエッセンス