%load_ext sympy.interactive.ipythonprinting
import numpy as np
from matplotlib import pyplot as plt
from sympy.matrices import *
from sympy import *
点の座標と変数a, bを初期化する.
points = [(4, -17), (15, -4), (30, -7), (100, 50)]
a, b = symbols("a b")
ここでの目的は,$$ J = \frac{1}{2} \sum_{\alpha=1}^{N} (y_{\alpha} - (ax_{\alpha} + b))^2 \to min$$ となるようなa, bを見つけることである.
この式をPythonで表現すると次のようになる.
J = sum([(y - (x*a + b))**2 for x, y in points]) / 2
J
ここから$ \frac{\partial J}{\partial a}$と $\frac{\partial J}{\partial b}$を計算する.
pda = diff(J, a)
pda
pdb = diff(J, b)
pdb
連立方程式 $$ 1141a + 149b - 4662 = 0$$ $$ 149a + 4b - 22 = 0 $$ を解いて、a,bを得る.
solutions = solve([pda, pdb], [a, b])
a = solutions[a].evalf()
a
b = solutions[b].evalf()
b
よって最小二乗法によって求められた曲線は $y = 0.687 x - 20.1$ である.
pxs, pys = zip(*points)
xs = np.linspace(0, 100, 1000)
ys = xs * a + b
plt.plot(pxs, pys, ".")
plt.plot(xs, ys,)
[Line2D(_line1)]