#はじめに
前回の記事では、回帰分析の１つである単回帰分析について書きました。
この記事では、変数同士の相関関係から変動を予測することができる重回帰分析について書きます。

#重回帰分析はどのような時に使われるのか。
重回帰分析は、マーケティングで使われる分析手法のひとつです。複数の変数を持つデータの関連性を明らかにします。
例えば、店舗面積/従業員数/商圏人口/広告費から、来客数を予測したり、
他にも、部屋の広さ/駅からの距離/犯罪発生率から、家賃を予測します。
このように、複数の要因から１つの値を予測するときに、重回帰分析を利用します。

また、重回帰分析では、それぞれの要因が予測したい値に与える影響度合を数値で表現できます
そのため、要因分析としても使われます。店舗面積/従業員数/商圏人口/広告費それぞれが、来店客数を増やす要因なのか、減らす要因なのかが分かります。
要因分析をすることで、来店客数を増やす有効なパラメータを定量的に評価することができ、マーケティングでも使われています。

重回帰分析は、様々な分野において、変動する値を予測するためや、要因分析のために使われています。
この記事では、重回帰分析の原理原則を理解することを目的とし、どのようなアルゴリズムで回帰係数を算出しているのかを紐解きます。

まずは、回帰分析の説明をします。

回帰分析とは

回帰分析とは、結果となる数値と要因となる数値の関係を数式(回帰式)で表す統計的手法です。
例えば、あるクラスのAさん〜Dさんの「勉強時間」と「テストの点数」のデータがあります。

このデータを視覚化すると、下のようなグラフになります。
横軸が勉強時間、縦軸がテストの点数を表しています。
勉強時間が長くなるほど、テストの点数が上がっている傾向がわかります。

4人のデータを用いて、「勉強時間」と「テストの点数」の関係を数式で表します。
「勉強時間」と「テストの点数」の関係が線形であると仮定して、データに当てはまる直線を下記のように引きます。

黒の直線は下の式(回帰式)で表せます。

y = 5.7x + 51

この式（回帰式）の意味は、


(点数)= 5.7 × (学習時間) + 51

となります。
学習時間が0のとき、51点となり、１時間勉強するごとに5.7点ずつ点数が上がっていることを意味しています。
この式（回帰式）を利用することで、勉強時間から、他の生徒がどの程度の点数を獲得できるかが予測できるようになります。
例えば、勉強時間が6時間のEさんの点数は、

\begin{align}
(Eさんの点数) &= 5.7× (学習時間) + 51 \\
&= 5.7 × 6 + 51 \\
&= 85(点)
\end{align}

と求めることができます。

このように、回帰分析は、結果となる数値と要因となる数値の関係を回帰式で表します。
また、要因となる数値を説明変数、結果となる数値を目的変数といいます。
今回の例でいうと、「勉強時間」から「点数」を予測したので、「勉強時間」が説明変数で、「点数」が目的変数となります。

回帰分析のうち、1つの目的変数を1つの説明変数で予測するものを単回帰分析と言います。
そして、1つの目的変数を２つ以上の説明変数で予測するものを重回帰分析と言います。

単回帰分析

y = a + b_1x_1 \\

重回帰分析

y = a + b_1x_1 + b_2x_2 + ・・・ + b_nx_n \\
y:目的変数 \\
a:切片 \\
b:回帰係数 \\
x:説明変数

では、どのようにしてインプットしたデータから、回帰式を算出しているのか。
それは、回帰式をコスト関数で算出したコストで評価し、最急降下法で、コストが小さい回帰係数を算出します。
コスト関数や、最急降下法を初めて聞く方は、前回の記事を一読した方がわかりやすいかもしれません。（参考：単回帰分析の原理原則）

本記事では、重回帰分析におけるコスト関数と最急降下法について述べていきます。
それでは、物理の点数、科学の点数のデータから数学の点数を予測するモデルケースを通して、重回帰分析を実装していきます。

#データの確認
CSVファイルを読み込み、データを確認します。

今回使用する学習データはこちらをクリック

`training_data_multipleLiner.csv`

