Yamamoto's Laboratory
フィッティング

SciPyフィッティング (fitting)

SciPy モジュール「scipy.optimize.curve_fit」を使ったフィッティング (fitthing) の方法を示します.

目次


はじめに

フィッティングとは

実際に観測されるデータは,誤差を持ちます.この誤差を持ったデータから,推定される曲線 (関数) を得る手法をフィティングと言います.具体例を示しましょう.図のようにバネが吊り下げられています.そこに,重りを1[kg]から10[kg]まで変化させて,バネの長さを測定します.この測定結果は以下のとおりです.

この測定結果から,重さがゼロの時のバネの長さ(初期長さ)とバネ定数を求めることを考えます.測定結果をプロットすると,図が得られます.このプロットから最も確からしい直線を引くことができれば,最も確からしい初期長さとバネ定数を得ることができます.バネの長さ\(y\)と重りの重さ\(x\)には, \begin{align} y=a+bx \end{align} の関係があります.\(a\)が初期長さ,\(b\)がバネ定数です.図が最も確からしい直線です.これは,誤差の二乗である\(\sum \left\[f(x_i)-y_i\right\]^2\)を最小にした直線です.

フィッティング方法

Python には,フィッティングのためのモジュール「scipy.optimize.curve_fit」があります.これを使うと容易に誤差を持つデータを任意の関数でフィッティングすることができます.これを使うためにのステップは,次のとおりです.

  1. モジュールをインポートします.よく使われるコマンドは,「from scipy.optimize import curve_fit」です.
  2. フィッティングする関数を定義します.関数の第一引数が独立変数,第二引数以降はフィッティングパラメーターです.例えば,「\(f(x,\, a\, b)\)」のようにです.
  3. データを準備します.例えば,\((x,\,y)\)のようにです.\(x\)と\(y\)のいずれも,シーケンスデータで同じ個数の実数から成ります.
  4. フィッティングを行います.コマンドは「popt, pcov = curve_fit(f, x, y)です.popt が最適推定値\((a,\,b)\),pcovは共分散です.

フィッティングに必要な手順を踏んでおり,この流れは簡単に理解できます.

使ってみよう

ここでは,誤差を持ったデータとして,\(y=\sin(x)\exp(-x/5)+\mathrm{Noise}\) を考えます.\(\mathrm{Noise}\) が観測で,それは標準偏差1を持った正規分布に0.05を乗じた値とします.これを,フィッティングパラメーター\((a,\,b,\,c)\)を持った関数 \(f(x)=a\sin(x)\exp(-bx)+c\) でフィッティングします.ノイズの無い元の関数と比較すると,\((a,\,b,\,c)\simeq(1,\,0.2,\,0)\) となることが分かります.サンプル数が増えると,この値になるはずです.

それでは,実際のプログラムを以下に示します.先に示したステップのモジュールのインポートは 004 行,フィッティング関数の定義は 056 行,データの準備は 065 と 067 行,フィッティングの実行は 023 行です.

このプログラムはメインルーチンで,(1) データ\((x,\,y)\)のデータを作成し,(2) インスタンス fit を生成し.(3) フィッティングの操作を行い,(4)プロットを描きます.メインルーチンを見れば分かるでしょう.

データを関数でフィッティングするプログラム(start.py)

