Pythonは使うけどライブラリは使わずに最小二乗法を３つの方法で解いてみた
機械学習の基礎中の基礎である単回帰分析をExcelで解く方法について、【6つもあった！】Excelで単回帰分析の最小二乗法を解く方法をすべて実演で６つの解法を紹介しました。
Pythonで単回帰分析をするときはscikit-learnを使うと一発でできますが、最小二乗法の理解のためにあえてライブラリを使わない３つの方法で解いてみました。
最後に比較対象として、scikit-learnを使うやり方も試してみました。
最小二乗法の公式
まずは単回帰分析を解くための公式を簡単に復習しておきます。
単回帰分析とは下記のようにxとyのデータがあるときに、未知のxからyを推定するための近似式を一次式で求めることです。
| 人数 x | 処理数 y |
| 17 | 173 |
| 4 | 60 |
| 6 | 80 |
| 10 | 127 |
| 5 | 63 |
| 9 | 111 |
一次式で求めるということは
$$y=ax+b$$
の式におけるaとbを決めるということですが、xとyのデータに中心化という処理を施すことによってaだけを求める問題に単純化することができます。
もう少し詳しく説明するために、先のデータxとyを散布図にしてみます。
赤線は管理人が適当に描いた近似直線ですが、このようにy切片（つまりb）が25近辺になることが予想されるため、bもちゃんと計算してやる必要があります。
ところが、これを中心化すると次のようになります。
原点を通る近似直線になりました。
中心化とは、xの各成分からxの平均を引いた変数をxcとして、yの各成分からyの平均を引いた変数をycとして、xcとycについて分析することですが、このようにすると近似式が原点を通る、つまりbが常にゼロとなるため、aだけを求めればよいということになります。
従って、
$$(y-yの平均)=a(x-xの平均)$$
におけるaを求める問題に単純化されます。
そしてこれを少し式変形すると、
\(y=a(x-xの平均)+yの平均\) ・・・ 式１
となります。
この式はxに未知の数が来たときにyを推定する式なので、推定値yを慣習的に\(\hat{y}\)と表します。
すると、この\(\hat{y}\)と真値であるyの差が最小になるようにaを求める問題に帰着します。
ここで微分を使って式変形をしていくと、
\(a=\frac{\sum_{i=1}^{N}x_iy_i}{\sum_{i=1}^{N}{x_i}^{2}}\) ・・・ 式２
ということが導かれます。
そしてこれを式１に代入すると、
\(b=(yの平均)-a(xの平均)\) ・・・ 式３
であることが導かれます。
これが単回帰分析を解く公式です。
最小二乗法の公式を使う２つの解法
それではこの公式にPythonを適用してみましょう。
次のように２通りのやり方があります。
そのままゴリゴリ計算する方法
式２、３をそのままPythonで計算する方法です。
式２についてはxベクトルの各要素の平均はx.mean()で計算でき、yベクトルの各要素の平均もy.mean()で計算できます。
また、xベクトルとyベクトルの第一成分同士の積＋第二成分同士の積＋、、、＋第ｎ成分同士の積はx*yという非常に簡単な式で計算できます。
従って、下記のようなコードで計算できます。
切片bは下記のように簡単に計算できます。
統計量からスマートに計算する方法
分散は偏差の二乗の平均です。
偏差とは平均値からのズレということで、10個のデータがあるとしたら、それぞれのデータが10個のデータの平均値からどれくらいズレているかを表します。
そしてそれら10個のズレをそれぞれ二乗してすべて足し合わせます。
最後にその合計値をデータ数である10で割った値が分散です。
従って、式３の分母はデータ数で割る前の分散です。
一方、共分散はXの偏差 × Y の偏差の平均です。
（特別な場合として、Xの偏差 × Xの偏差の平均は分散になります。従って、分散は共分散の特殊ケースともいえます）
式３の分子を見てみると、これはデータ数で割る前の共分散ということがわかります。
従って、式３の分母と分子をそれぞれデータ数で割ると、
$$a＝\frac{xとyの共分散}{xの分散}$$
であることがわかります。
このように考えると、aは共分散と分散という２つの統計量を計算することにより求められることになります。
共分散と分散はNumpyの関数で次のように簡単に計算できます。
気をつけないといけないのは、Numpyで計算される共分散は行列として計算されるということです。
２変数xとyの共分散行列は、
$$\begin{bmatrix} S_{xx} & S_{xy} \\ S_{yx} & S_{yy} \end{bmatrix}$$
となります。
ここで\(S_{xx}\)はxとxの共分散、つまりxの分散ということです。
確かにnp.var()で計算したxの分散と同じ値になっていますね。
また\(S_{xy}\)と\(S_{yx}\)は同じで、これが今知りたいxとyの共分散になります。
共分散行列から\(S_{xy}\)だけを取り出すには、行列の（１,２）成分、または（２,１）成分を抽出すればよいことになります。
１行目は０、そのうちの２列目は１を指定すればよいため、
で抽出できます。
これでxとyの共分散が求まったので、これを先に計算したxの分散で割ると、次のように\(a=9.04\)と計算されます。
切片bは先程と同じように計算できます。
行列を使う方法
これは【6つもあった！】Excelで単回帰分析の最小二乗法を解く方法をすべて実演の中で紹介した「行列を使う方法」と同じことをPythonを使って計算するだけです。
詳細は上の記事を参照していただくとして、次の行列方程式を解く問題に帰着します。
$$XA=Y$$
但し、
$$X=\begin{bmatrix} 17 &1\\ 4 & 1\\6 & 1\\10 & 1\\5 & 1\\9 &1 \end{bmatrix}$$
$$A=\begin{bmatrix} a \\ b \end{bmatrix}$$
$$Y=\begin{bmatrix} 173 \\ 60 \\ 80 \\ 127 \\ 63 \\ 111 \end{bmatrix}$$
通常、このような行列の方程式は\(X^{-1}\)を左から両辺に掛けることによって
$$A=X^{-1}Y$$
というようにして解くことができます。
ところが、このようにして解くことができるのは行列Xが正方行列で逆行列を持つことが前提です。
今回の場合は６行２列で正方行列ではありません。
どうしたらよいのでしょうか？
行と列をひっくり返した転置行列\(X^{T}\)を左からかけて正方行列にしてから逆行列を計算するという方法もあるのですが、Numpyには擬似逆行列として直接計算できる関数もあります。
疑似逆行列とは逆行列を一般化したもので、元の行列が正方行列でなくてもその擬似逆行列を左から掛けることによって単位行列にできるものです。
擬似逆行列は次のようにして計算できます。
このように計算した擬似逆行列\(X^{-1}\)を\(XA=Y\)の両辺に左から掛けることによって、\(A=X^{-1}Y\)によりAを計算することができます。
このように\(a=9.04\)、\(b=25.5\)と求まりました。
scikit-learnを使う方法
今までのやり方は公式を使うやり方で、公式を知らないとできません。
しかし、scikit-learnというライブラリーを使えば、一遍に求めることができてしまいます。
次のようにコーディングします。
前に計算した方法と同じ結果になりました。
実測データに対する単回帰式の当てはまり具合をグラフで見てみましょう。
当てはまり具合の良さも簡単にグラフで表示できました。