洛伦兹拟合两种编写代码的方式
Lorentzian fit two ways of writing a code
我现在正在为洛伦兹曲线拟合而苦苦挣扎。我会尽力解释我的问题。我需要为洛伦兹曲线拟合编写自己的代码,这样我就可以在方程式中添加一些东西。我用 model
和 def
实现了 Lorentzian 拟合,我写了类似的,但它不起作用。查看我的代码:
这是我的数据:
for dataset in [Bxfft]:
dataset = np.asarray(dataset)
freqs, psd = signal.welch(dataset, fs=266336/300, window='hamming', nperseg=16192, scaling='density')
plt.semilogy(freqs[30:-7000], psd[30:-7000]/dataset.size**0, color='r', label='Bx')
x = freqs[100:-7900]
y = psd[100:-7900]
这里是我定义的洛伦兹曲线拟合:
def lorentzian(x, amp, cen, sig):
return (amp/np.pi) * (sig/(x-cen)**2 + sig**2)
model = Model(lorentzian)
pars = model.make_params(amp=6, cen=5, sig=1)
pars['amp'].max = 6
result = model.fit(y, pars, x=x)
final_fit = result.best_fit
print(result.fit_report(min_correl=0.25))
plt.plot(x, final_fit, 'k--', linewidth=3)
此处由模型函数完成:
model2 = LorentzianModel()
params2 = model2.make_params(amplitude=6, center=5, sigma=1)
params2['amplitude'].value = 6
result2 = model2.fit(y, params2, x=x)
final_fit2 = result2.best_fit
print(result2.fit_report(min_correl=0.25))
plt.plot(x, final_fit2, 'k--', linewidth=3)
上图为 def
洛伦兹分布,下图为 model
洛伦兹分布。
结果是:
看起来像是括号问题。这个:
(amp/np.pi) * (sig/(x-cen)**2 + sig**2)
不是洛伦兹。这个:
(amp/np.pi) * (sig/((x-cen)**2 + sig**2))
是。此外,在罕见的情况下,您可能会遇到轻微的整数问题 cen,x,sig
都是整数。您可以使用 math.pow
来解决这个问题,或者他们在 lmfit
中所做的,然后将 x
乘以一个浮点数:1.0*x-cen
.
作为旁注,lmfit
出于某种原因等效地编写了此函数,但 a bit differently(在洛伦兹页面上找到)。不过我看不出这是什么原因。
我现在正在为洛伦兹曲线拟合而苦苦挣扎。我会尽力解释我的问题。我需要为洛伦兹曲线拟合编写自己的代码,这样我就可以在方程式中添加一些东西。我用 model
和 def
实现了 Lorentzian 拟合,我写了类似的,但它不起作用。查看我的代码:
这是我的数据:
for dataset in [Bxfft]:
dataset = np.asarray(dataset)
freqs, psd = signal.welch(dataset, fs=266336/300, window='hamming', nperseg=16192, scaling='density')
plt.semilogy(freqs[30:-7000], psd[30:-7000]/dataset.size**0, color='r', label='Bx')
x = freqs[100:-7900]
y = psd[100:-7900]
这里是我定义的洛伦兹曲线拟合:
def lorentzian(x, amp, cen, sig):
return (amp/np.pi) * (sig/(x-cen)**2 + sig**2)
model = Model(lorentzian)
pars = model.make_params(amp=6, cen=5, sig=1)
pars['amp'].max = 6
result = model.fit(y, pars, x=x)
final_fit = result.best_fit
print(result.fit_report(min_correl=0.25))
plt.plot(x, final_fit, 'k--', linewidth=3)
此处由模型函数完成:
model2 = LorentzianModel()
params2 = model2.make_params(amplitude=6, center=5, sigma=1)
params2['amplitude'].value = 6
result2 = model2.fit(y, params2, x=x)
final_fit2 = result2.best_fit
print(result2.fit_report(min_correl=0.25))
plt.plot(x, final_fit2, 'k--', linewidth=3)
上图为 def
洛伦兹分布,下图为 model
洛伦兹分布。
结果是:
看起来像是括号问题。这个:
(amp/np.pi) * (sig/(x-cen)**2 + sig**2)
不是洛伦兹。这个:
(amp/np.pi) * (sig/((x-cen)**2 + sig**2))
是。此外,在罕见的情况下,您可能会遇到轻微的整数问题 cen,x,sig
都是整数。您可以使用 math.pow
来解决这个问题,或者他们在 lmfit
中所做的,然后将 x
乘以一个浮点数:1.0*x-cen
.
作为旁注,lmfit
出于某种原因等效地编写了此函数,但 a bit differently(在洛伦兹页面上找到)。不过我看不出这是什么原因。