在 Python 3 中使用 Scipy 将多个洛伦兹拟合到布里渊光谱

Question

我正在尝试使用 scipy.optimize.curve_fit 拟合布里渊光谱（有几个峰）。我有多个具有多个峰的光谱，我正在尝试将它们与洛伦兹函数（每个峰一个洛伦兹）拟合。我正在尝试自动化批量分析的过程（即，使用 scipy 的峰值查找算法来获取峰值位置、峰值宽度和峰值高度，并将它们用作拟合的初始猜测）。我现在正在研究一个光谱，看看总体思路是否有效，然后我会将其扩展为自动并使用我拥有的所有光谱。到目前为止，我已经这样做了：

import numpy as np

import matplotlib.pyplot as plt

from scipy.signal import find_peaks

from scipy.optimize import curve_fit

#import y data from linked google sheet 

y_data = np.loadtxt( 'y_peo.tsv' )

#define x data 

x_data = np.arange( len( y_data ) )

#find peak properties (peak position, amplitude, full width half maximum ) to use as 

#initial guesses for the curve_fit function 

pk, properties = find_peaks(y_data, height = 3, width = 3, prominence=0.1 ) #pk returns peaks position, properties returns 

#other properties associated with the peaks

I = properties ['peak_heights'] #amplitude

fwhm = (properties['widths']) #full width half maximum 

#define sum of lorentzians

def func(x, *params): 

    y = np.zeros_like(x)

    for i in range(0, len(params), 3):

        x_0 = params[i]

        I = params[i+1]

        y_0 = params[i+2]

        y = y + (I*y_0**2)/(y_0**2 +(x-x_0)**2) 

    return y

#initial guess list, to be populated with parameters found above (pk, I, fwhm)

guess = [] 

for i in range(len(pk)): 

    guess.append(pk[i])

    guess.append(I[i])

    guess.append(fwhm[i]) 

#convert list to array

guess=np.array(guess)

#fit

popt, pcov = curve_fit(func, x_data, y_data, p0=guess, method = 'lm',  maxfev=1000000)

print(popt)

fit = func(x_data, *popt)

#plotting

fig= plt.figure(figsize=(10,5))

ax= fig.add_subplot(1,1,1)

ax.plot(pk, y_data[pk], 'o', ms=5)

ax.plot(x_data, y_data, 'ok', ms=1)

ax.plot(x_data, fit , 'r--', lw=0.5)

其中y_peo是因变量（这里是谷歌表格文件中的值：https ://docs.google.com/spreadsheets/d/1UB2lqs0Jmhthed7B0U9rznqcbpEMRmRX99c-iOBHLeE/edit?usp=sharing ）和像素是自变量（任意在代码中创建的值数组）。这就是我得到的：频谱拟合的结果。知道为什么最后一个三胞胎没有按预期安装吗？我检查并通过 scipy.signal.find_peaks 函数正确找到了所有峰值 - 因此我认为问题出在 scipy.optimize.curve_fit 上，因为我必须增加 maxfev 的数量才能“工作”。关于如何以更聪明的方式解决这个问题的任何想法？

Answer 1

通过一些小的修改，OP 的代码运行得很好。现在看起来像：import numpy as npimport matplotlib.pyplot as pltfrom scipy.optimize import leastsqfrom scipy.signal import find_peaksdef 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 diffy0 = np.loadtxt( 'y_peo.tsv' )yData = np.loadtxt( 'y_peo.tsv' )xData = np.arange( len( yData ) )yGround = min( yData )yData = yData - yGroundyAmp = max( yData )yData = yData / yAmp #initial properties of peaks pk, properties = find_peaks( yData, height = .05, width = 3 )#, prominence=0.1 )#extract peak heights and fwhm I = properties [ 'peak_heights' ]fwhm = properties[ 'widths' ]guess = [0]for i in range( len( pk ) ):     guess.append( pk[i] )    guess.append( I[i] )    guess.append( fwhm[i] ) guess=np.array( guess )popt, pcov = leastsq( res_multi_lorentz, x0=guess, args=( xData, yData ) )print( popt )testData = [ multi_lorentz( x, popt ) for x in xData ]fitData = [ yGround + yAmp * multi_lorentz( x, popt ) for x in xData ]fig= plt.figure( figsize=( 10, 5 ) )ax= fig.add_subplot( 2, 1, 1 )bx= fig.add_subplot( 2, 1, 2 )ax.plot( pk, yData[pk], 'o', ms=5 )ax.plot( xData, yData, 'ok', ms=1 )ax.plot( xData, testData , 'r--', lw=0.5 )bx.plot( xData, y0, ls='', marker='o', markersize=1 )bx.plot( xData, fitData )plt.show()给予[9.39456104e-03 6.55864388e+01 5.73522507e-02 5.46727721e+00 1.21329586e+02 2.97187567e-01 2.12738107e+00 1.76823266e+02 1.57244131e-01 4.55424037e+00 3.51556692e+02 3.08959955e-01 4.39114581e+00 4.02954496e+02 9.02677035e-01 2.27647259e+00 4.53994668e+02 3.74379310e-01 4.02432791e+00 6.15694190e+02 4.51943494e-01 4.06522919e+00 6.63307635e+02 1.05793098e+00 2.32168786e+00 7.10644233e+02 4.33194434e-01 4.17050014e+00 8.61276198e+02 1.79240633e-01 4.26617114e+00 9.06211127e+02 5.03070470e-01 2.10563379e+00 9.50973864e+02 1.78487912e-01 4.01254815e+00]和