如何使用 scipy.optimize.curve_fit 在 python 上拟合一个好的 Lorentzian？答案

【问题标题】：How can I fit a good Lorentzian on python using scipy.optimize.curve_fit?如何使用 scipy.optimize.curve_fit 在 python 上拟合一个好的 Lorentzian？
【发布时间】：2019-03-01 08:32:57
【问题描述】：

我正在尝试拟合具有多个吸收峰（穆斯堡尔光谱）的洛伦兹函数，但 curve_fit 函数无法正常工作，仅拟合几个峰。怎么装？

Figure: Trying to adjusting multi-Lorentzian

下面我展示了我的代码。请帮帮我。

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

def mymodel_hema(x,a1,b1,c1,a2,b2,c2,a3,b3,c3,a4,b4,c4,a5,b5,c5,a6,b6,c6):
    f =  160000 - (c1*a1)/(c1+(x-b1)**2) - (c2*a2)/(c2+(x-b2)**2) - (c3*a3)/(c3+(x-b3)**2) - (c4*a4)/(c4+(x-b4)**2) - (c5*a5)/(c5+(x-b5)**2) - (c6*a6)/(c6+(x-b6)**2)
    return f

def main():
    abre = np.loadtxt('HEMAT_1.dat')
    x = np.zeros(len(abre))
    y = np.zeros(len(abre))

    for i in range(len(abre)):
       x[i] = abre[i,0]
       y[i] = abre[i,1]

    popt,pcov = curve_fit(mymodel_hema, x, y,maxfev=1000000000)

我的数据 --> https://drive.google.com/file/d/1LvCKNdv0oBza_TDwuyNwd29PgQv22VPA/view?usp=sharing

【问题讨论】：

我没有你的数据，但我确实有一个将双洛伦兹峰方程拟合到碳纳米管拉曼光谱的示例：bitbucket.org/zunzuncode/RamanSpectroscopyFit - 该示例使用 scipy 的差分进化遗传算法模块，该模块使用拉丁超立方体算法来确保对参数空间进行彻底搜索，这需要在其中搜索的参数范围 - 在示例代码中，这些范围是根据数据的最大值和最小值确定的。
请提供一些我们可以轻松复制粘贴的数据。
您好，我用我的数据编辑了我的问题。只需下载即可。
您是否也尝试过为参数提供初步猜测？这通常很有帮助。

标签： python curve-fitting

【解决方案1】：

此代码使用leastsq 而不是curve_fit，因为后者需要固定数量的参数。在这里我不想要这个，因为我让代码“决定”那里有多少个峰值。请注意，我对数据进行了缩放以简化拟合。真正的拟合参数可以很容易地按比例缩小（和标准误差传播）

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

def lorentzian( x, x0, a, gam ):
    return a * gam**2 / ( gam**2 + ( x - x0 )**2)

def multi_lorentz( x, params ):
    off = params[0]
    paramsRest = params[1:]
    assert not ( len( paramsRest ) % 3 )
    return off + sum( [ lorentzian( x, *paramsRest[ i : i+3 ] ) for i in range( 0, len( paramsRest ), 3 ) ] )

def res_multi_lorentz( params, xData, yData ):
    diff = [ multi_lorentz( x, params ) - y for x, y in zip( xData, yData ) ]
    return diff

xData, yData = np.loadtxt('HEMAT_1.dat', unpack=True )
yData = yData / max(yData)

generalWidth = 1

yDataLoc = yData
startValues = [ max( yData ) ]
counter = 0

while max( yDataLoc ) - min( yDataLoc ) > .1:
    counter += 1
    if counter > 20: ### max 20 peak...emergency break to avoid infinite loop
        break
    minP = np.argmin( yDataLoc )
    minY = yData[ minP ]
    x0 = xData[ minP ]
    startValues += [ x0, minY - max( yDataLoc ), generalWidth ]
    popt, ier = leastsq( res_multi_lorentz, startValues, args=( xData, yData ) )
    yDataLoc = [ y - multi_lorentz( x, popt ) for x,y in zip( xData, yData ) ]

print popt
testData = [ multi_lorentz(x, popt ) for x in xData ]

fig = plt.figure()
ax = fig.add_subplot( 1, 1, 1 )
ax.plot( xData, yData )
ax.plot( xData, testData )
plt.show()

提供

[ 9.96855817e-01  4.94106598e+02 -2.82103813e-01  4.66272773e+00
  2.80688160e+01 -2.72449246e-01  4.71728295e+00  1.31577189e+02
 -2.29698620e-01  4.20685229e+00  4.01421993e+02 -1.85917255e-01
  5.57859380e+00  2.29704607e+02 -1.47193792e-01  3.91112196e+00
  3.03387957e+02 -1.37127711e-01  4.39571905e+00]

和

【讨论】：