フィッティング
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」があります．これを使うと容易に誤差を持つデータを任意の関数でフィッティングすることができます．これを使うためにのステップは，次のとおりです．
- モジュールをインポートします．よく使われるコマンドは，「from scipy.optimize import curve_fit」です．
- フィッティングする関数を定義します．関数の第一引数が独立変数，第二引数以降はフィッティングパラメーターです．例えば，「\(f(x,\, a\, b)\)」のようにです．
- データを準備します．例えば，\((x,\,y)\) のようにです．\(x\) と \(y\) のいずれも，シーケンスデータで同じ個数の実数から成ります．
- フィッティングを行います．コマンドは「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) がフィッティングパラメーターです．これらの３つの値を調整し，関数とデータの距離を最小にします．
クラス 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 です．
- 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 アルゴリズム．境界線やスパースヤコビアンを扱いません．通常，小さな拘束されていない問題の最も効率的な方法です．
- 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()
ページ作成情報
参考資料
- マニュアル「scipy.optimize.curve_fit — SciPy v0.19.1 Reference Guide」に詳細が書かれています．
更新履歴
| 2017年08月10日 | 新規作成 |