0

このトピックに関する質問を既に開いていましたが、そこに投稿する必要があるかどうかわからなかったので、ここに新しい質問を開きました。

2 つ以上のピークをフィッティングするときに、また問題が発生します。最初の問題は、計算されたサンプル関数で発生します。

xg = np.random.uniform(0,1000,500)
mu1 = 200
sigma1 = 20
I1 = -2

mu2 = 800
sigma2 = 20
I2 = -1

yg3 = 0.0001*xg
yg1 = (I1 / (sigma1 * np.sqrt(2 * np.pi))) * np.exp( - (xg - mu1)**2 / (2 * sigma1**2) )
yg2 = (I2 / (sigma2 * np.sqrt(2 * np.pi))) * np.exp( - (xg - mu2)**2 / (2 * sigma2**2) )
yg=yg1+yg2+yg3

plt.figure(0, figsize=(8,8))
plt.plot(xg, yg, 'r.')

以下に示すドキュメントで見つけた2つの異なるアプローチを試しました（私のデータ用に変更されています）が、どちらも間違ったフィッティングデータと乱雑なグラフの混乱をもたらします（フィッティングステップごとに1行だと思います）。

1 回目の試行:

import numpy as np
from lmfit.models import PseudoVoigtModel, LinearModel, GaussianModel, LorentzianModel
import sys
import matplotlib.pyplot as plt

gauss1 = PseudoVoigtModel(prefix='g1_')
pars.update(gauss1.make_params())

pars['g1_center'].set(200)
pars['g1_sigma'].set(15, min=3)
pars['g1_amplitude'].set(-0.5)
pars['g1_fwhm'].set(20, vary=True)
#pars['g1_fraction'].set(0, vary=True)

gauss2  = PseudoVoigtModel(prefix='g2_')
pars.update(gauss2.make_params())

pars['g2_center'].set(800)
pars['g2_sigma'].set(15)
pars['g2_amplitude'].set(-0.4)
pars['g2_fwhm'].set(20, vary=True)
#pars['g2_fraction'].set(0, vary=True)

mod = gauss1 + gauss2 + LinearModel()

pars.add('intercept', value=0, vary=True)
pars.add('slope', value=0.0001, vary=True)

init = mod.eval(pars, x=xg)
out = mod.fit(yg, pars, x=xg)

print(out.fit_report(min_correl=0.5))
plt.figure(5, figsize=(8,8))
out.plot_fit()

「fraction」パラメータを含めると、よく

'NameError: name 'pv1_fraction' is not defined in expr='<_ast.Module object at 0x00000000165E03C8>'.

定義する必要がありますが。このアプローチでも、実際のデータに対してこのエラーが発生します。

2 回目の試行:

import pandas as pd
import matplotlib.pyplot as plt
import numpy as np
import lmfit

def gauss(x, sigma, mu, A):
    return A*np.exp(-(x-mu)**2/(2*sigma**2))

def linear(x, m, n):
    return m*x + n

peak1 = lmfit.model.Model(gauss, prefix='p1_')
peak2 = lmfit.model.Model(gauss, prefix='p2_')
lin = lmfit.model.Model(linear, prefix='l_')
model = peak1 + lin + peak2
params = model.make_params()

params['p1_mu'] = lmfit.Parameter(value=200, min=100, max=250)
params['p2_mu'] = lmfit.Parameter(value=800, min=100, max=1000)
params['p1_sigma'] = lmfit.Parameter(value=15, min=0.01)
params['p2_sigma'] = lmfit.Parameter(value=20, min=0.01)
params['p1_A'] = lmfit.Parameter(value=-2, min=-3)
params['p2_A'] = lmfit.Parameter(value=-2, min=-3)
params['l_m'] = lmfit.Parameter(value=0)
params['l_n'] = lmfit.Parameter(value=0)

out = model.fit(yg, params, x=xg)

print out.fit_report()
plt.figure(8, figsize=(8,8))
out.plot_fit()

どちらの場合も、結果は次のようになります。すべてのフィッティング試行をプロットしているように見えますが、正しく解決することはありません。最適なパラメーターは、私が指定した範囲内にあります。

ここに画像の説明を入力

ここに画像の説明を入力

