【機械学習入門③】最小二乗法とは？手計算とPythonを使った実装を解説【初心者向け】

前回の記事『【機械学習入門②】線形回帰・単回帰モデルをPythonの実装付きで理解！
』では、線形回帰の単回帰モデルについて触れてきました。

今回は、単回帰モデルを具体的に導出する方法である「最小二乗法」について触れていきたいと思います。

この記事で学べる内容は、以下のとおりです。

最小二乗法は、式展開も含めて勉強になるので、ゆっくり読み進めていただけたらと思います(*’ω’*)

最小二乗法とは？【分かりやすく解説】

プログラムを使うと、単回帰モデルは一瞬で求めることができます。

ただ、プログラムを書けるだけだとぶっちゃけ微妙なので、その裏で使われている「最小二乗法」について学んでいきましょう。

1.1. 【前回のおさらい】線形回帰モデルの求め方

以前に学習した単回帰モデルを求めるには、「各点からなるべく近い位置を通る直線を求める」という解説をしていました。

図で確認するなら、以下のとおりです。

上の図のように、すべての点を通る直線は引けないので、各点からなるべく近い位置を通る直線を引いています。

そして、この各点からなるべく近い位置を通る直線を求めるのに必要なアクションは、\(y = w_0 + w_1 x_1\)の\(w_0\)と\(w_1\)を求めることでした。

ここで疑問になってくるのが、「どうやって、\(y = w_0 + w_1 x_1\)の\(w_0\)と\(w_1\)を求めるんだろう」です。

プログラムで書いてしまえば、上の図のように勝手に出力してくれますが、、、とはいえ、どうやって計算されているのか心残りになります。

そこで登場するのが、今回学習する最小二乗法です。

最小二乗法を使えば、「各点からなるべく近い位置を通る直線(=すなわち単回帰モデル)」のパラメータ(\(w_0\)と\(w_1\))を求められます。

1.2. 最小二乗法の式

実際のデータを\(y_i\)、予測値を\(\hat{y_i}\)とすると、最小二乗法の式は以下のようになります。

$$
\begin{eqnarray}
E(w, b) &=& \frac{1}{2} \sum_{i=1}^{N} (y_i – \hat{y_i})^2\\
&=& \frac{1}{2} \sum_{i=1}^{N} (y_i – (w x_i + b))^2
\end{eqnarray}
$$

つまり、以下の図のように、実測値(=青点)から予測値(=青の直線)を引いたものを、二乗してすべて足し合わせるということになります。

そして、この\(E(w, b)\)が最小になる。

すなわち、実際のデータと予測値の誤差を最小にできる\(w\)と\(b\)を求めれば、もともと欲しかった単回帰モデルを求めることができます。

なお、\(w\)と\(b\)は公式が決まっているので、必要な値を求めて、あとは代入するだけでOKです。

$$
\begin{eqnarray}
w &=& \frac{\bar{xy} – \bar{y} \times \bar{x}}{\left(\bar{x^2} – \bar{x} \times \bar{x} \right)}\\
b &=& \bar{y} – w \bar{x}
\end{eqnarray}
$$

※なお、\(\displaystyle \bar{x}\)は「\(x\)の平均」という意味で、\(\displaystyle \bar{x} = \frac{1}{N} \sum_{i=1}^{N} x_i\)のことです。

1.3. 最小二乗法を手計算で計算する

\(w\)と\(b\)を求めるには、最小二乗法の式\(\displaystyle E(w, b) = \frac{1}{2} \sum_{i=1}^{N} (y_i – (w x_i + b))^2\)を最小にする必要があります。

最小二乗法の式を”最小”にできれば、それは「実際のデータと予測値の誤差を最小にできて、より正確な単回帰モデルを作成できる」ということになるわけです。

そして、最小二乗法の良いところは「解析的に求まる(=数式の変形だけで解ける)こと」です。

なので、今回は実際に、\(w\)と\(b\)の値を、数式の変形だけで求めていきましょう。

\(E(w, b)\)を最小にする\(w\)と\(b\)を求めるには、\(E(w, b)\)が最小になる、すなわち「\(E(w, b)\)を微分したらゼロになる点を求める」ということになります。

「微分してゼロになる点が最小になる」という部分は、そういうものなんだ〜と認識して頂けたらと。

まずは、最小二乗法の式\(E(w, b)\)を\(b\)で微分しましょう。

そして、\(\displaystyle \frac{\partial E}{\partial b}=0\)を計算することで、微分してゼロになる\(b\)を求めていきます。

