の続編
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
だけでなく、パラメータの範囲はmin
やmax
として引数を与えることができる。
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_yerror
とscale_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()