概要

ここでは実験データがtxtファイルやdatファイルで得られた状態からデータをフィッティングしてパラメーターを決めることを目指します。関数が線形、非線形に関わらず使えるようにここではPythonのscipyパッケージに入っているscipy.optimize モジュールの一部のcurve_fitというモジュールを使う。

一次関数のフィッティング

まず必要なモジュールを読み込みます。

# フィッティングに必要
from scipy.optimize import curve_fit
import numpy as np
# 図示に必要
import matplotlib.pyplot as plt

以下のデータをフィッティングしてみる。

x = np.asarray([1,2,4,6,7])
y = np.asarray([1,3,3,5,4])
plt.scatter(x, y, c='k')
plt.show()

フィッティングする関数形を定義して実際にフィッティングを行う。

# フィッティングしたい関数式を定義する。ここでは一次関数。
def linear_fit(x, a, b):
    return a*x + b
# フィッティングを行う。
# prameter_initial = np.array([1,2])←ここで初期値を設定するがここでは飛ばす。
popt, pcov = curve_fit( linear_fit, x, y)#, p0= prameter_initial)

ここまででフィッティングが行われ、第1の返り値poptの中にパラメータa,bの結果が格納されている。今回は一次関数の傾きaと切片bの二つのパラメーターがpoptに[a,b]=[popt[0],popt[1]]と入っている。どんな値になったのか確認する。

print (popt[0],popt[1])

多分、a=0.4999…、b=1.2000…と求まるはずである。この値は妥当なのか確認するために重ねてプロットする。なんかいい感じに見えなくもない。

c=popt[0]
d=popt[1]
f = c*x + d #f = linear_fit(x,popt[0],popt[1])かf = linear_fit(x,*popt)でもOK。
plt.scatter(x, y, c='k')
plt.plot(x, f, c='r')
plt.show()

二次関数のフィッティング

同様に必要なモジュールを読み込みます。

# フィッティングに必要
from scipy.optimize import curve_fit
import numpy as np
# 図示に必要
import matplotlib.pyplot as plt

以下のデータをフィッティングしてみる。

x = np.asarray([1,2,2.5,4,5,6,7,8,9,10])
y = np.asarray([10,8,5,2,1,1.2,2.5,5,7,10])
plt.scatter(x, y, c='k')
plt.show()

フィッティングする関数形を定義して実際にフィッティングを行う。

# フィッティングしたい関数式を定義する。ここでは二次関数。
def two_fit(x,a,b):
    return  a*(x-b)**2

# フィッティングを行う。初期値も設定。
prameter_initial = np.array([1,5])
popt, pcov = curve_fit(two_fit, x, y, p0= prameter_initial)

ここまででフィッティングが行われ、第1の返り値poptの中にパラメータa,bの結果が格納されている。今回は一次関数の傾きaと切片bの二つのパラメーターがpoptに[a,b]=[popt[0],popt[1]]と入っている。どんな値になったのか確認する。

print (popt[0],popt[1])

多分、a=0.5404046977156955、b=5.488472528172034と求まるはずである。この値は妥当なのか確認するために重ねてプロットする。一次関数のところでは説明しなかったが、この時のプロットの方法は二通りある。まず、一次関数の時と同じ方法でやってみる。

c=popt[0]
d=popt[1]
f = c*(x-d)**2 #f = two_fit(x,popt[0],popt[1])#こんな感じでも書けるし、f = two_fit(x,*popt)でもOK。
plt.scatter(x, y, c='k')
plt.plot(x, f, c='r')
plt.show()

一次関数の時はわからなかったが、カクカクしているのがわかると思う。データ数が多ければ、これは気にならないがこのくらいだと目で見た感じがなんか嫌なので、二つ目の方法として滑らかに書いてみる。

# 描画範囲の指定
# x1 = np.arange(x軸の最小値, x軸の最大値, 刻み)
x1 = np.arange(1, 10, 0.000001)
c=popt[0]
d=popt[1]
f = c*(x1-d)**2
plt.scatter(x, y, c='k')
plt.plot(x1, f, c='r')
plt.show()

いい感じに滑らかになる。描画範囲としてx1を作って刻みを細かくする。
ここまでで基本的なフィッティングの方法と描画の方法がわかったので、最後にtxtファイルやdatファイルの実験データを読み込んでフィッティングしていく。

実験データのフィッティング

最後に以上を踏まえてtxtファイル、datファイルからデータを読み込んで、フィッティングを行いパラメーターを決めていく。一番最後にコード全体を載せておく。