Pysics	Science	Statistics	Math
64	67	69	68
74	74	63	67
60	59	55	57
84	88	89	91
80	88	78	82
75	68	68	63
66	59	55	65
77	71	66	73
70	71	60	67
89	84	80	85
73	76	69	82
78	74	63	67
67	63	67	57
80	81	83	74
72	75	73	83
68	63	76	71
81	89	80	85
64	66	59	69
79	81	79	87
70	76	69	71
76	80	68	80
71	64	75	66
81	77	84	73
85	83	87	90
61	67	70	61
74	71	75	73
85	82	78	84
60	60	57	52
83	88	80	92
85	80	76	87
84	86	86	80
82	79	88	83
68	66	60	71
89	85	88	91
85	90	85	98
68	66	62	61
74	76	82	77
76	75	84	81
62	68	70	60
64	67	68	73
65	66	59	58
64	57	56	57
73	74	63	69
65	65	69	69
87	91	97	97
62	60	61	57
77	69	62	70
77	76	75	70
89	97	94	99
63	56	59	62
74	69	62	75
72	67	67	67
69	77	81	73
88	86	95	86
77	84	75	87
83	80	77	80
77	77	80	74
84	83	86	87
77	79	69	83
66	74	71	74
74	69	79	75
67	75	78	73
83	85	85	93
72	65	73	69
70	76	73	83
71	70	74	75
89	94	82	86
82	86	85	89
86	80	88	75
85	89	100	97
89	83	79	76
77	74	74	71
79	83	93	89
75	71	79	78
68	66	68	67
64	56	57	58
66	60	55	57
81	76	80	80
68	72	71	79
71	70	75	72
62	54	56	54
80	85	93	85
80	84	86	83
66	62	67	58
73	77	69	74
70	70	71	67
62	55	62	56
62	70	76	76
74	67	65	65
80	83	78	78
90	84	83	90
84	76	78	70
84	87	96	92
78	81	80	78
90	92	81	89
60	68	55	65
86	85	80	93
72	64	59	69
69	63	62	57

import numpy as np
import pandas as pd

data = np.loadtxt('./training_data_multipleLiner.csv',delimiter=",",skiprows =1)
df = pd.DataFrame(data,columns=['物理','科学','統計','数学'])

	物理	科学	統計	数学
0	64.0	67.0	69.0	68.0
1	74.0	74.0	63.0	67.0
2	60.0	59.0	55.0	57.0
•	•	•	•	•
98	69.0	63.0	62.0	57.0

こちらのデータセットには、99人分の物理、科学、統計、数学のテストの点数の情報が入っています。
今回は、物理の点数と科学の点数から、数学の点数と予測する重回帰分析をしていきます。
つまり、説明変数が物理の点数と科学の点数として、目的変数を数学の点数とします。

初めに、説明変数と目的変数の関係を視覚化します。
x軸を物理の点数、y軸を科学の点数、z軸を数学の点数でデータをプロットします。

%matplotlib notebook
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
# 点数をプロットしてみる
fig = plt.figure()
ax = Axes3D(fig)
ax.scatter(df['物理'],df['科学'],df['数学'],color='#ef1234')
ax.set_xlabel('Pysics')
ax.set_ylabel('Science')
ax.set_zlabel('Math')
plt.show()

なんとなく、科学が高い人が数学も高そうな感じがしますが、正直よくわかりませんね。
このグラフからは、数学の点数を予測するのに、物理と科学、どちらの点数が重要なのか、どのくらい関係があるのかまではわかりません。

そこで重回帰分析により、下記の式の回帰係数と切片を算出し、物理の点数と科学の点数が、数学の点数に与える影響度合いを分析していきます。

(数学の点数) = 切片 + b_1×（物理の点数） + b_2（科学の点数） \\
b:回帰係数 \\

#データの前処理
元のデータをそのまま利用して、学習することも可能です。
（ここで言う学習は、データから回帰係数と切片を算出することを指しています。）
一般的に、各特徴量間でデータのスケールが異なる場合、事前にデータのスケールを調整したほうが、パラメータに偏りがなくなるので、より効率的に学習が行えるようになります。

今回の例で言うと、テストの点数なので最低点は０点、最高100点ですが、それぞれの教科の難易度により平均点や分散などが変わってきます。

そこでデータのスケールを合わせるだけでなく、各特徴量を標準正規分布に従うように変換する標準化を行います。

機械学習の手法は「入力のそれぞれの特徴量がどう変化したら、出力がどう変化するのか」を見ているので、各特徴量のばらつき度合い（標準偏差）を揃えておくことで、各特徴量の動きに対する感度を平等に見ることができるようになるわけです。

##データの標準化を行う。
標準化は下記に式で行います。
具体的には、各特徴量で、平均を０、標準偏差を１になるように調整します。

