#はじめに
重回帰分析モデルの自力実装にチャレンジしました。
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変数として採用したデータ間で単位が異なるので、本来はどちらも標準化してから重回帰のモデルを作った方が良さそうです。精度の上げ方などは別途勉強していきたいと思います。
##Next
次回は主成分分析の理論理解と自力実装に取り組んでいきます。