回帰分析の理解まとめ 〜 重回帰分析 〜

前回こちらの記事で回帰分析のうち単回帰分析について概要とパラメータの求め方をまとめした。
takarotoooooo.hatenablog.com

今回は重回帰分析に関してまとめてみようと思います。

重回帰分析とは？

ある目的変数 y を複数の説明変数 x_1, x_2, ... , x_M で表す式を求めることです。
例えば、目的変数を下記の式で近似できると仮定した時に、最適なパラメータ w_1, w_2, ... , w_M を求めます。

y = w_{1}x_{1} + w_{2}x_{2} + ... + w_{M}x_{M}
ここで、 w_{0}x_{0}をバイアス項（全ての xが0の時の目的変数の値）とし、
\boldsymbol{x} = \begin{pmatrix} x_0 \\ x_1 \\ \vdots \\ x_M \end{pmatrix}, \boldsymbol{w} = \begin{pmatrix} w_0 \\ w_1 \\ \vdots \\ w_M \end{pmatrix}
とすると
\begin{eqnarray} y = \boldsymbol{w}^T\boldsymbol{x} \end{eqnarray}
と表すことができます。

評価関数

重回帰分析においても、評価関数として二乗誤差を使います。
y_i, \hat{y_i}をそれぞれ、 i番目の実測値、予測値とすると
\begin{eqnarray} \mathcal{L}(\boldsymbol{w}) &=& \sum_{i=1}^{N}{ ( y_i - \hat{y_i} )^2 } \\ &=& \begin{pmatrix} {{y_1} - \hat{y_1}} & {{y_2} - \hat{y_2}} & \cdots & {{y_N} - \hat{y_N}} \end{pmatrix} \begin{pmatrix} {{y_1} - \hat{y_1}} \\ {{y_2} - \hat{y_2}} \\ \vdots \\ {{y_N} - \hat{y_N}} \end{pmatrix} \end{eqnarray}

\boldsymbol{y} = \begin{pmatrix} y_1 \\ y_2 \\ \vdots \\ y_N \end{pmatrix}, \hat{ \boldsymbol{y} } = \begin{pmatrix} \hat{y_1} \\ \hat{y_2} \\ \vdots \\ \hat{y_N} \end{pmatrix}
とすると

\begin{eqnarray} \mathcal{L}(\boldsymbol{w}) &=& (\boldsymbol{y} - \hat{\boldsymbol{y}})^T (\boldsymbol{y} - \hat{\boldsymbol{y}}) \end{eqnarray}

評価関数が最小値をとるパラメータを探す

まずは先ほどおいた \hat{\boldsymbol{y}} を下記のように展開します
\hat{\boldsymbol{y}} = \begin{pmatrix} { \hat{y_1} } \\ { \hat{y_2} } \\ \vdots \\ { \hat{y_N} } \end{pmatrix} = \begin{pmatrix} { \boldsymbol{x_1}^T{w_1} } \\ { \boldsymbol{x_2}^T{w_2} } \\ \vdots \\ { \boldsymbol{x_N}^T{w_N} } \end{pmatrix} = \begin{pmatrix} { \boldsymbol{x_1}^T } \\ { \boldsymbol{x_2}^T } \\ \vdots \\ { \boldsymbol{x_N}^T } \end{pmatrix} \boldsymbol{w} = \begin{pmatrix} { x_{10} } & { x_{11} } & { x_{12} } & \cdots & { x_{1M} } \\ { x_{20} } & { x_{21} } & { x_{22} } & \cdots & { x_{2M} } \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ { x_{N0} } & { x_{N1} } & { x _{N2} } & \cdots & { x_{NM} } \end{pmatrix} \boldsymbol{w}

ここで

\boldsymbol{X} = \begin{pmatrix} { x_{10} } & { x_{11} } & { x_{12} } & \cdots & { x_{1M} } \\ { x_{20} } & { x_{21} } & { x_{22} } & \cdots & { x_{2M} } \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ { x_{N0} } & { x_{N1} } & { x _{N2} } & \cdots & { x_{NM} } \end{pmatrix}

とすると、下記のように表せます。

\hat{\boldsymbol{y}} = \boldsymbol{X}\boldsymbol{w}

これを評価関数に代入し、式を展開すると下記のように表せます。

\begin{eqnarray} \mathcal{L}(\boldsymbol{w}) &=& (\boldsymbol{y} - \hat{\boldsymbol{y}})^T (\boldsymbol{y} - \hat{\boldsymbol{y}}) \\ &=& (\boldsymbol{y} - \boldsymbol{X}\boldsymbol{w})^T (\boldsymbol{y} - \boldsymbol{X}\boldsymbol{w}) \\ &=& (\boldsymbol{y}^T - \boldsymbol{w}^T\boldsymbol{X}^T) (\boldsymbol{y} - \boldsymbol{X}\boldsymbol{w}) \\ &=& \boldsymbol{y}^T \boldsymbol{y} - \boldsymbol{w}^T\boldsymbol{X}^T \boldsymbol{y} - \boldsymbol{y}^T \boldsymbol{X}\boldsymbol{w} + \boldsymbol{w}^T\boldsymbol{X}^T \boldsymbol{X}\boldsymbol{w} \\ &=& \boldsymbol{y}^T \boldsymbol{y} - 2\boldsymbol{y}^T \boldsymbol{X}\boldsymbol{w} + \boldsymbol{w}^T\boldsymbol{X}^T \boldsymbol{X}\boldsymbol{w} \end{eqnarray}

ここで評価関数を最小にする \boldsymbol{w}を求めます。

