線形回帰

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

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

<img width="326" alt="スクリーンショット 2016-05-03 11.42.20.png" src="https://qiita-image-store.s3.amazonaws.com/0/40159/c17412a1-c451-73d7-859d-b5e097be4663.png">

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

上図でいうと、赤点から青線へ向かって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()
```

<img width="232" alt="スクリーンショット 2016-05-03 12.07.01.png" src="https://qiita-image-store.s3.amazonaws.com/0/40159/c4bbf590-2914-c953-4e29-184b94ad2ac9.png">

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')
```

<img width="403" alt="スクリーンショット 2016-05-03 12.10.07.png" src="https://qiita-image-store.s3.amazonaws.com/0/40159/7ef950b4-fcce-1625-b747-3ec9f84ea19e.png">

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

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

```{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) 
```

<img width="367" alt="スクリーンショット 2016-05-03 12.22.38.png" src="https://qiita-image-store.s3.amazonaws.com/0/40159/5c0b15a2-44db-182f-9379-b16b101aedb2.png">

めっちゃ簡単。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での線形回帰の結果の図。赤線が得られた回帰式。

<img width="403" alt="スクリーンショット 2016-05-03 12.54.10.png" src="https://qiita-image-store.s3.amazonaws.com/0/40159/f77b1479-fa78-f1d1-1c8c-f4b69fb0dc14.png">

## やってみる（重回帰）
### 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')
```

誤差プロット結果

<img width="395" alt="スクリーンショット 2016-05-03 14.41.15.png" src="https://qiita-image-store.s3.amazonaws.com/0/40159/67362d61-33be-b104-3a80-2dd5f39627c3.png">

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



### 補足
#### 【※1】  npmpyで使う用のデータ整形について
単回帰の直線の式は、y = ax + bと表現できるが、これはベクトル表記に直すと、
y = Apと書き換えることができる。

Aとpはそれぞれ、以下となり、この内積がyになる。(y = a * x + b * 1)
<img width="97" alt="スクリーンショット 2016-05-03 13.54.20.png" src="https://qiita-image-store.s3.amazonaws.com/0/40159/d15bace4-94b2-b03b-ac90-d4bec7a14ad5.png">

<img width="77" alt="スクリーンショット 2016-05-03 13.54.28.png" src="https://qiita-image-store.s3.amazonaws.com/0/40159/7a9bf2c8-01b6-9d77-2e76-ff210224f203.png">

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

```{.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]を指定。