001   # -*- coding: utf-8 -*-
002   import numpy as np
003   import matplotlib.pyplot as plt
004   from scipy.optimize import curve_fit
005   
006   # ===========================================================
007   # フィッティングのクラス
008   # ===========================================================
009   class FITTING(object):
010       # -------------------------------------------------------
011       # コンストラクター
012       # -------------------------------------------------------
013       def __init__(self, function, x, y):
014           self.f = function
015           self.x = x
016           self.y = y
017   
018   
019       # -------------------------------------------------------
020       # フィッティングの実行
021       # -------------------------------------------------------
022       def do_fitting(self):
023           self.popt, self.pcov = curve_fit(self.f, self.x, self.y)
024   
025           print('a: {0:e}\nb: {1:e}\nc: {2:e}'.\
026                 format(self.popt[0], self.popt[1], self.popt[2]))
027   
028   
029       # -------------------------------------------------------
030       # フィッティング結果のプロット
031       # -------------------------------------------------------
032       def plot(self, Nx=65):
033           xmin, xmax = min(self.x), max(self.x)
034           xp = np.linspace(xmin, xmax, Nx)
035           fig = plt.figure()
036           plot = fig.add_subplot(1,1,1)
037           plot.set_xlabel("x", fontsize=12, fontname='serif')
038           plot.set_ylabel("y", fontsize=12, fontname='serif')
039           plot.tick_params(axis='both', length=10, which='major')
040           plot.tick_params(axis='both', length=5,  which='minor')
041           plot.set_xlim(xmin, xmax)
042           plot.set_ylim([-1.2,1.2])
043           plot.minorticks_on()
044           plot.plot(xp, self.f(xp, *self.popt), 'b-')
045           plot.plot(x, y, 'ro', markersize=10)
046           fig.tight_layout()
047           plt.show()
048           fig.savefig('result.pdf', orientation='portrait', \
049                            transparent=False, bbox_inches=None, frameon=None)
050           fig.clf()
051   
052   
053   # -------------------------------------------------------
054   # データをフィッティングする関数
055   # -------------------------------------------------------
056   def fit_func(x, a, b, c):
057       return a*np.sin(x)*np.exp(-b*x)+c
058       
059   
060   # -------------------------------------------------------
061   # メイン関数
062   # -------------------------------------------------------
063   if __name__ == "__main__":
064   
065       x = np.linspace(0, 6*np.pi, 32)
066       noise = 0.05*np.random.normal(size=x.size)
067       y = np.sin(x)*np.exp(-x/5) + noise
068   
069       fit = FITTING(fit_func, x, y)
070       fit.do_fitting()
071       fit.plot(Nx=257)
072       

データをフィッティング(当てはめ)する関数は,056 – 057 行で定義されています.fit_func(x, a, b, c) がデータをフィッティングする関数です.第一引数のxが関数の独立変数です.残りの (a, b, c) がフィッティングパラメーターです.これらの3つの値を調整し,関数とデータの距離を最小にします.

クラス FITTING には,フィッティングを行い結果を表示する機能 (関数) があります.コンストラクター __init__ はフィットする関数とデータをインスタンス変数に,do_fitting フィッティングを実行する,plot は結果を表示する関数です.コンストラクター __init__ は説明するまでも無いでしょう.結果を表示する関数 plot の動作については,Matplotlib を参照ください.

本ページのテーマであるフィッティングを行う関数は,do_fitting です.関数(メソッド) curve_fit(self.f, self.x, self.y) です.引数 self.x がフィッティングする関数,引数 (self.x, self.y) はデータです.戻り値の self.popt はパラメータの最適値です.これは配列で,フィットする関数 fit_func(x, a, b, c) のフィッティングパラメータと\((a,\,b,\,c)\)=(self.popt[0], self.popt[1], self.popt[2]) の関係があります.戻り値のself.pcov は共分散です.これについては,あとで説明します.

プログラムを実行すると,端末に popt と pconv の値が以下のように表示されます.popt は推定される値 \((1,\,0.2,\,0)\) に近いことが分かるでしょう.

popt =  [ 0.93104142  0.18006138  0.00347886] 

pcov =  [[  3.66792254e-03   7.50297351e-04  -1.80400603e-04]
 [  7.50297351e-04   2.65394736e-04  -3.83320805e-05]
 [ -1.80400603e-04  -3.83320805e-05   1.12286047e-04]]

合わせて,図1 に示すプロットも表示されます.

図1: 赤丸は元データ,実線はフィッティングの結果.

共分散とは

scipy.optimize.curve_fit は,あまり聞き慣れない共分散が計算結果として出力されます.

scipy.optimize.curve_fit

フィッティングのモジュール「scipy.optimize.curve_fit」の引数と戻り値,エラー/ワーニングについて説明します.

scipy.optimize.curve_fit(
    f, xdata, ydata, p0=None, sigma=None, absolute_sigma=False,
    check_finite=True, bounds=(-inf, inf), method=None, jac=None, **kwargs)
引数
fcallable
モデル関数 f(x, ...).第1引数は独立変数,第2引数以降はフィティングパラメーターです.
xdataAn M-length sequence or an (k,M)-shaped array for functions with k predictors
独立変数.シーケンスデータです.
ydataM-length sequence
従属変数.ydata = f(xdata, ...) となるシーケンスデータ.
p0None, scalar, or N-length sequence, optionalp0=None
パラメータの初期推定値.None の場合,初期値はすべて 1 となります(introspection を使用して関数のパラメータの数を決定できる場合は [私には理解できない],ValueError が発生します).
sigmaNone or M-length sequence or MxM array, optionalsigma=None
ydataの不確かさの設定です.残差を r = ydata - f(xdata, *popt). と定義すると,sigma の解釈は次元に依存します.
  • 1次元の sigma は,y データのエラーの標準偏差です.この場合,最適化関数 (フィッティングにより最小値にする) は chisq = sum((r/sigma)** 2) です.
  • 2次元の sigma は,y データの誤差の共分散行列です,最適化関数は chisq = r.T @ inv(sigma) @ r です.