\begin{eqnarray} \frac{\partial}{\partial{\boldsymbol{w}}}\mathcal{L}(\boldsymbol{w}) &=& \frac{\partial}{\partial{\boldsymbol{w}}}{({\boldsymbol{y}^T}{\boldsymbol{y}} - 2{\boldsymbol{y}^T}{\boldsymbol{X}}{\boldsymbol{w}} + {\boldsymbol{w}^T\boldsymbol{X}^T}{\boldsymbol{X}\boldsymbol{w}})} \\ &=& \frac{\partial}{\partial{\boldsymbol{w}}}{({\boldsymbol{y}^T}{\boldsymbol{y}})} - \frac{\partial}{\partial{\boldsymbol{w}}}{(2{\boldsymbol{y}^T}{\boldsymbol{X}}{\boldsymbol{w}})} + \frac{\partial}{\partial{\boldsymbol{w}}}{({\boldsymbol{w}^T\boldsymbol{X}^T}{\boldsymbol{X}\boldsymbol{w}})} \\ &=& \boldsymbol{0} - 2{\boldsymbol{X}^T}{\boldsymbol{y}} + ({\boldsymbol{X}^T}{\boldsymbol{X}}+({\boldsymbol{X}^T}{\boldsymbol{X}})^T)\boldsymbol{w} \\ &=& - 2{\boldsymbol{X}^T}{\boldsymbol{y}} + 2{\boldsymbol{X}^T}{\boldsymbol{X}}\boldsymbol{w} \end{eqnarray}

つまり

\begin{eqnarray} \frac{\partial}{\partial{\boldsymbol{w}}}\mathcal{L}(\boldsymbol{w}) &=& 0 \\ {- 2{ \boldsymbol{X}^T }{ \boldsymbol{y} } } + 2{ \boldsymbol{X}^T }{ \boldsymbol{X} }\boldsymbol{w} &=& 0 \\ { \boldsymbol{X}^T }{ \boldsymbol{X} }\boldsymbol{w} &=& { \boldsymbol{X}^T }{ \boldsymbol{y} } \\ \end{eqnarray}

つまり、下記が成り立ちます

\begin{eqnarray} \boldsymbol{w} &=& ({ \boldsymbol{X}^T }{ \boldsymbol{X} })^{-1}{ \boldsymbol{X}^T }{ \boldsymbol{y} } \end{eqnarray}

実際に求めてみる

パラメータの求め方がわかったので、実際に以下のサンプルデータを使って求めてみます。

x_1 x_2 x_3 x_4 y
0 0 0 0 -5.092
1 1 9 1 55.14354986
・・・ ・・・ ・・・ ・・・ ・・・
48 3 2 0 89.69
49 4 1 1 101.47423445

このデータを横軸を x_1、縦軸を x_2としてプロットしたものが下記です。

このデータをもとに \boldsymbol{w}を下記の式から求めると
\begin{eqnarray} ({ \boldsymbol{X}^T }{ \boldsymbol{X} })^{-1}{ \boldsymbol{X}^T }{ \boldsymbol{y} } \end{eqnarray}

w_0 w_1 w_2 w_3 w_4
-0.78 2.069 -4.28 5.41 11.69

と、なり
\hat{\boldsymbol{y}} = \boldsymbol{X}\boldsymbol{w}
を計算して予測データとしてプロットしたものが下記となり、おおよそ観測データを近似できているように見えます。


\boldsymbol{w}が示す意味

\boldsymbol{w}のそれぞれの値 w_1, ... , w_M は単回帰分析の時と同様に x_1, x_2, ... , x_M が y に与える影響度の大きさを表しています。
ただし、 x_1, x_2, ... , x_M それぞれの値のスケールが違うことで、 w_1 < w_2 の時に x_1 よりも x_2 の方が y により大きな影響を与えるということが言えません。

パラメータ \boldsymbol{w}の値を見るだけで、どの変数 x_i が y に対して、どの程度の影響があるのかを知れるようにするために、各 x_i の値を標準化した上でパラメータ \boldsymbol{w}を求めます。

標準化

各変数の標準化の方法は変数が従う確率分布によって異なりますが、ここでは各変数が正規分布に従うと仮定して標準化してみます。
各 x_i に関して平均を引き、標準偏差で割ることで平均0、分散1の値として標準化することができます。

x_1 x_2 x_3 x_4
-1.69 -1.41 -1.56 -1
-1.62 -0.70 1.56 1
・・・ ・・・ ・・・ ・・・
1.62 0.70 -0.87 -1
1.69 1.41 -1.21 1

この標準化された \boldsymbol{X}_{scaling} をもとに算出されたパラメータ \boldsymbol{w}_{scaling} が下記です。

w_0 w_1 w_2 w_3 w_4
71.54 29.85 -6.05 15.54 5.84

これを見ると w_1 の絶対値が大きいので、 x_1 が y への影響度が一番大きいことが判断できます。

まとめ

  • 重回帰分析は、ある目的変数 y を複数の説明変数 x_1, x_2, ... , x_M で表す式を求めること
  • 最適なパラメータ w は行列を用いて ({ \boldsymbol{X}^T }{ \boldsymbol{X} })^{-1}{ \boldsymbol{X}^T }{ \boldsymbol{y} } で求められる
  • パラメータ w は説明変数を標準化することで、説明変数毎の目的変数への影響度を表すことができる。

今回は重回帰分析に関して概要からパラメータの求め方までをまとめてみました。
次回はロジスティック回帰分析に関してまとめてみようと思います。