x = \frac{x - x_{mean}}{x_{std}} \\
x_{mean}:xの平均値 \\
x_{std}:xの標準偏差 \\

# 正規化(Z-score Normalization)
def norm(X):
    X_norm = np.zeros((X.shape[0],X.shape[1]))
    mean = np.zeros((1,X.shape[1]))
    std = np.zeros((1,X.shape[1]))
    for i in range(X.shape[1]):
        mean[:,i] = np.mean(X[:,i])
        std[:,i] = np.std(X[:,i])
        X_norm[:,i] = (X[:,i] - float(mean[:,i])) / float(std[:,i])
    return X_norm,mean,std

標準化を行う関数を定義したので、データの標準化を行います。
また、標準化ができているのか、平均と標準偏差を確認して見ましょう。

# 説明変数はX、目的変数はy
X=data[:,:2]
y=data[:,3]

# 標準化
X_norm,mean,std = norm(X)

# 標準化の確認
print("Xの平均値 = {}".format(X_norm.mean()))
print("Xの標準偏差 = {}".format(X_norm.std()))

ちゃんと標準化できていますね。

# Xの平均値 = 2.0634448134447354e-16
# Xの標準偏差 = 0.9999999999999999

##コスト関数
コスト関数(平均二乗誤差)は、ある回帰式が学習データとどれだけフィットしているかを定量的に評価します。
今回は、説明変数は科学の点数と物理の点数の２種類ですので、下記の式のようになります。

y = w_0 + w_1x_1 + w_2x_2 \\

\begin{align}
y&:目的変数 \\
w_0&:切片 \\
w_1,w_2&:回帰係数 \\
x&:説明変数 
\end{align}

例えば、(w₀,w₁,w₂) = (0,0,0)　のときは、下図のような平面になります。

コスト関数を数式で表すと下記の式になります。

\begin{align}
J(w_0,w_1,w_2) &= \frac{1}{2m}\sum_{i=1}^{m}(\hat{y}-y)^2 
\end{align}

また、回帰式をベクトルで表すと、

\begin{align}
\hat{y} &= w_0 + w_1x_1 + w_2x_2 \\
&= XW
\end{align}

X=\left[
\begin{matrix} x_0 & x_1 & x_2 \end{matrix}
\right]
,
W=\left[
\begin{matrix} w_0 \\ w_1 \\ w_2 \end{matrix}
\right]
(x_0 = 1)

と表せます。

数式で見ると難しそうに見えますが、やっていることは単純です。
まず、yは、学習データの目的変数（数学の点数）です。下の図で言うと、赤点のことです。
次に、y^は、学習データの説明変数（物理・科学の点数）に回帰係数をかけた予測値です。下の図で言うと、青い平面です。

コスト関数は、それぞれの点における予測値と、実際の答えの差（残差）を計算し、全ての残差を合計した値となります。（正確には、個数の２倍で割っています。）

ちなみに上図のコストを計算すると、2882となります。
次に、(w₀,w₁,w₂) = (50,3,3)　のときを考察して見ます。

先ほどよりは、若干データに近づいていることがわかります。
この時のコストを計算すると、552となります。
先ほどが2882だったのに対して、552と、コストが大幅に減ったことがわかります。
このように、コスト関数を用いることで、ある回帰式が学習データとどれだけフィットしているかを定量的に評価することができます。

それでは実装していきます。
まずは、切片の計算を工夫するために、先ほど標準化した説明変数に切片項（x₀）を追加します。
切片項（x₀）は１です。

また、(w₀,w₁,w₂)の初期値を(0,0,0)と定義します。

# 切片項を追加
X_padded = np.column_stack((np.ones((len(y),1)),X_norm))

# モデル式のパラメータ初期化
weight_int = np.zeros((3,1))

次に、コスト関数を定義します。

# コスト関数
def cost(X,y,weight):
    m = len(y)
    J = 0
    y_hut = X.dot(weight)
    diff = np.power((y_hut - np.transpose([y])),2)
    J = (1.0 /(2*m)) * diff.sum(axis = 0)
    return J

これでコスト関数を定義できました。
例えば先ほど(w₀,w₁,w₂) = (50,3,3)のときのコストは、下記のように求められます。

# コスト関数
print(cost(X_padded,y,[50,3,3])[0])
# 552

これで、ある回帰式が、学習データにどれほどフィットしているかを定量的に評価できるようになりました。
あとは、コストが最小となるような回帰係数を、最急降下法で求めていきます。