None(デフォルト)は1で埋められた1-dシグマと同等です.
absolute_sigmabool, optionalabsolute_sigma=False
Trueの場合,sigma は絶対的な意味で使用され,推定されたパラメータの共分散 pcov はこれらの絶対値を反映します.Falseの場合,sigma は相対値で評価されます. 返されるパラメータ共分散行列 pcov は,スケーリングされた sigma を基にします. 最適なパラメータ popt の chisq が 1 に等しくなるように,スケーリング係数は設定される.すなわち,フィット後の残差の標本分散に一致するように,sigma はスケーリングされる. 数学的に,pcov(absolute_sigma=False) = pcov(absolute_sigma=True) * chisq(popt)/(M-N) です.
check_finitebool, optionalcheck_finite=Truecheck_finite=True
True の場合には,入力配列に inf の nan が含まれているとValueErrorを送出します. False の場合には,入力配列に nans が含まれていると,無意味な結果が黙って生成されることがあります.
bounds2-tuple of array_like, optionalbounds=(-inf, inf)
フィッティングパラメーターの下限と上限を指定します.デフォルトは無限です.タプルの各要素は,パラメータの数に等しい長さの配列か,スカラ(すべてのパラメータで同じになるようにバインドされている)でなければなりません.すべてまたは一部のパラメータの境界を,np.infに適切な符号を付けて無効にすることもできます.
method{'lm', 'trf', 'dogbox'}, optional
最適化のためにアルゴリズムの指定.以下の指定が可能です(筆者は詳細がわからないので,「scipy.optimize.least_squares」を翻訳).
  • 'trf': Trust Region Reflective アルゴリズム.特に,境界のある大きなスパース問題に適しています.般的に頑強な方法.
  • 'dogbox': 矩形の信頼領域を持つ dogleg アルゴリズムで,典型的な使用例は境界の小さな問題です.ランクが不足しているヤコビ行列の問題にはお勧めできません.
  • 'lm': MINPACK で実装されている Levenberg-Marquardt アルゴリズム.境界線やスパースヤコビアンを扱いません.通常,小さな拘束されていない問題の最も効率的な方法です.
拘束されていない問題の場合,デフォルトは 'lm'で,境界が指定されている場合は 'trf'です. 観測の数が変数の数よりも少ない場合,メソッド 'lm'は機能しません.この場合, 'trf'または 'dogbox'を使用してください.
jaccallable, string or None, optionaljac=None
モデル jac(x, ...) の署名付き関数.パラメータに対するモデル関数のヤコビ行列を密な配列_構造として計算します.提供された sigma に従ってスケーリングされます.None(デフォルト)の場合,ヤコビアンは数値で推定されます.'trf' と 'dogbox' メソッドの文字列キーワードを使用して,有限差分スキームを選択できます.(申し訳ない,筆者には理解できない)
kwargs
method='lm'または 'least_squares' の場合は leastsq に渡されるキーワード引数.
戻り値
poptarray
二乗残差の和が最小になるようにパラメータの最適値.フィッテング結果です.
pcov2d array
推定されたpoptの共分散. 対角線は,パラメータ推定値の分散を提供する. パラメータの標準偏差エラーを計算するには,perr = np.sqrt(np.diag(pcov))を使用します.上記のように,シグマパラメータが推定共分散にどのように影響するかは,absolute_sigma引数に依存します.ソリューションのヤコビ行列がフルランクを持たない場合,lmメソッドはnp.infで埋められた行列を返し,他方ではtrfとdogboxメソッドはMoore-Penrose pseudoinverseを使用して共分散を計算します マトリックス.
エラー/ワーニング

サンプルプログラム

フィッティング結果をプロット中に

データを関数でフィッティングするプログラム(fit.py)

