[lmfit] 14. 指数関数的減衰モデルによるカーブフィッティング

はじめに

lmfitは非線形最小二乗法を用いてカーブフィットするためのライブラリであり、Scipy.optimize.curve_fitの拡張版に位置する。ここでは、lmfitでデータを指数関数的減衰モデルによりカーブフィッティングする方法について説明する。

コード

	import matplotlib.pyplot as plt
	import numpy as np
	from numpy.random import default_rng
	from lmfit.models import ExponentialModel
	x = np.linspace(0,5,50)
	def f(x,a,t):
	return a*np.exp(-x/t)
	rng = default_rng()
	noise = rng.random(50)
	y = f(x,10,1) + noise
	fig, ax = plt.subplots(dpi=130)
	ax.plot(x,y,"o")
	ax.set(xlabel="x",ylabel="y")
	plt.savefig("lmfit14_1.png",dpi=130)
	plt.show()
	model = ExponentialModel()
	params = model.guess(y, x=x)
	params
	"""
	name value initial value min max vary
	amplitude 8.25485286 None -inf inf True
	decay 1.51239229 None -inf inf True
	"""
	result = model.fit(y, params, x=x)
	print(result.fit_report())
	"""
	[[Model]]
	Model(exponential)
	[[Fit Statistics]]
	# fitting method = leastsq
	# function evals = 16
	# data points = 50
	# variables = 2
	chi-square = 6.67349486
	reduced chi-square = 0.13903114
	Akaike info crit = -96.6939658
	Bayesian info crit = -92.8699198
	[[Variables]]
	amplitude: 9.87943792 +/- 0.19812223 (2.01%) (init = 8.254853)
	decay: 1.23987140 +/- 0.03684039 (2.97%) (init = 1.512392)
	[[Correlations]] (unreported correlations are < 0.100)
	C(amplitude, decay) = -0.680
	"""
	fig, ax = plt.subplots(dpi=130)
	result.plot_fit(ax=ax)
	plt.savefig("lmfit14_2.png",dpi=130)
	plt.show()
	fig = plt.figure(dpi=130,figsize=(6,6))
	result.plot(fig=fig)
	plt.savefig("lmfit14_3.png",dpi=130)
	plt.show()

view raw lmfit14.py hosted with ❤ by GitHub

解説

モジュールのインポートなど

	import matplotlib.pyplot as plt
	import numpy as np
	from numpy.random import default_rng
	from lmfit.models import ExponentialModel

view raw lmfit14.py hosted with ❤ by GitHub

バージョン

	#version
	import matplotlib
	print(matplotlib.__version__)
	3.3.4
	import lmfit
	print(lmfit.__version__)
	1.0.2
	print(np.__version__)
	1.20.1

view raw lmfit14.py hosted with ❤ by GitHub

データの生成

	x = np.linspace(0,5,50)
	def f(x,a,t):
	return a*np.exp(-x/t)
	rng = default_rng()
	noise = rng.random(50)
	y = f(x,10,1) + noise

view raw lmfit14.py hosted with ❤ by GitHub

yデータはf(x,a,t)で定義した関数で作成する。rng = default_rng()とし、rng.random(50)でランダムデータを作成しyデータのノイズとする。xとyの関係を図で示すと以下のようになる。

	fig, ax = plt.subplots(dpi=130)
	ax.plot(x,y,"o")
	ax.set(xlabel="x",ylabel="y")
	plt.savefig("lmfit14_1.png",dpi=130)
	plt.show()

view raw lmfit14.py hosted with ❤ by GitHub

モデルの定義

model = ExponentialModel()

view raw lmfit14.py hosted with ❤ by GitHub

lmfit.models の ExponentialModelをモデル関数として用いる。

初期パタメータの推定

params = model.guess(y, x=x)

view raw lmfit14.py hosted with ❤ by GitHub

model.guess(y, x=x)により、上図のデータを指数関数的減衰関数で近似するためのフィッティングパラメータについて、初期値を推定する。パラメータ(params)は以下のようになる。

	params
	"""
	name value initial value min max vary
	amplitude 8.25485286 None -inf inf True
	decay 1.51239229 None -inf inf True
	"""

view raw lmfit14.py hosted with ❤ by GitHub

カーブフィット

result = model.fit(y, params, x=x)

view raw lmfit14.py hosted with ❤ by GitHub

model.fit(y, params, x=x)により、カーブフィッティングを実行する。

フィッティング結果の表示

print(result.fit_report())により、フィッティングの結果を見ることができる。

	print(result.fit_report())
	"""
	[[Model]]
	Model(exponential)
	[[Fit Statistics]]
	# fitting method = leastsq
	# function evals = 16
	# data points = 50
	# variables = 2
	chi-square = 6.67349486
	reduced chi-square = 0.13903114
	Akaike info crit = -96.6939658
	Bayesian info crit = -92.8699198
	[[Variables]]
	amplitude: 9.87943792 +/- 0.19812223 (2.01%) (init = 8.254853)
	decay: 1.23987140 +/- 0.03684039 (2.97%) (init = 1.512392)
	[[Correlations]] (unreported correlations are < 0.100)
	C(amplitude, decay) = -0.680
	"""

view raw lmfit14.py hosted with ❤ by GitHub

result.plot_fit()によりデータとフィッティングカーブが表示される。

	fig, ax = plt.subplots(dpi=130)
	result.plot_fit(ax=ax)
	plt.savefig("lmfit14_2.png",dpi=130)
	plt.show()

view raw lmfit14.py hosted with ❤ by GitHub

result.plot()とすることで残差とともにフィッティング結果が表示される。

	fig = plt.figure(dpi=130,figsize=(6,6))
	result.plot(fig=fig)
	plt.savefig("lmfit14_3.png",dpi=130)
	plt.show()

view raw lmfit14.py hosted with ❤ by GitHub

コードをダウンロード(.pyファイル)

コードをダウンロード(.ipynbファイル)

参考

Built-in Fitting Models in the models module — Non-Linear Least-Squares Minimization and Curve-Fitting for Python

lmfit.github.io

Random Generator — NumPy v2.2 Manual

numpy.org