#最急降下法
最急降下法は、関数を最小にするようなパラメーターを求めるためのアルゴリズムです。
今回はコストが最小になる回帰係数（w₀,w₁,w₂)を求めます。

コスト関数は下記の式でした。

\begin{align}
J(w_0,w_1,w_2) &= \frac{1}{2m}\sum_{i=1}^{m}(\hat{y}-y)^2 \\
 &= \frac{1}{2m}\sum_{i=1}^{m}(w_0x_0 + w_1x_1 + w_2x_2 - y)^2
\end{align}

\begin{align}
m &: データ数  \\
\hat{y} &: 予測値 (\hat{y}=w_0x_0 + w_1x_1 + w_2x_2) \\
y &: 実測値
\end{align}

コスト関数をw₀で偏微分すると、

\begin{align}
\frac{\partial }{\partial w_0}J(w_0,w_1,w_2) &= \frac{\partial }{\partial w_0}\frac{1}{2m}\sum_{i=1}^{m}(w_0x_0 + w_1x_1 + w_2x_2 - y)^2 \\
&= \frac{1}{m}(w_0x_0 + w_1x_1 + w_2x_2 - y)
\end{align}

コスト関数をw₁で偏微分すると、

\begin{align}
\frac{\partial }{\partial w_1}J(w_0,w_1,w_2) &= \frac{\partial }{\partial w_1}\frac{1}{2m}\sum_{i=1}^{m}(w_0x_0 + w_1x_1 + w_2x_2 - y)^2 \\
&= \frac{1}{m}(w_0x_0 + w_1x_1 + w_2x_2 - y)x_1
\end{align}

コスト関数をw₂で偏微分すると、

\begin{align}
\frac{\partial }{\partial w_1}J(w_0,w_1,w_2) &= \frac{\partial }{\partial w_1}\frac{1}{2m}\sum_{i=1}^{m}(w_0x_0 + w_1x_1 + w_2x_2 - y)^2 \\
&= \frac{1}{m}(w_0x_0 + w_1x_1 + w_2x_2 - y)x_2
\end{align}

と求めることができます。
偏微分をすると、コスト関数の傾きがわかります。
コストが小さくなるように、w₀,w₁を更新したいので、元々の回帰係数から偏微分で得られた値を引きます。

\begin{align}
w_0 &:= w_0 - \alpha\frac{1}{m}(w_0x_0 + w_1x_1 + w_2x_2 - y) \\
w_1 &:= w_1 - \alpha\frac{1}{m}(w_0x_0 + w_1x_1 + w_2x_2 - y)x_1 \\
w_2 &:= w_1 - \alpha\frac{1}{m}(w_0x_0 + w_1x_1 + w_2x_2 - y)x_2
\end{align}

\begin{align}
\alpha &:学習率 \\
m &:データ数 \\
x &:説明変数 \\
y &:目的変数 \\
\end{align}

αは学習率と呼ばれ、回帰係数を更新する値の大きさを調整します。
学習率は、大きすぎると回帰係数が収束しない可能性があり、また、小さすぎると、回帰係数の収束に時間がかかってしまいます。
今回は学習率=0.01でと決めて、進めていきます。

では、最急降下法を実装していきます。

# 最急降下法
def gradientDescent(X,y,weight,alpha,iterations):
    m=len(y)
    J_history = np.zeros((iterations,1))
    weight_history = []
    
    for i in range(iterations):
        weight = weight - alpha * (1.0/m) * np.transpose(X).dot(X.dot(weight) - np.transpose([y]))
        weight_history.append(weight)
        J_history[i] = cost(X,y,weight)
        print(alpha * (1.0/m) * np.transpose(X).dot(X.dot(weight) - np.transpose([y])))
    return weight,J_history,weight_history

今回は、500回の回帰係数の更新する中で、コストがどのように推移するかを可視化するために、J_historyを追加しています。
また、更新された回帰係数を保持するために、weight_historyを追加しています。
次にα、学習回数の初期値を設定します。

# 学習率と学習回数(イテレーション)
alpha = 0.01
num_iters = 500

それではコスト関数と最急降下法を用いて、モデルの学習をしていきます。

weight,J_history,weight_history = gradientDescent(X_padded,y,weight_int,alpha,num_iters)
print(weight)
# [[74.57745626]
# [ 3.29981383]
# [ 7.16228608]]

