Pythonを使って単回帰分析における最小二乗法を4つの方法で解いてみた

Photo by charlesdeluvio on Unsplash

機械学習の基礎中の基礎である単回帰分析をExcelで解く方法について、【6つもあった!】Excelで単回帰分析の最小二乗法を解く方法をすべて実演で6つの解法を紹介しました。

今回はPythonで解く4つの解法を紹介したいと思います。

まだPythonを勉強し始めたばかりなので、管理人が気づかないだけで他にもあると思いますが、現時点での知識ということでまとめました。

 

【沸騰アジア!】

親日国カンボジアであなたの会社の商品を販売してみませんか?

食品/飲料ヘルスケアコスメ商流と物流を管理人がサポートします。

Naga World in Phnom Penh (2022)

単回帰分析の公式

まずは単回帰分析を解くための公式を簡単に復習しておきます。

単回帰分析とは下記のように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の平均\) ・・・ 式1

となります。

この式は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}}\) ・・・ 式2

ということが導かれます。

そしてこれを式1に代入すると、

\(b=(yの平均)-a(xの平均)\)   ・・・ 式3

であることが導かれます。

これが単回帰分析を解く公式です。

 

詳しくはこちらのUdemyの講座でとてもわかりやすく解説されています。

【キカガク流】人工知能・機械学習 脱ブラックボックス講座 – 初級編 –

1000人以上が受講している(株)キカガクの『脱ブラックボックスセミナー』が遂に登場!機械学習の参考書を「閉じてしまった人」への再入門に最適な講座です。微分・線形代数といった数学の基礎からPythonでの実装まで短時間で習得しましょう。

 

公式を使う2つの解法

それではこの公式にPythonを適用してみましょう。

次のように2通りのやり方があります。

 

そのまま計算する方法

式2、3をそのままPythonで計算する方法です。

式2についてはxベクトルの各要素の平均はx.mean()で計算でき、yベクトルの各要素の平均もy.mean()で計算できます。

また、xベクトルとyベクトルの第一成分同士の積+第二成分同士の積+、、、+第n成分同士の積x*yという非常に簡単な式で計算できます。

従って、下記のようなコードで計算できます。

切片bは下記のように簡単に計算できます。

 

統計量として計算する方法

分散偏差の二乗の平均です。

偏差とは平均値からのズレということで、10個のデータがあるとしたら、それぞれのデータが10個のデータの平均値からどれくらいズレているかを表します。

そしてそれら10個のズレをそれぞれ二乗してすべて足し合わせます。

最後にその合計値をデータ数である10で割った値が分散です。

従って、式3の分母はデータ数で割る前の分散です。

 

一方、共分散Xの偏差 × Y の偏差の平均です。

(特別な場合として、Xの偏差 × Xの偏差の平均は分散になります。従って、分散は共分散の特殊ケースともいえます)

式3の分子を見てみると、これはデータ数で割る前の共分散ということがわかります。

 

従って、式3の分母と分子をそれぞれデータ数で割ると、

$$a=\frac{xとyの共分散}{xの分散}$$

であることがわかります。

このように考えると、aは共分散分散という2つの統計量を計算することにより求められることになります。

 

共分散と分散はNumpyの関数で次のように簡単に計算できます。

 

気をつけないといけないのは、Numpyで計算される共分散は行列として計算されるということです。

2変数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}\)だけを取り出すには、行列の(1,2)成分、または(2,1)成分を抽出すればよいことになります。

1行目は0、そのうちの2列目は1を指定すればよいため、

で抽出できます。

これで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が正方行列逆行列を持つことが前提です。

今回の場合は6行2列で正方行列ではありません。

どうしたらよいのでしょうか?

 

行と列をひっくり返した転置行列\(X^{T}\)を左からかけて正方行列にしてから逆行列を計算するという方法もあるのですが、Numpyには擬似逆行列として直接計算できる関数もあります。

疑似逆行列とは逆行列を一般化したもので、元の行列が正方行列でなくてもその擬似逆行列を左から掛けることによって単位行列にできるものです。

擬似逆行列は次のようにして計算できます。

 

このように計算した擬似逆行列\(X^{-1}\)を\(XA=Y\)の両辺に左から掛けることによって、\(A=X^{-1}Y\)によりAを計算することができます。

このように\(a=9.04\)、\(b=25.5\)と求まりました。

 

scikit-learnを使う解法

今までのやり方は公式を使うやり方で、公式を知らないとできません。

しかし、scikit-learnというライブラリーを使えば、一遍に求めることができてしまいます。

次のようにコーディングします。

 

前に計算した方法と同じ結果になりました。

実測データに対する単回帰式の当てはまり具合をグラフで見てみましょう。

 

当てはまり具合の良さも簡単にグラフで表示できました。

 

まとめ

実際にPythonで単回帰分析をするときはscikit-learnを使うのだと思いますが、理論の理解のために4つの方法で解いてみました。

理論の部分がわかりにくかった人は、Udemyの【キカガク流】人工知能・機械学習 脱ブラックボックス講座 – 初級編 –の講座でやさしく解説されていますので受講してみてはいかがでしょうか。

 

Udemy【キカガク流】人工知能・機械学習 – 初級編 – を受講してみた感想