はじめに

多変量解析について理解を深めたいと考え、まず重回帰分析について勉強しました。
細かい式展開含めてほぼ完全に理解できたのではないか思ったので備忘録として投稿させていただきます。

参考

重回帰分析の理解に当たって下記を参考にさせていただきました。
どれも非常にわかりやすく、これから勉強される方にオススメです。

多変量解析法入門　著者; 永田靖, 棟近雅彦　出版社; サイエンス社
多変量解析のはなし　著者; 有馬哲, 石村貞夫　出版社; 東京図書
Chainer Tutorial-Chainerのチュートリアルサイト
【キカガク流】人工知能・機械学習脱ブラックボックス講座 - 中級編 - Udemyのオンライン講座

重回帰分析の理論

重回帰分析とは

重回帰分析を一言でいうと、ある結果（目的変数）を複数の原因（説明変数）から予測するモデルです。
数式で表すと下記のような感じ。

$$\hat{y}=w_0x_0+w_1x_1 + w_2x_2 + \cdots + w_nx_n$$

それぞれの説明変数を$x_1, x_2, x_3, \cdots, x_n$とする。
予測値を$\hat{y}$とする。
回帰係数（各説明変数の重み）を$w_1, w_2, w_3, \cdots, w_n$をとする。(これが求まれば勝ち)
$w_0$は切片。後々の計算を楽にする為、$x_0$(値は常に1)を掛け合わせた状態で記載する。

重回帰分析の最終的なゴールはできるだけ正確な予測モデルをつくることです。
つまり、重回帰分析が導き出した予測値と実測値の差が小さければ小さいほど良いのです。
これを数式に直すと、下記が最小になればOKということになります。

$$\sum_{n=1}^{n} (y_{n} -\hat{y}_{n} )^2$$

それではこのタスクを行列とベクトルを用いて解いていきます。

重回帰分析モデルの導出

まず各変数を下記のように行列・ベクトルの形式で表します。

\begin{split}\begin{aligned}
\boldsymbol{y}=
\begin{bmatrix}
y_{1} \\
y_{2} \\
\vdots \\
y_{n}
\end{bmatrix} \\
\end{aligned}\end{split}

\begin{split}\begin{aligned}
\boldsymbol{\hat{y}}=
\begin{bmatrix}
\hat{y}_{1} \\
\hat{y}_{2} \\
\vdots \\
\hat{y}_{n}
\end{bmatrix} \\
\end{aligned}\end{split}

\begin{split}\begin{aligned}
\boldsymbol{w}=
\begin{bmatrix}
w_{0} \\
w_{1} \\
\vdots \\
w_{n}
\end{bmatrix} \\
\end{aligned}\end{split}

\begin{split}\begin{aligned}
X=
\begin{bmatrix}
x_{10} & x_{11} & x_{12} & \cdots & x_{1m} \\
x_{20} & x_{21} & x_{22} & \cdots & x_{1m} \\
\vdots \\
x_{n0} & x_{n1} & x_{n2} & \cdots & x_{nm}
\end{bmatrix} \\
\end{aligned}\end{split}

$\boldsymbol{y}$は目的変数の実測値をベクトル化したもの
$\boldsymbol{\hat{y}}$は回帰式によって導き出された予測値をベクトル化したもの
$\boldsymbol{w}$は重回帰式を作成した際の回帰係数をベクトル化したもの
$X$はサンプル数$n$個、変数の数$m$個の説明変数の実測値を行列化したもの

この時、重回帰式によって導き出される予測値は下記のような数式で表されます。


\begin{split}\begin{aligned}
\boldsymbol{\hat{y}}=
\begin{bmatrix}
\hat{y}_{1} \\
\hat{y}_{2} \\
\vdots \\
\hat{y}_{n}
\end{bmatrix} \\
\end{aligned}\end{split}
=

\begin{split}\begin{aligned}
\begin{bmatrix}
x_{10} & x_{11} & x_{12} & \cdots & x_{1m} \\
x_{20} & x_{21} & x_{22} & \cdots & x_{1m} \\
\vdots \\
x_{n0} & x_{n1} & x_{n2} & \cdots & x_{nm}
\end{bmatrix} \\
\end{aligned}\end{split}

\begin{split}\begin{aligned}
\begin{bmatrix}
w_{0} \\
w_{1} \\
\vdots \\
w_{n}
\end{bmatrix} \\
\end{aligned}\end{split}

では、前述した解きたいタスクを行列とベクトルの形に変換していきます。

\sum_{n=1}^{n} (y_{n} -\hat{y}_{n} )^2\\
=(\boldsymbol{y}-\boldsymbol{\hat{y}})^T(\boldsymbol{y}-\boldsymbol{\hat{y}})\\
=(\boldsymbol{y}-X\boldsymbol{w})^T(\boldsymbol{y}-X\boldsymbol{w})\\
=(\boldsymbol{y}^T-(X\boldsymbol{w})^T)(\boldsymbol{y}-X\boldsymbol{w})\\
=(\boldsymbol{y}^T-\boldsymbol{w}^TX^T)(\boldsymbol{y}-X\boldsymbol{w})\\
=\boldsymbol{y}^T\boldsymbol{y}-\boldsymbol{w}^TX^T\boldsymbol{y}-\boldsymbol{y}^TX\boldsymbol{w}+\boldsymbol{w}^TX^TX\boldsymbol{w}\\
=\boldsymbol{y}^T\boldsymbol{y}-2\boldsymbol{y}^TX\boldsymbol{w}+\boldsymbol{w}^TX^TX\boldsymbol{w}

すっきりとした形になりました。
ちなみに上記の式展開で使用した公式はこちらです。

