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

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

optimize

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

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

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
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()

14~18行目:近似する関数形を記述する.a,bが決めたい係数である.

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

23,24行目:得られた結果をa_fit,b_fitに格納.

得られたa,bの値は,

(4.9868427708984067, 10.188239312198032)

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

fitted line by optimize

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

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
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であった.これにとても近い数値が得られている. グラフを描けば以下の通り.最小二乗法により得られた1次関数がデータによく沿っている.

fitted parabora

Table Of Contents

Previous topic

フィッティング2(最小二乗法)

Next topic

周波数解析

This Page