このタイプのエラーを知っている人はいますか？または、これに対する解決策はありますか？そして、これらのアプローチNameErrorからモデル関数を呼び出すときに回避する方法を知っている人はいますか?lmfit

4

2 に答える 2

2

私はあなたのために多少許容できる解決策を持っています。あなたのデータがどのように変化するか分からないので、一般的な意味でうまくいくとは言えませんが、あなたが始めるべきです. データが 0 ～ 1000 に沿っており、示したように線に沿って 2 つのピークまたはディップがある場合、それは機能するはずです。

scipy curve_fit を使用して、関数のすべてのコンポーネントを 1 つの関数にまとめました。開始位置をcurve_fit関数に渡すことができます。(おそらく、使用しているライブラリでこれを行うことができますが、私はそれに慣れていません) mu パラメーターを変化させて、二乗誤差が最小のものを見つけるループ内ループがあります。何度もデータを適合させる必要がある場合、またはリアルタイムのシナリオでデータを適合させる必要がある場合、これは適していませんが、一部のデータを適合させる必要がある場合は、このコードを起動してコーヒーを飲みます。

from scipy.optimize import curve_fit
import numpy as np
import matplotlib.pyplot as plt
import pylab
from matplotlib import cm as cm
import time

def my_function_big(x, m, n, #lin vars
                       sigma1, mu1, I1,  #gaussian 1
                       sigma2, mu2, I2):  #gaussian 2

  y = m * x + n + (I1 / (sigma1 * np.sqrt(2 * np.pi))) * np.exp( - (x - mu1)**2 / (2 * sigma1**2) ) + (I2 / (sigma2 * np.sqrt(2 * np.pi))) * np.exp( - (x - mu2)**2 / (2 * sigma2**2) )
  return y


#make some data
xs = np.random.uniform(0,1000,500)
mu1 = 200
sigma1 = 20
I1 = -2

mu2 = 800
sigma2 = 20
I2 = -1

yg3 = 0.0001 * xs
yg1 = (I1 / (sigma1 * np.sqrt(2 * np.pi))) * np.exp( - (xs - mu1)**2 / (2 * sigma1**2) )
yg2 = (I2 / (sigma2 * np.sqrt(2 * np.pi))) * np.exp( - (xs - mu2)**2 / (2 * sigma2**2) )
ys = yg1 + yg2 + yg3

xs = np.array(xs)
ys = np.array(ys)

#done making data



#start a double loop...very expensive but this is quick and dirty
#it would seem that the regular optimizer has trouble finding the minima so i 
#found that having the near proper mu values helped it zero in much better

start = time.time()

serr = []
_x = []
_y = []
for x in np.linspace(0, 1000, 61):
  for y in np.linspace(0, 1000, 61):
    cfiti = curve_fit(my_function_big, xs, ys, p0=[0, 0, 1, x, 1, 1, y, 1], maxfev=20000000)
    serr.append(np.sum((my_function_big(xs, *cfiti[0]) - ys) ** 2))
    _x.append(x)
    _y.append(y)

serr = np.array(serr)
_x = np.array(_x)
_y = np.array(_y)

print 'done loop in loop fitting'
print 'time: %0.1f' % (time.time() - start)

gridsize=20
plt.subplot(111)
plt.hexbin(_x, _y, C=serr, gridsize=gridsize, cmap=cm.jet, bins=None)
plt.axis([_x.min(), _x.max(), _y.min(), _y.max()])
cb = plt.colorbar()
cb.set_label('SE')
plt.show()   

ix = np.argmin(serr.ravel())
mustart1 = _x.ravel()[ix]
mustart2 = _y.ravel()[ix]
print mustart1
print mustart2

cfit = curve_fit(my_function_big, xs, ys, p0=[0, 0, 1, mustart1, 1, 1, mustart2, 1], maxfev=2000000000)

xp = np.linspace(0, 1000, 1001)

plt.figure()
plt.scatter(xs, ys) #plot synthetic dat
plt.plot(xp, my_function_big(xp, *cfit[0]), '-', label='fit function')  #plot data evaluated along 0-1000
plt.legend(loc=3, numpoints=1, prop={'size':12})
plt.show()    
pylab.close()

二乗誤差がミュー空間でどのように振る舞うかを視覚化します。

フィット関数

幸運を！

于 2015-12-23T20:21:55.507 に答える