Python
scipy
numpy

カーブフィッティング手法 scipy.optimize.curve_fit の使い方を理解する

Pythonを使ってカーブフィッティング(曲線近似)する方法として、 scipy.optimize.curve_fit を使う方法がありますが、使い方が少し理解しにくいと思ったので整理してみました。


用いる実験値

Numpy.polyfit を使ったカーブフィッティング で用いたのと同じデータを用いて、比較してみましょう。

x_observed = [9, 28, 38, 58, 88, 98, 108, 118, 128, 138, 148, 158, 168, 178, 188, 198, 208, 218, 228, 238, 278, 288, 298]

y_observed = [51, 80, 112, 294, 286, 110, 59, 70, 56, 70, 104, 59, 59, 72, 87, 99, 64, 60, 74, 151, 157, 57, 83]

Python の List 型から、Numpy の Array 型に変換しておきましょう。

import numpy as np

x_observed = np.array(x_observed)
y_observed = np.array(y_observed)


1次式近似

まず、1次式に近似してみましょう。関数 func1 を定義します。a と b が、求めたい値になります。

def func1(X, a, b): # 1次式近似

Y = a + b * X
return Y

関数 func1 の使用例はこちらです。

func1(x_observed, 2, 2)

array([ 20,  58,  78, 118, 178, 198, 218, 238, 258, 278, 298, 318, 338,

358, 378, 398, 418, 438, 458, 478, 558, 578, 598])

さて、scipy.optimize.curve_fit を使って近似してみましょう。

from scipy.optimize import curve_fit  

popt, pcov = curve_fit(func1,x_observed,y_observed) # poptは最適推定値、pcovは共分散
popt

array([125.77023172,  -0.1605313 ])

ここで得られた popt が最適推定値を格納しています。Numpy.polyfit を使ったカーブフィッティング で得られた推定値と比較してみましょう。


2次式近似

次に、2次式に近似してみましょう。関数 func2 を定義します。a と b と c が、求めたい値になります。

def func2(X, a, b, c): # 2次式近似

Y = a + b * X + c * X ** 2
return Y

関数 func2 の使用例はこちらです。

func2(x_observed, 2, 2, 2)

array([   182,   1626,   2966,   6846,  15666,  19406,  23546,  28086,

33026, 38366, 44106, 50246, 56786, 63726, 71066, 78806,
86946, 95486, 104426, 113766, 155126, 166466, 178206])

さて、scipy.optimize.curve_fit を使って近似してみましょう。

from scipy.optimize import curve_fit  

popt, pcov = curve_fit(func2,x_observed,y_observed) # poptは最適推定値、pcovは共分散
popt

array([ 1.39787684e+02, -4.07809516e-01,  7.97183226e-04])

ここで得られた popt が最適推定値を格納しています。Numpy.polyfit を使ったカーブフィッティング で得られた推定値と比較してみましょう。

次のようにすれば、最適推定値を用いた近似曲線が描けます。

func2(x_observed, 1.39787684e+02, -4.07809516e-01,  7.97183226e-04)

array([136.1819702 , 128.9940092 , 125.44205497, 118.81645644,

110.07383349, 107.47849913, 105.04260142, 102.76614035,
100.64911593, 98.69152815, 96.89337701, 95.25466253,
93.77538468, 92.45554348, 91.29513893, 90.29417102,
89.45263976, 88.77054514, 88.24788717, 87.88466585,
88.02614699, 88.46010889, 89.05350743])


3次式近似

同様に、3次式に近似してみましょう。関数 func3 を定義します。a と b と c と d が、求めたい値になります。

def func3(X, a, b, c, d): # 3次式近似

Y = a + b * X + c * X ** 2 + d * X ** 3
return Y

さて、scipy.optimize.curve_fit を使って近似してみましょう。

from scipy.optimize import curve_fit  

popt, pcov = curve_fit(func3,x_observed,y_observed) # poptは最適推定値、pcovは共分散
popt

array([ 7.84214107e+01,  1.88213104e+00, -1.74165777e-02,  3.89638087e-05])

ここで得られた popt が最適推定値を格納しています。Numpy.polyfit を使ったカーブフィッティング で得られた推定値と比較してみましょう。

