洛伦兹拟合的半参数全宽与python中的scipy curve_fit

Question

我想确定数据的洛伦兹拟合的半参数全宽 (FWHM)，我正在使用 scipy 中的

curve_fit

函数。当我使用高斯拟合时，FWHM 通过

2.35482*sigma

计算。一般来说，在这种情况下，洛伦兹（参见 MWE）的 FWHM 必须确定为

2*sigma

，但这并不适合通过将其显示到绘图中（绘图显示

1*sigma

）。我做错了什么？

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


x = np.linspace(0,5,10)
y = [0.0,0.05,0.16,0.3,0.6,0.8,0.6,0.3,0.16,0.0]


def fit_Lorenz(x, y, f0, init):
    
    def Lorentzian(f, amp, cen, wid, Offset):
        return np.sqrt(amp/((f - cen)**2 + wid**2 / 4 )) + Offset
     
    init_vals = init
    popt, cov = curve_fit(Lorentzian, x, y ,p0=init_vals, maxfev= 10000)
    
    Amp = popt[0]
    cen = np.abs(popt[1])
    wid = np.abs(popt[2])
    Offset = popt[3]
    
    x_fit = np.linspace(min(x), max(x), 1000)   
    y_fit = np.zeros(len(x_fit))
    
    for i in range(len(x_fit)):
        y_fit[i] = Lorentzian(x_fit[i],Amp,cen,wid,Offset)
    
    return x_fit, y_fit, wid, cen, Amp

init_vals= [1,2,1,1]
x_fit, y_fit, sigma, x0, Amp = fit_Lorenz(x,y, max(y), init_vals)

FWHM = sigma
FWHM_val = (max(y_fit) - min(y_fit))/2 + min(y_fit)


plt.plot(x, y, '.', label='data', color='blue')
plt.plot(x_fit, y_fit, '-', label='Lorentzian fit', color='red')
plt.xlim(0,5)
f1 = x0-FWHM/2
f2 = x0+FWHM/2
plt.hlines(FWHM_val,f1,f2,linestyle='solid', color='red')
plt.legend(loc='upper left', fontsize=13)

Answer 1

我快速搜索了一下，发现了洛伦兹方程的一个不同的方程。使用 this site，我将您的方程更改为其形式，将半高计算为

0.5*Amp/sigma

（您将值称为

sigma

，但该网站使用 gamma），同时考虑偏移量，并绘制结果。

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

x = np.linspace(0, 5, 10)
y = [0.0, 0.05, 0.16, 0.3, 0.6, 0.8, 0.6, 0.3, 0.16, 0.0]


def fit_Lorenz(x, y, f0, init):

    def Lorentzian(f, amp, cen, wid, Offset):
        return amp*wid/((f-cen)**2 + wid**2) + Offset

    init_vals = init
    popt, cov = curve_fit(Lorentzian, x, y, p0=init_vals, maxfev=10000)

    Amp = popt[0]
    cen = np.abs(popt[1])
    wid = np.abs(popt[2])
    Offset = popt[3]

    x_fit = np.linspace(min(x), max(x), 1000)
    y_fit = np.zeros(len(x_fit))

    for i in range(len(x_fit)):
        y_fit[i] = Lorentzian(x_fit[i], Amp, cen, wid, Offset)

    return x_fit, y_fit, wid, cen, Amp, Offset


init_vals = [1, 2, 1, 1]
x_fit, y_fit, sigma, x0, Amp, Offset = fit_Lorenz(x, y, max(y), init_vals)

# peak is Amp/sigma high, so half the height is half that
half_height = 0.5*Amp/sigma + Offset
FWHM = 2*sigma
f1 = x0-FWHM/2
f2 = x0+FWHM/2


plt.plot(x, y, '.', label='data', color='blue')
plt.plot(x_fit, y_fit, '-', label='Lorentzian fit', color='red')
plt.xlim(0, 5)
plt.hlines(half_height, f1, f2, linestyle='solid', color='green')
plt.legend(loc='upper left', fontsize=13)

Answer 2

我建议

lmfit

（https://lmfit.github.io/lmfit-py/，

pip install lmfit

，披露我是主要作者）使这变得更容易一些，因为大部分工作都需要是内置的。使用 lmfit，它看起来像（不是缺少缩进）：

import numpy as np
import matplotlib.pyplot as plt
from lmfit.models import LorentzianModel, ConstantModel

x = np.linspace(0, 5, 10)
y = np.array([0.0, 0.05, 0.16, 0.3, 0.6, 0.8, 0.6, 0.3, 0.16, 0.0])

lmodel = LorentzianModel() + ConstantModel() 
params = lmodel.make_params(amplitude=1, center=2, sigma=1, c=0)
result = lmodel.fit(y, params, x=x)

print(result.fit_report())

parvals = result.params.valuesdict()
height = parvals['height']
offset = parvals['c']
center = parvals['center']
hwhm   = parvals['fwhm']/2.0

xfine = np.linspace(x.min(), x.max(), 500)
plt.plot(x,  y, 'o', label='data')
plt.plot(x, result.best_fit, '+', label='fit at data points')
plt.hlines(offset+height/2, center-hwhm, center+hwhm, color='red')
plt.plot(xfine, result.eval(x=xfine), label='smoothed fit')

plt.legend()
plt.show()

打印出报告

[[[Model]]
    (Model(lorentzian) + Model(constant))
[[Fit Statistics]]
    # fitting method   = leastsq
    # function evals   = 48
    # data points      = 10
    # variables        = 4
    chi-square         = 0.00325547
    reduced chi-square = 5.4258e-04
    Akaike info crit   = -72.3000261
    Bayesian info crit = -71.0896857
    R-squared          = 0.99542520
[[Variables]]
    amplitude:  3.05607907 +/- 0.23483180 (7.68%) (init = 1)
    center:     2.77164969 +/- 0.02190759 (0.79%) (init = 2)
    sigma:      1.05036602 +/- 0.06440647 (6.13%) (init = 1)
    c:         -0.12587146 +/- 0.02700091 (21.45%) (init = 0)
    fwhm:       2.10073203 +/- 0.12881294 (6.13%) == '2.0000000*sigma'
    height:     0.92613452 +/- 0.02737682 (2.96%) == '0.3183099*amplitude/max(1e-15, sigma)'
[[Correlations]] (unreported correlations are < 0.100)
    C(amplitude, c)     = -0.9573
    C(amplitude, sigma) = +0.9328
    C(sigma, c)         = -0.8587][1]][1]

并制作出

的图形