実装

前回は数学的な議論に徹してしまい、自分でも何を目的に書いていたのかを見失ってしまいました。今回は、Pythonで重回帰の実装します。

自力で作る

自力で作ることは、Pythonの文法やNumpyに慣れるために必要不可欠です。その1の実装を参考にするとより理解が深まるかと思います。

重みは$\beta = (X^TX)^{-1}X^Ty$と与えられますが、これは正規方程式$X^Ty = X^TX\beta$の解です。$y=X\beta$を解かないのは、一般的に$X$は正方行列ではないので$X$が逆行列を持たないことがあるからです。$X^TX$が正則でないと、エラーが出てしまいますが、ひとまずそのまま$\beta = (X^TX)^{-1}X^Ty$で計算させることにしましょう。

model = LinearRegression()
model.fit(X_train, y)
model.predict(X_test)

を骨格に作ってみましょう。行列と行列、または行列とベクトルの積の計算はNumpyのnumpy.dot(A, B)を使用すると大変楽です。

import numpy as np
import matplotlib.pyplot as plt

はお馴染みで、回帰モデルのクラスを作りましょう。
Pythonの文法に慣れていないせいで、なかなか苦労しました。
なお、簡単のため、バイアスは重みを表す行列の左の列に加えています。
それに伴って、説明変数の一番左に新しく要素1を追加しています。

class LinearRegression():
        
    def fit(self, X, y):
        
        X_ = np.concatenate([np.ones(X.shape[0]).reshape(-1, 1), X], axis=1)
        # X_ = np.c_[np.ones(X.shape[0]), X]とも書ける
        self.beta = np.linalg.inv(np.dot(X_.T, X_)).dot(X_.T).dot(y)
        
    def predict(self, X):
        X_ = np.concatenate([np.ones(X.shape[0]).reshape(-1, 1), X], axis=1)
        # X_ = np.c_[np.ones(X.shape[0]), X]とも書ける
        return np.dot(X_, self.beta)

例として、$y=2x_1+3x_2+4$に近い分布をもつサンプルを使用。

n = 100   # 100個のサンプルを用意
np.random.seed(0)
X_train = (np.random.random((n, 2)) - 0.5)  * 20 #要素が乱数となるようにn＊2の行列を生成（範囲を-10から10までとした）
y = 2*X_train[:, 0] + 3*X_train[:, 1] + 4 + np.random.randn(n) * 0.1

model = LinearRegression()
model.fit(X_train, y)
print("beta={}".format(model.beta))

これを実行すると

beta=[3.9916141 1.9978685 3.00014788]

と出力されます。beta[0]がバイアスです。よい近似になっていることが分かります。

可視化しましょう。手始めにX_trainを散布図に表してみます。

from mpl_toolkits.mplot3d import Axes3D

fig = plt.figure()
ax = Axes3D(fig)

ax.plot(X_train[:, 0], X_train[:, 1], y, marker="o", linestyle="None")

これを踏まえて、予測された平面を図示してみましょう。

fig = plt.figure()
ax = Axes3D(fig)

ax.plot(X_train[:, 0], X_train[:, 1], y, marker="o", linestyle="None")

xmesh, ymesh = np.meshgrid(np.linspace(-10, 10, 20), np.linspace(-10, 10, 20))
zmesh = (model.beta[1] * xmesh.ravel() + model.beta[2] * ymesh.ravel() + model.beta[0]).reshape(xmesh.shape)
ax.plot_wireframe(xmesh, ymesh, zmesh, color="r")

scikit-learnのライブラリ（LinearRegression）で実装する

こちらは実務向けです。一気に実装します！

from sklearn.linear_model import LinearRegression

n = 100   # 100個のサンプルを用意
np.random.seed(0)
X_train = (np.random.random((n, 2)) - 0.5)  * 20 #要素が乱数となるようにn＊2の行列を生成（範囲を-10から10までとした）
y = 2*X_train[:, 0] + 3*X_train[:, 1] + 4 + np.random.randn(n) * 0.1

model = LinearRegression()
model.fit(X_train, y)
print("coef_, intercept_ = {}, {}".format(model.coef_, model.intercept_))

X_test = np.c_[np.linspace(-10, 10, 20), np.linspace(-10, 10, 20)]
model.predict(X_test)

fig = plt.figure()
ax = Axes3D(fig)
ax.plot(X_train[:, 0], X_train[:, 1], y, marker="o", linestyle="None")
xmesh, ymesh = np.meshgrid(X_test[:, 0], X_test[:, 1])
zmesh = (model.coef_[0] * xmesh.ravel() + model.coef_[1] * ymesh.ravel() + model.intercept_).reshape(xmesh.shape)
ax.plot_wireframe(xmesh, ymesh, zmesh, color="r")

ちゃんと同じ図が得られましたか？

残差をプロットしてみましょう。

# 残差プロット
y_train_pred = model.predict(X_train)
plt.scatter(y_train_pred, y_train_pred - y, marker="o", label="training data")
plt.scatter(y_test_pred, y_test_pred - (2*X_test[:, 0]+3*X_test[:, 1]+4), marker="^", label="test data")
plt.xlabel("predicted values")
plt.ylabel("residuals")
plt.hlines(y=0, xmin=-40, xmax=55, color="k")
plt.xlim([-40, 55])

重回帰についてはもっと色々勉強すべきことが多いですが、実務はこのぐらいで十分でしょう。
次回はRidge回帰、Lasso回帰の予定です。

参考文献

加藤公一『機械学習のエッセンス』SB Creative、2019年
S. Raschka, V. Mirjalili『Python機械学習プログラミング[第2版]』インプレス、2020年