フィッティング３(最小二乗法，optimizeの利用)

　前項では，最小二乗法がうまく機能して，測定データによく沿う傾き，切片を決めることができた． ところでpythonでは様々な機能がモジュールとして提供されているはずであった．最小二乗法もありそうである．

optimize

検索などをして調べてみると，scipyのoptimizeにleastsqという機能があるのが見つかる． 　>> SciPy Optimize

optimize.leastsqを使うとコードは以下のようになる．

import numpy as np
import matplotlib.pyplot as plt
from scipy import optimize

def f(x):
    return 5.*x + 10.

#Pusedo measured data by random number
xdata = np.linspace(-10, 10, num=101)
np.random.seed(1234)
ydata = f(xdata) + 5.*np.random.randn(xdata.size)

#Least squares method with scipy.optimize
def fit_func(parameter,x,y):
    a = parameter[0]
    b = parameter[1]
    residual = y-(a*x+b)
    return residual

parameter0 = [0.,0.]
result = optimize.leastsq(fit_func,parameter0,args=(xdata,ydata))
print(result)
a_fit=result[0][0]
b_fit=result[0][1]

print(a_fit,b_fit)

#PLot
plt.figure(figsize=(8,5))
plt.plot(xdata,ydata,'bo', label='Exp.')
plt.plot(xdata,a_fit*xdata+b_fit,'k-', label='fitted line', linewidth=10, alpha=0.3)
plt.plot(xdata,f(xdata),'r-', label='mother line', linewidth=3)
plt.xlabel('Time')
plt.ylabel('Verocity')
plt.legend(loc='best',fancybox=True, shadow=True)
plt.grid(True)
plt.show()

１４～１８行目：近似する関数形を記述する．a,bが決めたい係数である．

２１行目：主要部分．a,bの初期値parameter0から初めて，(xdata,ydata)に対してfit_funcが最小になるような，a,bを返してくれる．

２３，２４行目：得られた結果をa_fit,b_fitに格納．

得られたa,bの値は，

(4.9868427708984067, 10.188239312198032)

やはり測定データによく一致した直線が得られている．また当然だが，測定データを生成するのにつかった　\(f(x) = 5x + 10\)　にもぴたりと一致している．

optimize.leastsqは，１次関数以外にも使える．「データの描画」項で読み込んだデータはいかにも２次関数で近似できそうである． 　\(f(x) = ax^2 + bx + c\)

import numpy as np
import matplotlib.pyplot as plt
from scipy import optimize

file_id = 'sample_data.csv'
#file_path = '/Users/kentaro/work_space/python/'
file_path = './'

rfile = file_path + file_id
data = np.loadtxt(rfile, comments='#' ,delimiter=',')

x_csv = data[:,0]
y_csv = data[:,1]

#Least squares method with scipy.optimize
def fit_func(parameter,x,y):
    a = parameter[0]
    b = parameter[1]
    c = parameter[2]
    residual = y-(a*x**2+b*x+c)
    return residual

parameter0 = [0.,0.,0.]
result = optimize.leastsq(fit_func,parameter0,args=(x_csv,y_csv))
a_fit=result[0][0]
b_fit=result[0][1]
c_fit=result[0][2]
print(a_fit,b_fit,c_fit)

#PLot
plt.figure(figsize=(8,5))
plt.plot(x_csv,y_csv,'bo', label='Exp.')
plt.plot(x_csv,a_fit*x_csv**2+b_fit*x_csv+c_fit,'k-', label='fitted parabora', linewidth=10, alpha=0.3)
plt.xlabel('X-axis')
plt.ylabel('Y-axis')
plt.legend(loc='best',fancybox=True, shadow=True)
plt.grid(True)
plt.show()

(0.88720506883977268, 0.41856724361772535, 5.7365571795694734)

もともと測定データを作るのにつかった一次関数の傾きa=5,切片b=10であった．これにとても近い数値が得られている． グラフを描けば以下の通り．最小二乗法により得られた１次関数がデータによく沿っている．