はじめに
どうも。今回は最小二乗法をPythonで実装する方法を記事にしてみました。先に述べておきますが、以下の数学の知識がないとチンプンカンプンだと思いますのであしからず。
必要な知識
・1次関数や2次関数などの基本的な関数
・微分と偏微分
・総数（Σ）、平均について
高校生までの数学の知識があれば、特に問題はないと思います。では、早速どんなものか見てみましょうか。
最小二乗法とは
こんなデータがあったとします。（下のグラフ）
「それぞれの点からの誤差を最小にした1次関数を描いて下さい」と言われても、無理ですよね。これを可能にする方法が最小二乗法（最小自乗法）です。測定データ$y$がモデル関数（元となる関数）$f(x)$と誤差$\varepsilon$の和であるときに使用できます。式で表すと下のようになります。
y = f(x) + \varepsilon
ランダムなデータでは、最小二乗法を使うことは出来ないということです。
さて、今回のデータは1次関数に近似できそうなのでそれをモデル関数として、近似してみましょう。なお、1次関数の式は以下のようにおいておきます。
f(x) = ax + b
最小二乗法で$a$と$b$の適切なパラメータを求めます。まずは、名前にもなっている実際のデータとモデル関数の差の二乗の総和を求めます。式にすると下のようになります。
J=\sum_{i=1}^n(y_i-f(x_i))^2=\sum_{i=1}^n(y_i-ax_i-b)^2
誤差の二乗を求める理由は、誤差が正負両方あった場合お互いを打ち消し合って正確に求めることができなくなるため、二次関数のようにして防ぐために行います。
誤差の二乗の総和$J$を$a$と$b$で偏微分します。偏微分したものが以下のようになります。
\frac{\partial J}{\partial a}=-2\sum_{i=1}^nx_i(y_i-ax_i-b)=-2(\sum_{i=1}^n x_iy_i - a\sum_{i=1}^n x_i^2 - b\sum_{i=1}^n x_i)\tag{1}
\frac{\partial J}{\partial b}=-2\sum_{i=1}^n(y_i-ax_i-b)=-2(\sum_{i=1}^n y_i - a\sum_{i=1}^n x_i - nb)\tag{2}
この2式が0となるときの$a$と$b$を求めます。まず、式(2)を変形して$b=$の式にします。
b=\frac{1}{n}\sum_{i=1}^n y_i - \frac{a}{n}\sum_{i=1}^n x_i\tag{3}
式(3)を式(1)に代入し、変数を$a$だけにして$a=$の式にします。
a=\frac{n\sum_{i=1}^n x_iy_i-\sum_{i=1}^n x_i \sum_{i=1}^n y_i}{n\sum_{i=1}^n x_i^2-(\sum_{i=1}^n x_i)^2}
あとは、これを式(3)に代入して$a$と$b$のパラメータを求めることが出来ますが...、Σが入りすぎてごちゃごちゃになっていますね。少し見やすいように文字を置いてみましょう。
a=\frac{nXY_{sum}-X_{sum}Y_{sum}}{nX^2_{sum}-(X_{sum})^2}\tag{4}
b=Y_{ave}-aX_{ave}\tag{5}
どうでしょうか。少しは見やすくなったでしょうか。文字の意味はそれぞれ下のようになります。
- $XY_{sum}$ $x$と$y$の積の総和
- $X_{sum}$ $x$の総和
- $Y_{sum}$ $y$の総和
- $X^2_{sum}$ $x^2$の総和
- $Y_{ave}$ $y$の平均
- $X_{ave}$ $x$の平均
ちなみに、式(4)の分子と分母を$n^2$で割ると、全てそれぞれの平均で求めることも出来ます。
a=\frac{XY_{ave}-X_{ave}Y_{ave}}{X^2_{ave}-(X_{ave})^2}\tag{6}
では、原理もわかったことでPythonで計算していきましょう。
Pythonで実装
データは最初に示したグラフのデータを使用しました。下にソースコードを示します。なお、数値計算にはnumpyを使用しました。
import numpy as np
import matplotlib.pyplot as plt
#使用するデータ
X = np.array([1, 2, 3, 4, 5, 6, 7, 8, 9])
Y = np.array([2, 5, 3, 7, 12, 11, 15, 16, 19])
plt.scatter(X, Y, color = 'red')
#データの数
n = Y.size
#XとYの積
XY = X * Y
#Xの各要素の二乗
X_2 = X * X
#パラメータの導出
a = (n * np.sum(XY) - np.sum(X) * np.sum(Y)) / (n * np.sum(X_2) - (np.sum(X)**2))
b = np.average(Y) - a * np.average(X)
print('a', a)
print('b', b)
Y = a * X + b
plt.plot(X, Y)
plt.show()
実行結果
コマンドに出力された結果を下に示します。
a 2.15
b -0.75
傾きは2.15、切片は-0.75となったようです。また、グラフを下に示します。
見ての通り、うまく1次関数を描けました。ちなみに、誤差の変化は$J(a,b)$の2変数関数として表すことが出来ます。3Dのグラフに表したものを下に示します。
赤点の座標$(x, y, z)=(2.15, -0.75, 16.6)$が今回求めた$a$と$b$の値と、誤差$J$の値を示しています。グラフでも分かるように、この座標が極小値となっています。一応、ソースコードを載せておきますので、ぜひ実行して一度自分の目で確かめてみて下さい。
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
fig = plt.figure()
ax = Axes3D(fig)
#J関数の定義
def Er(a, b, y_2, xy, y, x_2, x, n):
return y_2-2*a*xy-2*b*y+a**2*x_2+2*a*b*x+n*b**2
#使用するデータ
X = np.array([1, 2, 3, 4, 5, 6, 7, 8, 9])
Y = np.array([2, 5, 3, 7, 12, 11, 15, 16, 19])
#データの数
n = Y.size
#XとYの積
XY = X * Y
#Xの各要素の二乗
X_2 = X * X
#パラメータの導出
a = (n * np.sum(XY) - np.sum(X) * np.sum(Y)) / (n * np.sum(X_2) - (np.sum(X)**2))
b = np.average(Y) - a * np.average(X)
print(a)
print(b)
A = np.arange(a, a+0.1, 0.00001)
B = np.arange(b, b+0.1, 0.00001)
A, B = np.meshgrid(A, B)
J = Er(A, B, np.sum(Y*Y), np.sum(XY), np.sum(Y), np.sum(X_2), np.sum(X), n)
ax.plot_wireframe(A, B, J)
j = Er(a, b, np.sum(Y*Y), np.sum(XY), np.sum(Y), np.sum(X_2), np.sum(X), n)
print(j)
ax.scatter(a, b, j, color='red')
ax.set_xlabel("A")
ax.set_ylabel("B")
ax.set_zlabel("J")
plt.show()
まとめ
今回は最小二乗法をPythonで実装してみたという話でした。最小二乗法は統計学の内の一つであり、AIなどを学ぶ際には必ず必要な分野です。私も現在勉強していますので、またメモ程度に残していこうと思います。ご質問などありましたら、コメントまでよろしくお願いします。