numpyで最小二乗法フィッティングを行う。

Pythonの行列演算ライブラリnumpyを用いて,単一変数xによって説明されるyについて
最小二乗フィッティングを行いたいと思います。

まず、numpyを用いて確率的に3次のグラフを描画します。

グラフの描画

# モジュールのインポート
import numpy as np
import matplotlib.pyplot as plt

# 説明変数(1次元)
x = np.arange(-3,7,0.5)
# 応答変数(説明変数の3次元関数とし、正規分布に基づいて乱数的に生成した。
y = 10*np.random.rand()+x * np.random.rand() + 
    2*np.random.rand()*x**2  +x**3

# 描画
plt.scatter(x,y)
plt.show()

最小二乗法では、データと予測値のL2ノルムを
誤差(Error)と解釈し、それを最小となるように回帰線の係数を求めます。
予測したいデータをy、回帰線を$p$とすると、Errorは

Error = \sum_{i=1}^{N}(y_i - p_i)

となります。
このErrorを最小化することが最小二乗法回帰のゴールです。

また、$p_i$は以下の様にn次式で表されます。

1次式\\
p_i = a_1 x + a_0 \\
2次式\\
p_i = a_2 x^2 + a_1 x + a_0 \\
3次式 \\
p_i = a_3 x^3 +  a_2 x^2 + a_1 x + a_0 \\
n次式 \\
p_i = a_n x^n + ...  a_2 x^2 + a_1 x + a_0\\

今回は、numpyのpolyfit関数を用いてフィッティング後の式の係数$A=(a_0,a_1,..a_N)$を求めたいと思います。
polyfit関数により回帰式の係数を求めます。
その後、係数をn次式に当てはめて
回帰式を求めますが、複雑になってきたら
polyfit1d関数などを用いるのも便利です。

フィッティングとそれにより得られたn次式の描画

# 1次式
coef_1 = np.polyfit(x,y,1) #係数
y_pred_1 = coef_1[0]*x+ coef_1[1] #フィッティング関数

# 2次式
coef_2 = np.polyfit(x,y,2) 
y_pred_2 = coef_2[0]*x**2+ coef_2[1]*x + coef_2[2] 

# 3次式
coef_3 = np.polyfit(x,y,3) 
y_pred_3 = np.poly1d(coef_3)(x) #np.poly1d,求めた係数coef_3を自動で式に当てはめる。

# 描画

plt.scatter(x,y,label="raw_data") #元のデータ
plt.plot(x,y_pred_1,label="d=1") #1次式
plt.plot(x,y_pred_2,label="d=2") #2次式
plt.plot(x,y_pred_3,label = "d=3") #3次式
plt.legend(loc="upper left")
plt.title("least square fitting")
plt.show()

今回は3次式によるフィッティングが良さそうです。
高次になるほど誤差は小さくなりがちですが、
得られたデータセットにのみ依存する過学習
には注意しましょう。