$$
\begin{eqnarray}
\frac{\partial E}{\partial b} &=& \frac{1}{2} \sum_{i=1}^{N} 2 \times (y_i – (w x_i + b)) \times (-1)\\
&=&\sum_{i=1}^{N} (- y_i + w x_i + b)\\
&=& 0
\end{eqnarray}
$$

計算すると、上記のようになります。

\(\displaystyle \sum_{i=1}^{N}\)は、\(i\)に関係しない部分は、単純に\(N\)個足し合わせるだけなので、結局以下のような式変形になります。

$$
\begin{eqnarray}
\sum_{i=1}^{N} (- y_i + w x_i + b) &=& – \sum_{i=1}^{N} y_i + w \sum_{i=1}^{N} x_i + bN \\
&=& 0
\end{eqnarray}
$$

\(\displaystyle \sum_{i=1}^{N}\)の式を右辺に移行しましょう。すると、以下のようになります。

$$bN = \sum_{i=1}^{N} y_i – w \sum_{i=1}^{N} x_i $$

ここで両辺を\(N\)で割れば、以下の式になります。

$$b = \frac{1}{N} \sum_{i=1}^{N} y_i – w \frac{1}{N} \sum_{i=1}^{N} x_i \tag{1}$$

「\(b=\)」の式を導出できました。

ここで、\(\displaystyle \frac{1}{N} \sum_{i=1}^{N} y_i\)は、\(N\)個足し合わせて、\(N\)で割っているので、要するに「平均」になります。

平均\(\displaystyle \frac{1}{N} \sum_{i=1}^{N} y_i\)を、\(\bar{y}\)で表せば、式(1)は以下のようになります。

$$b = \bar{y} – w \bar{x} \tag{2}$$

これで\(b\)が求まる、、、と思いきや、式の中には\(w\)が入っています。

\(w\)はこれから求めないと分からない部分なので、\(b\)のときと同様に、最小二乗法の式\(E(w, b)\)を\(w\)で微分しましょう。

$$
\begin{eqnarray}
\frac{\partial E}{\partial w} &=& \frac{1}{2} \sum_{i=1}^{N} 2 \times (y_i – (w x_i + b)) \times (- x_i) \\
&=& \sum_{i=1}^{N} (- x_i y_i + w {x_i}^2 + b x_i) \\
&=& 0
\end{eqnarray}
$$

ちょっと余計な部分を前に出してあげましょう。

$$
\begin{eqnarray}
\sum_{i=1}^{N} (- x_i y_i + w {x_i}^2 + b x_i) &=& – \sum_{i=1}^{N} x_i y_i + w\sum_{i=1}^{N} {x_i}^2 + b \sum_{i=1}^{N} x_i \tag{3}\\
&=& 0
\end{eqnarray}
$$

ここで、先ほど求めた式(2)を、式(3)に代入してあげましょう。
$$
\begin{eqnarray}
– \sum_{i=1}^{N} x_i y_i + w\sum_{i=1}^{N} {x_i}^2 + b \sum_{i=1}^{N} x_i &=& – \sum_{i=1}^{N} x_i y_i + w\sum_{i=1}^{N} {x_i}^2 + (\bar{y} – w \bar{x}) \sum_{i=1}^{N} x_i \\
&=& – \sum_{i=1}^{N} x_i y_i + w\sum_{i=1}^{N} {x_i}^2 + \bar{y} \sum_{i=1}^{N} x_i – w \bar{x} \sum_{i=1}^{N} x_i \\
&=& – \sum_{i=1}^{N} x_i y_i + \bar{y} \sum_{i=1}^{N} x_i + w\sum_{i=1}^{N} {x_i}^2 – w \bar{x} \sum_{i=1}^{N} x_i \\
&=& – \sum_{i=1}^{N} x_i y_i + \bar{y} \sum_{i=1}^{N} x_i + w \left( \sum_{i=1}^{N} {x_i}^2 – \bar{x} \sum_{i=1}^{N} x_i \right) \\
&=& 0
\end{eqnarray}
$$

代入して変形していくと、上記のようになります。

この最後の式を\(N\)で割ってあげましょう。そうすると、ふたたび平均\(\displaystyle \bar{y} = \frac{1}{N} \sum_{i=1}^{N} y_i\)の形で表せます。