\begin{split}\begin{aligned}
&\left( 1\right) \ \left( {\bf A}^{\rm T} \right)^{\rm T} = {\bf A} \\
&\left( 2\right) \ \left( {\bf A}{\bf B} \right)^{\rm T} = {\bf B}^{\rm T}{\bf A}^{\rm T}\\
&\left( 3\right) \ \left( {\bf A}{\bf B}{\bf C} \right)^{\rm T} = {\bf C}^{\rm T}{\bf B}^{\rm T}{\bf A}^{\rm T}
\end{aligned}\end{split}

上記の式で求めたいのは、重み$\boldsymbol{w}$の値です。
なので、それ以外は定数項としてまとめてしまいます。

\boldsymbol{y}^T\boldsymbol{y}-2\boldsymbol{y}^TX\boldsymbol{w}+\boldsymbol{w}^TX^TX\boldsymbol{w}\\
=c+b\boldsymbol{w}+\boldsymbol{w}^TA\boldsymbol{w}\\

定数項は下記のようにまとめています。

A=X^TX\\
b=-2\boldsymbol{y}^TX\\
c=\boldsymbol{y}^T\boldsymbol{y}

先ほど導き出した式を下記のように置きます。
$$L=c+b\boldsymbol{w}+\boldsymbol{w}^TA\boldsymbol{w}$$
解きたいタスクは、$L$が最小となるような重み$\boldsymbol{w}$を求めるというものです。
これは$L$を$\boldsymbol{w}$で微分した式を0と置き、それを解くことで求めることができます。
なぜ微分を使って導出できる理由は、簡単に言うと目的関数（ここでいう$L$）が$\boldsymbol{w}$の各値について二次関数になっているからです。そのため、微分して0と置くことで最小値を求めることができるのです。より詳しい解説が知りたい方は「Chainer Tutorial」をご確認ください。

それでは$L$を$\boldsymbol{w}$で微分していきます。


\begin{split}\begin{aligned}
\frac{\partial}{\partial {\boldsymbol{w}}} L
&= \frac{\partial}{\boldsymbol{w}} (c + b\boldsymbol{w} + \boldsymbol{w}^T{A}\boldsymbol{w}) \\
&=\frac{\partial}{\partial {\boldsymbol{w}}} (c) + \frac{\partial}{\partial {\boldsymbol{w}}} ({b}{\boldsymbol{w}}) + \frac{\partial}{\partial {\boldsymbol{w}}} ({\boldsymbol{w}}^{T}{A}{\boldsymbol{w}}) \\
&={0} + {b} + {w}^{T}({A} + {A}^{T})
\end{aligned}\end{split}

ちなみに上記の微分で使用した公式はこちらです。
(1)(2)の導出方法は「Chainer Tutorial」
(3)の導出方法は「ベクトルと行列による微分」が詳しいです。

\begin{split}\begin{aligned}
&\left( 1 \right) \ \frac{\partial}{\partial {\bf x}} \left( c \right) = {\bf 0} \\
&\left( 2 \right) \ \frac{\partial}{\partial {\bf x}} \left( {\bf b}{\bf x} \right) = {\bf b} \\
&\left( 3 \right) \ \frac{\partial}{\partial {\bf x}} \left( {\bf x}^{\rm T}{\bf A}{\bf x} \right) = {\bf x}^{\rm T} \left( {\bf A} + {\bf A}^{\rm T} \right)
\end{aligned}\end{split}

それでは、微分した結果を0と置いて式を展開していきます。


\begin{split}\begin{aligned}
-2{\boldsymbol{y}}^{T}{X} + {\boldsymbol{w}}^{T}({X}^{T}{X} + ({X}^{T}{X})^T)
&= {0} \\
-2{\boldsymbol{y}}^{T}{X} + 2{\boldsymbol{w}}^{T}{X}^{T}{X}
&= {0} \\
{\boldsymbol{w}}^{T}{X}^{T}{X}
&= {\boldsymbol{y}}^{T}{X}  \\
({\boldsymbol{w}}^{T}{X}^{T}{X})^T 
&= ({\boldsymbol{y}}^{T}{X})^T  \\
{X}^{T}{X}{\boldsymbol{w}}
&= X^T{\boldsymbol{y}}  \\
({X}^{T}{X})^{-1}{X}^{T}{X}{\boldsymbol{w}}
&=({X}^{T}{X})^{-1}X^T{\boldsymbol{y}}  \\
{\boldsymbol{w}}
&=({X}^{T}{X})^{-1}X^T{\boldsymbol{y}}  \\

\end{aligned}\end{split}

これで求めたかった重み$\boldsymbol{w}$を導出できました。
ただし下記注意点があります。

上記数式ですが、最後両辺に逆行列$({X}^{T}{X})^{-1}$を掛けています。従って、${X}^{T}{X}$の逆行列が計算可能な時のみ$\boldsymbol{w}$の解が定まります。そのためには行列が正則である必要があり、そうならない場合は大きく分けて下記2つの理由が考えられます。

①説明変数の数がサンプル数よりも大きい
→重回帰分析は、要するに連立方程式を解いてるので未知数の数が式の数よりも多いと解が求まりません。
②多重共線性が存在する（説明変数同士の相関が強く、ある説明変数が他の説明変数の線形結合で表される状態）
→このとき、正則の条件である行列のランク数とサンプル数の一致($rankA=n$)が実現されません。

次はこちらの理論をPythonで自力実装していきます。

重回帰分析を理解する（理論編）

はじめに

参考

重回帰分析の理論

重回帰分析とは

重回帰分析モデルの導出

Next