線形回帰モデル¶
Pythonの機械学習用ライブラリscikit-learnを使って，線形回帰モデルを使って簡単な回帰問題にチャレンジしてみましょう.
0.ライブラリのインポート¶
import pandas as pd
import numpy as np
import sklearn
import seaborn as sns
import matplotlib
import matplotlib.pyplot as plt
%matplotlib inline
np.set_printoptions(precision=4)
print("numpy :", np.__version__)
print("pandas :", pd.__version__)
print("sklearn :", sklearn.__version__)
print("seaborn :", sns.__version__)
print("matplotlib :", matplotlib.__version__)
1. データの読込・整形¶
sklearn.datasetsからBostonデータセットを読み込みましょう．
# make data samples
from sklearn.datasets import load_boston
boston = load_boston()
次に，pandas DataFrame()クラスのインスタンスとして，変数df_data, df_target, dfを定義します．
df_data = pd.DataFrame(boston.data, columns = boston.feature_names)
df_target = pd.DataFrame(boston.target, columns = ['TARGET'])
df = pd.concat([df_target, df_data], axis = 1)
df = df.loc[:, ['TARGET', 'RM', 'DIS', 'INDUS', 'AGE', 'LSTAT', 'CRIM']]
df.head(10)
データの要約統計量(サンプル数, 平均, 標準偏差, 四分位数, 中央値, 最小値, 最大値など)をみましょう．
df.describe().T
データの共分散行列を描画します．
対角成分は自分との共分散(相関)を表すため常に1.0となります．
df.corr()
データの散布図行列を描画します．
相関が大きい説明変数のペアについては, 多重共線性を考えるべきです.
sns.set()
sns.pairplot(df, height=2.0, diag_kind='hist', markers='+')
plt.show()
2. データの分割¶
pandas DataFrame()クラスの変数dfから，説明変数と目的変数に相当するデータをそれぞれ取り出し，numpy.ndarray()クラスの変数x, yへ格納します．
x = df.iloc[:, 1:] # set x to explanatory variable
y = df.iloc[:, 0] # set y to response variable
全データをtrainデータとtestデータに分割します．
すなわち，変数xをx_trainとx_testに，
変数yをy_trainとy_testに分けます．
# split data by Hold-out-method
from sklearn.model_selection import train_test_split
x_train, x_test, y_train, y_test = train_test_split(x, y, test_size = 0.3, random_state = 0)
print()で配列の形状を確認してみましょう．
print("x_train: ", x_train.shape)
print("y_train: ", y_train.shape)
print("x_test: ", x_test.shape)
print("y_test: ", y_test.shape)
- x_train: 6次元データが354コ格納されている．
- y_train: 1次元データが354コ格納されている．
- x_test: 6次元データが152コ格納されている．
- y_test: 1次元データが152コ格納されている．
3. モデルの作成¶
ボストン市の住宅価格TARGETを目的変数にして線形回帰モデル$f$を作ります.$\beta$はモデルパラメータです．
$$ \begin{align} y &= f({\bf x} ; \beta) + u \\ &= f(x_1, \dots, x_6 ; \beta) + u \\ &= \beta_0 + \beta_1 x_1 + \cdots + \beta_6 x_6 + u \\ \\ u &\sim \mathcal{N}(0, \sigma^2) \end{align} $$- 目的変数 (1コ)
- $y~~$ : TARGET - 住宅価格の中央値（単位: 千ドル）
- 説明変数 (6コ)
- $x_1$ : RM - 住居の平均部屋数
- $x_2$ : DIS - 5つのボストン市の雇用施設からの重み付き距離
- $x_3$ : INDUS - 非小売業の土地割合（単位: ％）
- $x_4$ : AGE - 1940年より前に建てられた物件割合（単位: ％）
- $x_5$ : LSTAT - 低所得労働者の割合（単位: ％）
- $x_6$ : CRIM - 人口あたり犯罪発生率（単位: ％）
- 誤差項
- $u$ ： 平均0の正規分布に従うと仮定．
# Linear Regression
from sklearn.linear_model import LinearRegression
lr = LinearRegression()
4. モデルへデータを適合させる¶
データ$\mathcal{D} := \left\{ \left( \begin{array}{c} {\bf x}_1 \\ y_1 \end{array} \right), \dots, \left( \begin{array}{c} {\bf x}_n \\ y_n \end{array} \right) \right\}$を用いて計算される平均二乗誤差(Mean Squared Error, MSE)をパラメータ$\beta$の目的関数$J$とすれば，$J$を最小化するようなパラメータ$\beta^*$を求めることができます．
$$ \begin{align} \beta^* &= \underset{\bf \beta}{\rm argmin} ~ J(\beta) \\ J(\beta) &:= \sum_{i=1}^{n} {( y_i - f({\bf x}_i ; \beta) )}^2 \end{align} $$実装上は，sklearn.linear_model.LinearRegression()クラスのインスタンスlrを作成し，trainデータx_train,y_trainを与えることで最適なモデルパラメータ$\beta^* = [\beta_0^* \dots,\beta_6^*]$を求めます．
lr.fit(x_train, y_train)
trainデータから推定された最適なモデルパラメータ$\beta^* = [\beta_0^* \dots,\beta_6^*]$の値を確認してみます．
print(' beta0 : {:.4e}'.format(lr.intercept_))
print(' beta1 : {:.4e}'.format(lr.coef_[0]))
print(' beta2 : {:.4e}'.format(lr.coef_[1]))
print(' beta3 : {:.4e}'.format(lr.coef_[2]))
print(' beta4 : {:.4e}'.format(lr.coef_[3]))
print(' beta5 : {:.4e}'.format(lr.coef_[4]))
print(' beta6 : {:.4e}'.format(lr.coef_[5]))
5. モデルの評価¶
モデルlrに，入力としてx_train,x_testを与えて，y_train,y_testの値をそれぞれ予測させてみましょう．
# predictions
y_train_pred = lr.predict(x_train)
y_test_pred = lr.predict(x_test)
評価指標(Metrics)として，決定係数(R^2)，平均二乗誤差(MSE), 平均絶対誤差(MAE)の値を計算してみましょう．
# model evaluation
from sklearn.metrics import mean_absolute_error, mean_squared_error
print(' R^2 (train) : {:.4f}'.format(lr.score(x_train, y_train)))
print(' R^2(test) : {:.4f}'.format(lr.score(x_test, y_test)))
print('')
print(' MSE (train) : {:.4e}'.format(mean_squared_error(y_train, y_train_pred)))
print(' MSE (test) : {:.4e}'.format(mean_squared_error(y_test, y_test_pred)))
print('')
print(' MAE (train) : {:.4e}'.format(mean_absolute_error(y_train, y_train_pred)))
print(' MAE (test) : {:.4e}'.format(mean_absolute_error(y_test, y_test_pred)))
目的関数$J$をみれば明らかですが，trainデータに適合する(fitする)ようにモデルlrは設計されています．
（モデルにとって）未知のデータであるtestデータに対する評価指標の値は，trainデータに対するものと比べて低くなります．
モデルがtrainデータに対して過度に適合した結果，testデータに対する性能が悪化する現象を「過学習(Over fitting)」といいます．