次のようにすれば、最適推定値を用いた近似曲線が描けます。

func3(x_observed, 7.84214107e+01,  1.88213104e+00, -1.74165777e-02,  3.89638087e-05)

array([ 93.97825188, 118.32181643, 126.93087413, 136.59795028,

135.72770915, 132.27386543, 127.62777811, 122.02323006,
115.69400413, 108.87388316, 101.79665001, 94.69608754,
87.80597859, 81.36010602, 75.59225267, 70.73620141,
67.02573508, 64.69463654, 63.97668864, 65.10567422,
92.76660851, 106.63700433, 123.75703075])


何次式になっても使える汎用的な関数を作りたい

今まで1次式、2次式、3次式と別々の関数を作りましたが、これだと面倒すぎるので、何次式になっても使える汎用的な関数を作りましょう。何次式に近似するかはパラメータの数で自動的に判定するようにします。

import numpy as np

def func(X, *params):
Y = np.zeros_like(X)
for i, param in enumerate(params):
Y = Y + np.array(param * X ** i)
return Y

次の2つのコードを実行して、上記の func2 による計算結果と、func3 による計算結果と同じであることを確認しましょう。

func(x_observed, 1.39787684e+02, -4.07809516e-01,  7.97183226e-04)

array([136.1819702 , 128.9940092 , 125.44205497, 118.81645644,

110.07383349, 107.47849913, 105.04260142, 102.76614035,
100.64911593, 98.69152815, 96.89337701, 95.25466253,
93.77538468, 92.45554348, 91.29513893, 90.29417102,
89.45263976, 88.77054514, 88.24788717, 87.88466585,
88.02614699, 88.46010889, 89.05350743])

func(x_observed, 7.84214107e+01,  1.88213104e+00, -1.74165777e-02,  3.89638087e-05)

array([ 93.97825188, 118.32181643, 126.93087413, 136.59795028,

135.72770915, 132.27386543, 127.62777811, 122.02323006,
115.69400413, 108.87388316, 101.79665001, 94.69608754,
87.80597859, 81.36010602, 75.59225267, 70.73620141,
67.02573508, 64.69463654, 63.97668864, 65.10567422,
92.76660851, 106.63700433, 123.75703075])


何次式になっても使える汎用的な関数による多項式近似

ここまでくれば、何次式にでも多項式近似できます。ただし、パラメータの数を指定するために p0= で初期値を必要数だけ入力する必要が生じます。

from scipy.optimize import curve_fit  

popt, pcov = curve_fit(func,x_observed,y_observed, p0=[1, 1])
popt

array([125.77023172,  -0.1605313 ])

from scipy.optimize import curve_fit  

popt, pcov = curve_fit(func,x_observed,y_observed, p0=[1, 1, 1])
popt

array([ 1.39787684e+02, -4.07809516e-01,  7.97183226e-04])

from scipy.optimize import curve_fit  

popt, pcov = curve_fit(func,x_observed,y_observed, p0=[1, 1, 1, 1])
popt

array([ 7.84214107e+01,  1.88213104e+00, -1.74165777e-02,  3.89638087e-05])

from scipy.optimize import curve_fit  

popt, pcov = curve_fit(func,x_observed,y_observed, p0=[1, 1, 1, 1, 1])
popt

array([-7.54495613e+01,  1.13179580e+01, -1.50591743e-01,  7.02970546e-04,

-1.07313869e-06])


Numpy.polyfit を使ったカーブフィッティングとの比較

上記の計算結果を、Numpy.polyfit を使ったカーブフィッティング で得られた推定値と比較してみましょう。

Numpy.polyfit は、多項式近似するだけなら、便利で使いやすいですが、多項式近似しかできません。もっと他の関数で近似したい場合は、scipy.optimize.curve_fit の使い方を理解するのが良いでしょう。


少ない観測値を補間してから、正規分布の線形和で近似する

少ない観測値を補間してから、正規分布の線形和で近似する では、scipy.optimize.curve_fitを使って、正規分布(ガウス関数)の線形和で近似する方法について記述していますが、ここで説明した scipy.optimize.curve_fit の使い方を踏まえれば、理解が容易になるのではないかと期待します。