物理の駅 Physics station by 現役研究者

テクノロジーは共有されてこそ栄える

Python lmfitを使って、マルチガウシアンでフィッティングする

phst.hateblo.jp

の続編

mzfit が使っている tensorflow が重すぎなのと、モジュール読み込みに時間がかかるのと、機能的に不足する点があったのと、マルチプラットフォームであるべきという思想に反するので、lmfitを使ってマルチガウシアンでのフィッティングを書いてみた。

lmfitの公式ドキュメントのFitting Multiple Peaks用コードを参考にした。

データ生成部分。

import numpy as np

np.random.seed(seed=0)
signal1 = np.random.normal(5, 1, 2000)
signal2 = np.random.normal(10, 2, 10000)
signal3 = np.random.normal(15, 0.5, 5000)
data = np.concatenate([signal1, signal2, signal3])

hist, bin_edges = np.histogram(data,bins=100,range=(2,18))
bin_center = (bin_edges[1:]+bin_edges[:-1])*0.5
arr_x = []
arr_y = []
arr_yerror = []
for x,y in zip(bin_center,hist):
    if y==0:continue
    arr_x.append(float(x))
    arr_y.append(float(y))
    arr_yerror.append(y**0.5)
arr_x = np.array(arr_x)
arr_y = np.array(arr_y)
arr_yerror = np.array(arr_yerror)

3つのガウシアンモデルで関数形を用意し、初期値を与える。初期値を全く与えなくてもfittingはうまくいったが、作法として初期値は与えておくべきである。set関数には、valueだけでなく、パラメータの範囲はminmaxとして引数を与えることができる。

from lmfit.models import ExponentialModel, GaussianModel

gauss1 = GaussianModel(prefix='g1_')
pars = gauss1.make_params()
gauss2 = GaussianModel(prefix='g2_')
pars.update(gauss2.make_params())
gauss3 = GaussianModel(prefix='g3_')
pars.update(gauss3.make_params())

pars['g1_center'].set(value=3)
pars['g1_amplitude'].set(value=100)
pars['g2_center'].set(value=10)
pars['g2_sigma'].set(value=1)
pars['g2_amplitude'].set(value=100)
pars['g3_center'].set(value=16)
pars['g3_sigma'].set(value=1)
pars['g3_amplitude'].set(value=100)

mod = gauss1 + gauss2 + gauss3

evalで初期パラメータを使ったarr_yの値を取得する。fitでフィッティングを行う。先の記事のとおり、分散共分散行列を正しく求めるため、引数にweights=1/arr_yerrorscale_covar=False を与える。

init = mod.eval(pars, x=arr_x)
result = mod.fit(arr_y, pars, x=arr_x, weights=1/arr_yerror, scale_covar=False)

result.result.params

グラフ化

import matplotlib.pyplot as plt

fig, axes = plt.subplots(1, 2, figsize=(12.8, 4.8))
axes[0].errorbar(arr_x, arr_y, yerr=arr_yerror, marker='.', linestyle="")
axes[0].plot(arr_x, init, '--', label='initial fit')
axes[0].plot(arr_x, result.best_fit, '-', label='best fit')
axes[0].legend()

comps = result.eval_components(x=arr_x)
axes[1].errorbar(arr_x, arr_y, yerr=arr_yerror, marker='.', linestyle="")
axes[1].plot(arr_x, comps['g1_'], '--', label='Gaussian component 1')
axes[1].plot(arr_x, comps['g2_'], '-', label='Gaussian component 2')
axes[1].plot(arr_x, comps['g3_'], ':', label='Gaussian component 3')
axes[1].legend()

plt.show()

マルチガウシアンでフィット