線形回帰

線形回帰とは

説明変数(予測に利用するデータ)を使って、目的変数(予測するデータ)を予測する線形関数(予測モデル)を作る手法。
予測するためには予測モデルを作るための学習が必要。=教師あり学習

下の図の横軸が説明変数、縦軸が目的変数、赤点が得られたデータだとすると、青線が回帰式(学習から得られた予測モデル)となる。
この回帰式が作れたら、未知の説明変数が入ってきたら、その時の目的変数がどうなりそうか予測できるよねって感じ。

回帰式の作成の仕方は、最小二乗法(実際の目的変数の値と予測値との差の二乗（二乗誤差）を最小にする方法)を使う。

上図でいうと、赤点から青線へ向かってx軸の垂線を引くと、その長さが誤差となる。これをプラスマイナスの影響を考慮して二乗し、全て足した値を最小にするイメージ。

利用データ用意

ボストンの住宅価格のデータを使う。

{get_boston_house_price_dataset.py}

import numpy as np
import pandas as pd
from pandas import Series, DataFrame

import matplotlib.pyplot as plt
import seaborn as sns
sns.set_style('whitegrid')
%matplotlib inline

from sklearn.datasets import load_boston
boston = load_boston() # ボストンの住宅価格データの取得

# データフレームに変換
boston_df = DataFrame(boston.data)
boston_df.columns = boston.feature_names
boston_df['Price'] = boston.target

データの概要

{describe_boston_house_price_dataset.py}

boston_df[['RM', 'Price']].describe()

RMは、部屋の数。Priceが価格。
価格(目的変数)を予測するために、部屋の数を説明変数にしてみることとする。

簡単に部屋の数と価格の関係を見てみる。

{scatter_plot_boston_house_price_dataset.py}

# RM：部屋の数を横軸、Price:価格を縦軸に散布図を書いてみる
plt.scatter(boston_df['RM'], boston_df['Price'])
plt.xlabel('Number of rooms')
plt.ylabel('Price($1,000)')
plt.title('Scatter Plot')

散布図は上の通り。正の相関(部屋の数が増えれば、価格も増える関係)がありそう。

そこで相関係数を出してみる。

{calculate_correration.py}

np.corrcoef (boston_df['RM'], boston_df['Price'])[0,1]
> 0.69535994707153925

相関係数0.695 強め。
なので、説明変数にRMを使うのは、筋は悪くはなさそう。

やってみる（単回帰）

１．seabornのlmplotで散布図と同時に回帰式を見てみる

{make_lmplot.py}

sns.lmplot('RM', 'Price', data = boston_df) 

めっちゃ簡単。seabornのlmplotのオプション:dataにDataFrameを渡し、X軸とY軸にDataFrameのカラム名を指定してあげるだけで線を作ってくれる。

２．numpyでやってみる

{calcurate_single_regression.py}

# numpyに喰わせるためにデータを変形　【※1】
X = np.vstack(boston_df.RM) 
X = np.array([[value, 1] for value in X])
# 目的変数を用意
Y = boston_df.Price

# 線形回帰(np.linalg.lstsq)を実行してa:傾き、b:切片を取得。【※2】
a, b = np.linalg.lstsq(X, Y)[0] 

# 散布図を書いて、そこに、得られた結果から回帰式を追記してみる。
plt.plot(boston_df.RM, boston_df.Price, 'o')
x = boston_df.RM
# 回帰式は、y = a * x + bで表現できるので、x軸にxを、y軸に a * x + bを指定。
plt.plot(x, a * x + b, 'r')
# xとy軸のラベル付与
plt.xlabel('Number of Room')
plt.ylabel('Price($1,000)')

numpyでの線形回帰の結果の図。赤線が得られた回帰式。

やってみる（重回帰）

sckit-learnを利用

単回帰が説明変数1つを利用するのに対し、重回帰は複数の説明変数を利用する。
boston_dfには、RM:部屋の数を含め13個の変数が存在するので、これらを使って価格を予測するモデルを作ってみる。

予測モデルを作る際は、学習(training)用とテスト(test)用に分ける方法をとる。
分ける理由は、「学習データだけで精度を見ると、未知のデータの予測に使えるかの判断ができない」から。

{do_multiple_regression.py}

import sklearn
from sklearn.linear_model import LinearRegression # 線形回帰用