$$
\begin{eqnarray}
– \sum_{i=1}^{N} x_i y_i + \bar{y} \sum_{i=1}^{N} x_i + w \left( \sum_{i=1}^{N} {x_i}^2 – \bar{x} \sum_{i=1}^{N} x_i \right) &=& 0\\
– \frac{1}{N} \sum_{i=1}^{N} x_i y_i + \bar{y} \frac{1}{N} \sum_{i=1}^{N} x_i + w \left(\frac{1}{N} \sum_{i=1}^{N} {x_i}^2 – \bar{x} \frac{1}{N} \sum_{i=1}^{N} x_i \right) &=& 0 \\
-\bar{xy} + \bar{y} \times \bar{x} + w \left(\bar{x^2} – \bar{x} \times \bar{x} \right) &=& 0 \\
\end{eqnarray}
$$

かなり簡潔な形になってきましたね。あとはこれを、\(w=\)の形にするだけで完了です。

$$
\begin{eqnarray}
– \bar{xy} + \bar{y} \times \bar{x} + w \left(\bar{x^2} – \bar{x} \times \bar{x} \right) &=& 0 \\
w \left(\bar{x^2} – \bar{x} \times \bar{x} \right) &=& \bar{xy} – \bar{y} \times \bar{x}\\
w &=& \frac{\bar{xy} – \bar{y} \times \bar{x}}{\left(\bar{x^2} – \bar{x} \times \bar{x} \right)} \tag{4}
\end{eqnarray}
$$

これで\(w\)が完成しました。今回は右辺に\(b\)がないので、データ\((x, y)\)があれば手計算できますね。

\(b\)を求めるときは、

上記の手順で求めればOKです。

1.4. 最小二乗法を、さらに簡潔にする

最小二乗法で求めたい\(w\)と\(b\)は、\(E(w,b)\)を微分してゼロになることを利用して、以下のように求めることができました。

$$
\begin{eqnarray}
w &=& \frac{\bar{xy} – \bar{y} \times \bar{x}}{\left(\bar{x^2} – \bar{x} \times \bar{x} \right)}\\
b &=& \bar{y} – w \bar{x}
\end{eqnarray}
$$

この式ですが、確率の知識を使うと、さらに式を簡略化できます。

分散は、データのばらつき具合を表す指標で、「分散が大きい」とデータのばらつきが大きいということになります。

なお、分散は「ばらつき具合」なので、マイナスになることはありません。必ず\(0\)以上の値を取ります。

共分散\(Cov(x, y)\)は、\(x\)と\(y\)の関連の強さを表す指標です。

上記のような意味になります。

この分散と共分散を使えば、最終的に最小二乗法で求めたい\(w\)と\(b\)は以下のようになります。

$$
\begin{eqnarray}
w &=& \frac{Cov(x, y)}{V(x)}\\
b &=& \bar{y} – w \bar{x}
\end{eqnarray}
$$

Pythonを使った最小二乗法の実装方法【Pandas】

前回の単回帰モデルで、以下の式を導出しました。

$$y = -36.476 + 9.3764 x$$

今回学習した「最小二乗法」を使って、重みとバイアスの部分を求めてみましょう。

使うコードは以下のとおりです。

# ### 2.1 ライブラリの読み込み
import pandas as pd
from sklearn.datasets import load_boston

# ### 2.2 データの読み込み
boston = load_boston()
df = pd.DataFrame(boston.data, columns=boston.feature_names)
df['MEDV'] = boston.target

print(df.head())

# ### 2.3 データの準備
x = df['RM']
y = df['MEDV']

print('RMのデータ数:', x.shape)
print('MEDVのデータ数:', y.shape)

# ### 2.4 実際に計算する
calc_df = pd.DataFrame({'x':x, 'y':y})
print(calc_df.shape)
calc_df.head()

calc_df = calc_df.assign(
    x2 = lambda df:df.x*df.x,
    xy = lambda df:df.x*df.y
)

x_mean = calc_df.x.mean()
y_mean = calc_df.y.mean()
x2_mean = calc_df.x2.mean()
xy_mean = calc_df.xy.mean()

print(x_mean)
print(y_mean)
print(x2_mean)
print(xy_mean)

w = (xy_mean - x_mean*y_mean)/(x2_mean - x_mean * x_mean)
print(w)

b = y_mean - w * x_mean
print(b)

STEP① : ライブラリの読み込み

import pandas as pd
from sklearn.datasets import load_boston

まずは、今回必要になるライブラリをインポートします。

STEP② : データの読み込み

boston = load_boston()
df = pd.DataFrame(boston.data, columns=boston.feature_names)
df['MEDV'] = boston.target

print(df.head())

単回帰モデルのときと同じ重み\(w\)・バイアス\(b\)になるか確認したいので、今回も「ボストンの住宅価格」データを取得していきます。