001   '''
002   測定データとそのフィッティング結果を表示する.
003   '''
004   import matplotlib.pyplot as plt
005   from matplotlib.backends.backend_pdf import PdfPages
006   import numpy as np
007   from scipy import interpolate
008   import scipy.optimize as optimize
009   
010   
011   # =================================================================
012   # クラス: データを読み込んでフィッティングパラメーターを決める.
013   # =================================================================
014   class Fit_Fuction():
015       '''
016       データからフィッティング関数を決めるクラス.
017       '''
018       def __init__(self, data_file, x_column=0, y_column=1):
019           '''初期化を実行する.'''
020           print('data file: {0:s}'.format(data_file))
021           self.data = np.genfromtxt(data_file)
022           self.x, self.y = self.data[:,x_column], self.data[:,y_column]
023           print(self.x)
024           print('\n\n')
025           print(self.y)
026           # フィッティングの実行.sol[0]がパラメーターリスト
027           self.sol = optimize.curve_fit(self.f, self.x, self.y)
028   
029       def f(self, x, a, s, x0):
030           '''関数の定義: 関数の値を返す.'''
031           return a*np.exp(-(x-x0)**2/(2*s**2))
032   
033       def fitted_f(self, x):
034           '''フィッティング結果関数の値を返す.'''
035           return self.f(x, *self.sol[0])
036   
037       def get_data(self):
038           '''フィッティングの元データを返す.'''
039           return self.x, self.y
040   
041       def get_fit_results(self):
042           '''フィッティングのパラメーターを返す.'''
043           return self.sol    # sol[0] がフィッティング結果
044   
045   
046   # =================================================================
047   # クラス: データをフィッティング関数をプロット
048   # =================================================================
049   class Plot_Data():
050       '''
051       データとフィッティング関数をプロットするクラス.
052       '''
053       def __init__(self, fit_class):
054           '''プロットのクラスの初期化をおこなう.'''
055           self.x, self.y = fit_class.get_data()
056           self.xmin, self.xmax = 0, 1.0
057           self.ymin, self.ymax = 0, 1.2
058           self.fit_x = np.linspace(self.xmin, self.xmax, 512)
059           self.fit_y = fit_class.fitted_f(self.fit_x)
060           title_fn = '$f(x) = a_0\exp[-x^2/(2\sigma^2)]$\n'
061           title_res = '  $a_0={0:g}$\n  $\sigma={1:g}$\n  $x_0={2:g}$'.\
062                           format(*fit_class.get_fit_results()[0])
063           self.title = title_fn + title_res   # プロット中のタイトル
064           self.xl, self.yl = 0.05, 0.8        # タイトルの座標
065   
066       def mk_plot(self):
067           '''元データとフィッティング結果をプロットする'''
068           pp = PdfPages('plot_data.pdf')
069           fig = plt.figure()
070           ax1 = fig.add_subplot(2,2,1)
071           ax1.text(self.xl, self.yl, self.title, color='red', fontsize=10)
072           ax1.set_xlabel("x", fontsize=12, fontname='serif')
073           ax1.set_ylabel("y", fontsize=12, fontname='serif')
074           ax1.tick_params(direction='in', axis='both', length=10, which='major',
075                           bottom=True, top=True, left=True, right=True)
076           ax1.tick_params(direction='in', axis='both', length=5, which='minor',
077                           bottom=True, top=True, left=True, right=True)
078           ax1.set_xlim([self.xmin, self.xmax])
079           ax1.set_ylim([self.ymin, self.ymax])
080           ax1.minorticks_on()
081           ax1.plot(self.fit_x, self.fit_y, color='red', linestyle='solid', linewidth=1.0)
082           ax1.scatter(self.x, self.y, color='blue', marker='o', s=20)
083   
084           plt.tight_layout()
085           plt.subplots_adjust(top=0.9)
086           plt.draw()
087           plt.show()
088           fig.savefig(pp, format='pdf', orientation='portrait', \
089                       transparent=False, bbox_inches=None, frameon=None)
090           fig.clf()
091           pp.close()
092   
093   
094   # =================================================================
095   # メインルーチン
096   # =================================================================
097   if __name__ == '__main__':
098       fit_func = Fit_Fuction('sample.dat', x_column=0, y_column=1)
099       plot = Plot_Data(fit_func)
100       plot.mk_plot()

ページ作成情報

参考資料

  1. マニュアル「scipy.optimize.curve_fit — SciPy v0.19.1 Reference Guide」に詳細が書かれています.

更新履歴

2017年08月10日 新規作成


no counter