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