はじめに

重回帰分析モデルの自力実装にチャレンジしました。
sklearnの結果と同様のものを導き出すことをゴールとしています。
前回の投稿はこちら（重回帰分析を理解する（理論編））です。数学的な理論の方はそちらにまとめています。

参考

多変量解析法入門　著者; 永田靖, 棟近雅彦　出版社; サイエンス社
多変量解析のはなし　著者; 有馬哲, 石村貞夫　出版社; 東京図書
Chainer Tutorial-Chainerのチュートリアルサイト
【キカガク流】人工知能・機械学習脱ブラックボックス講座 - 中級編 - Udemyのオンライン講座

今回使用するデータ

今回はsklearnの中にデータセットとして用意されている、ボストン住宅価格のデータセットを用いて検証を行います。
まず下記のように、データを読み込みます。


from sklearn.datasets import load_boston
import pandas as pd

#データの読み込み
boston = load_boston()

#一旦、pandasのデータフレーム形式に変換
df = pd.DataFrame(boston.data, columns=boston.feature_names)

#目的変数（予想したい値）を取得
target = boston.target

df['target'] = target

各説明変数のデータのカラムはこのようになっています。
今回は「INDUS」と「CRIM」の2変数を用いて、住宅価格を予測するモデルを作成して検証しようと思います。

カラム	説明
CRIM	町ごとの一人当たりの犯罪率
ZN	宅地の比率が25,000平方フィートを超える敷地に区画されている。
INDUS	町当たりの非小売業エーカーの割合
CHAS	チャーリーズ川ダミー変数（川の境界にある場合は1、それ以外の場合は0）
NOX	一酸化窒素濃度（1000万分の1）
RM	1住戸あたりの平均部屋数
AGE	1940年以前に建設された所有占有ユニットの年齢比率
DIS	5つのボストンの雇用センターまでの加重距離
RAD	ラジアルハイウェイへのアクセス可能性の指標
TAX	10,000ドルあたりの税全額固定資産税率
PTRATIO	生徒教師の比率
B	町における黒人の割合
LSTAT	人口当たり地位が低い率
MEDV	1000ドルでの所有者居住住宅の中央値

理論編の復習

\begin{split}\begin{aligned}
\boldsymbol{y}=
\begin{bmatrix}
y_{1} \\
y_{2} \\
\vdots \\
y_{n}
\end{bmatrix} \\
\end{aligned}\end{split}

\begin{split}\begin{aligned}
\boldsymbol{\hat{y}}=
\begin{bmatrix}
\hat{y}_{1} \\
\hat{y}_{2} \\
\vdots \\
\hat{y}_{n}
\end{bmatrix} \\
\end{aligned}\end{split}

\begin{split}\begin{aligned}
\boldsymbol{w}=
\begin{bmatrix}
w_{0} \\
w_{1} \\
\vdots \\
w_{n}
\end{bmatrix} \\
\end{aligned}\end{split}

\begin{split}\begin{aligned}
X=
\begin{bmatrix}
x_{10} & x_{11} & x_{12} & \cdots & x_{1m} \\
x_{20} & x_{21} & x_{22} & \cdots & x_{1m} \\
\vdots \\
x_{n0} & x_{n1} & x_{n2} & \cdots & x_{nm}
\end{bmatrix} \\
\end{aligned}\end{split}

$\boldsymbol{y}$は目的変数の実測値をベクトル化したもの
$\boldsymbol{\hat{y}}$は回帰式によって導き出された予測値をベクトル化したもの
$\boldsymbol{w}$は重回帰式を作成した際の回帰係数をベクトル化したもの
$X$はサンプル数$n$個、変数の数$m$個の説明変数の実測値を行列化したもの

理論編では色々式変形を記載しましたが、結局のところ上記のような定義のもとで下記式で求めたい回帰係数を導出できるということでした。


\begin{split}\begin{aligned}
{\boldsymbol{w}}
&=({X}^{T}{X})^{-1}X^T{\boldsymbol{y}}  \\

\end{aligned}\end{split}

こちらをPythonで実装していきます。

重回帰モデルの実装

下記が重回帰モデルの自力実装を行った結果です。

import numpy as np

class LinearReg_made:

    def __init__(self):
        self.coef_ = None
        self.intercept_ = None

    def fit(self, X, y):
        #切片の計算を含めるために、説明変数の行列の1行目に全て値が1の列を追加
        X = np.insert(X, 0, 1, axis=1)
        #ここで上記数式の計算を実施
        temp = np.dot(np.dot(np.linalg.inv(np.dot(X.T, X)), X.T),y)
        #これは回帰係数の値
        self.coef_ = temp[1:]
        #これは切片の値
        self.intercept_ = temp[0]

    def predict(self, X):
        #重回帰モデルによる予測値を返す
        return (np.dot(X, self.coef_) + self.intercept_)

    def score(self, X, y):
        #残差の分散を求める
        u = ((y - self.predict(X))**2).sum()
        #実測値の分散を求める
        v = ((y - y.mean())**2).sum()
        #決定係数を返す       
        return 1 - u/v

検証

上記の自力実装のモデルと、sklearnの結果が一致しているか検証します。

sklearnの結果

from sklearn import linear_model

clf = linear_model.LinearRegression()

clf.fit(df[['INDUS', 'CRIM']], df['target'])

print(clf.coef_)
print(clf.intercept_)
print(clf.score(df_x[['INDUS', 'CRIM']], df_x['target']))

出力は上から、「回帰係数」「切片」「決定係数」です。

[-0.52335108 -0.24547819]
29.248292706457804
0.27798487095009117

自力実装の結果


#自力実装のモデルがNumpyの形式にしか対応していないので、そのように変換。
X = df[['INDUS', 'CRIM']].values
y = df['target'].values

linear = LinearReg_made()

linear.fit(X,y)

print(linear.coef_)
print(linear.intercept_)
print(linear.score(X,y))

出力はこちらです。数値は小数点の計算領域の差のみで、ほぼ完全一致していることがわかります。
※出力後のデータ形式が異なるので、その部分は異なっています。

[[-0.52335108]
 [-0.24547819]]
[29.24829271]
0.27798487095009117

決定係数の値がかなり低いので、モデルの精度としてはかなりイケてない感じになっています。
2変数として採用したデータ間で単位が異なるので、本来はどちらも標準化してから重回帰のモデルを作った方が良さそうです。精度の上げ方などは別途勉強していきたいと思います。

次回は主成分分析の理論理解と自力実装に取り組んでいきます。

重回帰分析を理解する（自力実装編）