行列を使って回帰分析してみる - 統計を学ぶ化学系技術者の記録

統計や機械学習について調べると、行列を使った解説が多数出てくる。何かの役に立つことを期待して行列を使った計算方法も調べておく。ちなみに線形代数がすこぶる苦手なので、計算はpythonに頼りっぱなしである。

データの準備

データはいつものように下記のものを使用する。

$x$ (説明変数)	1	2	3	4	5
$y$ （被説明変数）	2	6	6	9	6

エクセルで得られる回帰直線がこちら。これと同じ結果を、行列の計算から出すことが出来れば成功だ。

f:id:Chemstat:20200804164539j:plain

データを行列にしてみる

　データの集団が回帰直線 $y=\hat{b}x+\hat{a}$ で近似されると、個々のデータは残差項 $\varepsilon$ を加えて $y_{i}=\hat{b}x_{i}+\hat{a}+\varepsilon_{i}$ で表される。

f:id:Chemstat:20200806191422j:plain

今回だと5つのデータがあり、それぞれ下記左のような5つの式になる。これを行列としてまとめると下記右のような式が出来上がる。行列の計算をすると、ちゃんと対応していることが分かる。

定数項 $\hat{a}$ を表すために $\boldsymbol{X}$ の2列目がすべて1になっているのだが、pythonのstatsmodelsを使う際に、定数項を足す操作がこれに相当する。（たぶん。勉強してて気づいただけ。）

f:id:Chemstat:20201118062800j:plain

最小二乗法とは残差の平方和 $\small\begin{align*}\sum_{i=1}^n\varepsilon_i^2\end{align*}$ を最小にする $\hat{b}$ と $\hat{a}$ を求める作業である。これも行列で書くとこんな感じになる。

$\boldsymbol{\varepsilon}^{T}$ というのは転置行列で、元の行列のXYを入れ替えたものだ。

f:id:Chemstat:20201118071448j:plain

この辺の行列の基礎知識がなかったので、最初に数式を見た時はずいぶん警戒してしまったのだが、中身をみたら非常に単純だった。大学で線形代数の授業を受けてたはずなのだが、そのころの記憶はどこに行ってしまったのだろうか。

残差平方和が最小になる解を求める

残差平方和を最小にするということはつまり $\boldsymbol{\|\varepsilon\|^2=\|Y-Xb\|^2}$ を最小にする $\boldsymbol{b}$ を求めることに等しい。

この解を求めるには正規方程式という便利なものがあり、 $\boldsymbol{X^TXb=X^TY}$ を解くことで得られる。ということで早速解いていく。

f:id:Chemstat:20201118091039j:plain

このようにして無事 $\hat{b}=1.1$ 、 $\hat{a}=2.5$ が得られた。ちなみに実際の計算はpythonでやってるので気になる方は下記の参考コードをご覧ください。多分自力でもそんなに難しくないけど。

参考：正規方程式の導出

上のほうでいきなり正規方程式が出てきたが、その導出をここに記載する。ちなみに参考サイトに丁寧に解説されているので、詳しくはそちらを見て頂ければと思う。

f:id:Chemstat:20201118084731p:plain

最後の式変形で $\boldsymbol{{Y}^{T}Xb}=\boldsymbol{{b}^{T}{X}^{T}Y}$ をつかっている。これはこの値が１×１の行列になるので、転置しても同じという事らしい。一応計算してみた。

f:id:Chemstat:20201118084835j:plain
ちゃんときれいに定数になっている。

$\boldsymbol{X}$ と $\boldsymbol{Y}$ はデータとして与えられている定数なので、残差平方和 $\small\begin{align*}\sum_{i=1}^n\varepsilon_i^2\end{align*}$ は $\boldsymbol{b}$ の関数になっている。 $\small\begin{align*}\sum_{i=1}^n\varepsilon_i^2\end{align*}$ が最小になる $\boldsymbol{b}$ を探すには、偏微分して極値を求めるとそこが最小値となる（らしい）。

f:id:Chemstat:20201118085820p:plain

偏微分するとこのような式になり、これがゼロになるとき、つまり

f:id:Chemstat:20201118085802p:plain

を満たす $\boldsymbol{b}$ を求めればよいことになる。

ちなみに行列の偏微分の公式（偏微分されたら転置されるとか、対象行列の偏微分とか）がないと混乱するのでこちらも参考サイトをご覧ください。ちなみに私は調べて初めて知りました。

参考：pythonコード

import numpy as np

#Xのデータ
x_data = np.matrix([[1, 2, 3, 4, 5]]).T

#切片項
const = np.matrix([[1, 1, 1, 1, 1]]).T

#Yデータ
y_data = np.matrix([[2, 6,6, 9, 6]]).T

#Xと切片の行列
x_mat = np.hstack([x_data,const])

#傾きbと切片aの計算
print(x_mat.T @ x_mat).I @ (x_mat.T @ y_data)