# 利用データの確認(CRIM〜Priceまでの14個のデータがある。Priceは目的変数なので、それ以外の13個を説明変数として利用してみる)
boston_df.info()
> <class 'pandas.core.frame.DataFrame'>
> Int64Index: 506 entries, 0 to 505
> Data columns (total 14 columns):
> CRIM       506 non-null float64
> ZN         506 non-null float64
> INDUS      506 non-null float64
> CHAS       506 non-null float64
> NOX        506 non-null float64
> RM         506 non-null float64
> AGE        506 non-null float64
> DIS        506 non-null float64
> RAD        506 non-null float64
> TAX        506 non-null float64
> PTRATIO    506 non-null float64
> B          506 non-null float64
> LSTAT      506 non-null float64
> Price      506 non-null float64
> dtypes: float64(14)

# Price以外の説明変数用DataFrame作成
X_multi = boston_df.drop('Price', 1) 
X_multi.shape
> (506, 13)

# データをtrain用とtest用に分ける。
# sklearn.cross_validation.train_test_split(X, Y) でXを説明変数、Yを目的変数で指定するとサクッとわけてくれる。
X_train, X_test, Y_train, Y_test = sklearn.cross_validation.train_test_split(X_multi, boston_df.Price) 

# 分割されたデータの数を確認
print(X_train.shape, X_test.shape, Y_train.shape, Y_test.shape)
> (379, 13) (127, 13) (379,) (127,)
# train用に379個、test用に127個に分割したことがわかる

# インスタンス生成
lreg = LinearRegression()
# 重回帰のモデル作成実施
lreg.fit(X_train, Y_train)
# 学習データとテストデータに対する予測実施
pred_train = lreg.predict(X_train)
pred_test = lreg.predict(X_test) 

# 学習データとテストデータの平均二乗誤差を確認[average(実際の価格 - 予測値]^2)]
np.mean((Y_train - pred_train) ** 2)
> 20.5592370941859
np.mean((Y_test - pred_test) ** 2) 
> 28.169312238554202

# 実際の正解と予測の残差をグラフで確認する。
# 横軸にtrainの結果、縦軸に誤差を表示
train = plt.scatter(pred_train, (pred_train - Y_train), c= 'b', alpha=0.5)  
# 横軸にtestの結果、縦軸に誤差を表示
test = plt.scatter(pred_test, (pred_test - Y_test), c ='r', alpha=0.5)

# グラフの整え
# 誤差0の線を引く。
plt.hlines(y = 0, xmin = -1.0, xmax = 50)
# 凡例
plt.legend((train, test), ('Training', 'Test'), loc = 'lower left') # 凡例
# タイトル
plt.title('Residual Plots')

誤差プロット結果

黒の横線が誤差=0を表していて、それより上に点があれば誤差プラス、下にあれば誤差マイナス。この誤差が上下に均一に分布していれば、線形回帰でこの予測を解くのに適していると言える。
しかし、この誤差が偏った結果となっていると、線形回帰を適応すること自体を見直す方が良い。

補足

【※1】 npmpyで使う用のデータ整形について

単回帰の直線の式は、y = ax + bと表現できるが、これはベクトル表記に直すと、
y = Apと書き換えることができる。

Aとpはそれぞれ、以下となり、この内積がyになる。(y = a * x + b * 1)

この形に整形をしている。

{.py}

# boston_df.RMだと、1次元配列。
X = boston_df.RM
# Xの行列の形を確認(506行だけの1次元配列)
X.shape
> (506,)
X
>0     6.575
>1     6.421
>2     7.185
>...
>504    6.794
>505    6.030
>Name: RM, Length: 506, dtype: float64

numpyでこれを利用するためには、二次元配列に変えてあげないといけない。

{.py}

# そこで二次元配列に変換するためにnp.vstackを使う。
X = np.vstack(boston_df.RM)
# Xの行列の形を確認(506行,1列の2次元配列)
X.shape
>(506, 1)
X
>array([[ 6.575],
>       [ 6.421],
>       [ 7.185],
>       ・・・
>       [ 6.794],
>       [ 6.03 ]])

{.py}

# さらにこれをA = [x, 1]に変換。Xの値を1次元目に、2次元目に固定値で1をセット
X = np.array([[value, 1] for value in X])
# Xの行列の形を確認(506行,2列の2次元配列)
X.shape
> (506, 2)
X
>array([[array([ 6.575]), 1],
>       [array([ 6.421]), 1],
>       [array([ 7.185]), 1],
>       ..., 
>       [array([ 6.976]), 1],
>       [array([ 6.794]), 1],
>       [array([ 6.03]), 1]], dtype=object)

これで説明変数Xが[x 1]の形となり、めでたくnumpyで扱える形になる。

【※2】 np.linalg.lstsqについて

linalgがlinear Algebra(線形代数)の略称で、lstsqがLeast Square(最小二乗法)の略称。

返り値は、0行目:傾きと切片、1行:残差、2行目:説明変数のRank、3行目:説明変数の特異値(Singlar value)が入る。今回は傾きと切片を求めたいので[0]を指定。