回帰係数を求めることができました。
つまり、数学の点数は下記の式で表すことが出来ます。

(数学の点数)= 3.29×(標準化した物理の点数) +　7.16×(標準化した科学の点数)　＋ 74.5

どちらの回帰係数も正であるので、物理や科学の点数が高い方が数学の点数が高い傾向がわかります。
また、数学の点数には、物理の点数よりも科学の点数の方が、２倍以上関係があることを示しています。

また、回帰式と学習データを視覚化すると、

# 実測値の描画
fig = plt.figure()
ax = Axes3D(fig)
ax.scatter(data[:,0],data[:,1],data[:,3],color='red')

# 回帰平面の描画
X,Y = np.meshgrid(np.arange(data[:,0].min(),data[:,0].max(),1),np.arange(data[:,1].min(),data[:,1].max(),1))
X_norm = (X-float(mean[:,0])) / float(std[:,0])
Y_norm = (Y-float(mean[:,1])) / float(std[:,1])
Z = weight[0] + weight[1] * X_norm + weight[2] * Y_norm
ax.plot_surface(X,Y,Z,alpha=0.5)

# ラベルを追加
ax.set_xlabel('Pysics')
ax.set_ylabel('Science')
ax.set_zlabel('Math')
plt.show()

このように、学習データにフィットした平面を描くことができました。

また、コストの推移も確認してみましょう。

# コストと学習回数のグラフ
plt.plot(range(J_history.size),J_history,"-b",linewidth=0.5)
plt.xlabel("Number of Iterations")
plt.ylabel("Cost J")
plt.grid(True)
plt.show()

縦軸がコストで、横軸が学習回数を表しています。
学習を繰り返すごとに、コストが下がっているのがわかります。

##学習の推移を確認してみる。
500回、回帰係数を更新しました。その履歴がweight_historyに格納されているので、確認してみます。

nums = [0,99,199,299,399,499]

for num in nums:
    weight = weight_history[num]
    # 実測値の描画
    fig = plt.figure()
    ax = Axes3D(fig)
    ax.scatter(data[:,0],data[:,1],data[:,3],color='red')

    # 回帰平面の描画
    X,Y = np.meshgrid(np.arange(data[:,0].min(),data[:,0].max(),1),np.arange(data[:,1].min(),data[:,1].max(),1))
    X_norm = (X-float(mean[:,0])) / float(std[:,0])
    Y_norm = (Y-float(mean[:,1])) / float(std[:,1])
    Z = weight[0] + weight[1] * X_norm + weight[2] * Y_norm
    ax.plot_surface(X,Y,Z,alpha=0.5)

    # ラベルを追加
    ax.set_xlabel('Pysics')
    ax.set_ylabel('Science')
    ax.set_zlabel('Math')
    plt.savefig('{}回目.png'.format(num))

学習回数が１回、100回、200回、300回、400回、500回のときを描画します。

学習回数が１回のとき、(w₀,w₁,w₂) = (0.75,0.09,0.10)となり、下記のようになります。

まだ学習を1度しかしていないため、初期の平面から動きがほとんどありませんね。

学習回数が１00回のとき、(w₀,w₁,w₂) = (47.5,3.95,4.93)となり、下記のようになります。

100回の学習で、平面が大きく傾いてきました。
傾きはなんとなく良さそうですが、高さ（切片）が足りてないようです。

学習回数が200回のとき、(w₀,w₁,w₂) = (65.0,4.18,6.03)となり、下記のようになります。

パッとみた感じだと、これで良さそうに見えます。
学習回数が200回のときを見ると、コストは大幅に減少しています。

学習回数が300回のとき、(w₀,w₁,w₂) = (71.3,3.91,6.51)となり、下記のようになります。

4.18あった物理の回帰係数が、3.91に下がりました。
ここまで上昇傾向にあった係数が下がりました。微調整が行われています。

学習回数が400回のとき、(w₀,w₁,w₂) = (73.7,3.59,6.86)となり、下記のようになります。

学習回数が500回のとき、(w₀,w₁,w₂) = (74.5,3.29,7.16)となり、最終的に下記の通りに落ち着きます。

#まとめ
本記事では、重回帰分析についてまとめました。
重回帰分析の魅力は、関連性がわからないデータから、回帰係数という形で正・負の関連度合を算出できる点だと思います。
次は、重回帰分析を発展させたリッジ回帰やLasso回帰、Elastic Netにも取り組んでみたいと思ってます。