前回の単回帰モデルでは、RM(=住宅に含まれる平均部屋数)とMEDV(=住宅価格)を使ってい他ので、今回もこの2つだけ使っていきます。

STEP③ : データの準備

x = df['RM']
y = df['MEDV']

print('RMのデータ数:', x.shape)
print('MEDVのデータ数:', y.shape)

それでは、単回帰モデルの実装と同じ変数だけ抽出します。

特徴量(=機械学習で使う変数)はRM、予測対象はMEDVでしたね。

STEP④ : 実際に計算する

calc_df = pd.DataFrame({'x':x, 'y':y})
print(calc_df.shape)
calc_df.head()

calc_df = calc_df.assign(
    x2 = lambda df:df.x*df.x,
    xy = lambda df:df.x*df.y
)

x_mean = calc_df.x.mean()
y_mean = calc_df.y.mean()
x2_mean = calc_df.x2.mean()
xy_mean = calc_df.xy.mean()

print(x_mean)
print(y_mean)
print(x2_mean)
print(xy_mean)

w = (xy_mean - x_mean*y_mean)/(x2_mean - x_mean * x_mean)
print(w)

b = y_mean - w * x_mean
print(b)

それでは、実際に\(w\)と\(b\)に必要な値を求めていきます。

\(w\)と\(b\)に必要な値は、以下の式に含まれるものです。

$$
\begin{eqnarray}
w &=& \frac{\bar{xy} – \bar{y} \times \bar{x}}{\left(\bar{x^2} – \bar{x} \times \bar{x} \right)}\\
b &=& \bar{y} – w \bar{x}
\end{eqnarray}
$$

結局のところ、上記の値を求めれば、\(w\)と\(b\)は導出できますね。

ただ、現状だと\(x^2\)と\(xy\)のカラムがないので、Pandasを使って追加してあげましょう。

calc_df = calc_df.assign(
    x2 = lambda df:df.x*df.x,
    xy = lambda df:df.x*df.y
)

assign()を使うと、新しい列を追加できます。

その後、今回必要になる4つの値をmean()メソッドを使って計算します。

x_mean = calc_df.x.mean()
y_mean = calc_df.y.mean()
x2_mean = calc_df.x2.mean()
xy_mean = calc_df.xy.mean()

print(x_mean)
print(y_mean)
print(x2_mean)
print(xy_mean)

これで、最小二乗法で必要になる4つの値を取得できました。最後に、以下の式に各値を代入すれば完了です。

$$
\begin{eqnarray}
w &=& \frac{\bar{xy} – \bar{y} \times \bar{x}}{\left(\bar{x^2} – \bar{x} \times \bar{x} \right)}\\
b &=& \bar{y} – w \bar{x}
\end{eqnarray}
$$

w = (xy_mean - x_mean*y_mean)/(x2_mean - x_mean * x_mean)
print(w)

b = y_mean - w * x_mean
print(b)

計算すると、以下のような結果になるかと思います。

単回帰モデルをプログラムで作成したときは、\(y = -36.476 + 9.3764 x$\)でしたが、手計算で計算すると若干の誤差が生じますね。

これはコンピューターの問題なので、気にしなくて大丈夫です。

このような面倒な計算を、機械学習ライブラリScikit-learnを使えば自動的に計算してくれます。ライブラリは偉大ですね。

【機械学習入門③】最小二乗法まとめ

というわけで、今回は以上です。

単回帰モデルの裏側で動いている最小二乗法について、深く理解できたのではないでしょうか。

次回は、単回帰モデルと同じ線形回帰の中でも、「重回帰モデル」について学習していきます。

【機械学習入門④】線形回帰・重回帰モデルをPythonの実装付きで理解！ 前回の記事『【機械学習入門③】最小二乗法とは？手計算とPythonを使った実装を解説【初心者向け】』では、単回帰モデルの導出方法...

単回帰モデルでは使う特徴量が1つだけでしたが、重回帰モデルでは複数の特徴量をモデルに組み込めます。

なので、より機械学習っぽくなって、勉強していて楽しいはずです。ぜひ取り組んでみてください！

なお、『機械学習入門シリーズ』は長いので、休みながら学習していきましょう(*’ω’*)

休憩するときは、記事をツイートしておくと便利ですよ！

※メンションに気が付いたら、拡散する可能性高いです(｀・ω・´)！

それでは、この辺で。

おすすめの記事

【機械学習入門④】線形回帰・重回帰モデルをPythonの実装付きで理解！ 前回の記事『【機械学習入門③】最小二乗法とは？手計算とPythonを使った実装を解説【初心者向け】』では、単回帰モデルの導出方法...