Подгонка нескольких лоренцевых к спектру Бриллюэна с помощью Scipy in Python 3

Question

Я пытаюсь подогнать Brillouin Spectra (с несколькими пиками) с помощью scipy.optimize.curve_fit. У меня есть несколько спектров с несколькими пиками, и я пытаюсь подогнать их с помощью лоренцевых функций (по одному лоренцеву на пик). Я пытаюсь автоматизировать процесс массового анализа (то есть, используя алгоритм поиска пиков scipy, чтобы получить положения пиков, ширину пиков и высоту пиков и использовать их в качестве начальных предположений для подгонки). Сейчас я работаю над одним спектром, чтобы увидеть, работает ли общая идея, затем я расширю его до автоматического c и буду работать со всеми имеющимися у меня спектрами. До сих пор я делал это:

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 - зависимая переменная (вот значения в виде файла таблицы Google: 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 np
import matplotlib.pyplot as plt
from scipy.optimize import leastsq
from scipy.signal import find_peaks

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

y0 = np.loadtxt( 'y_peo.tsv' )
yData = np.loadtxt( 'y_peo.tsv' )
xData = np.arange( len( yData ) )

yGround = min( yData )
yData = yData - yGround
yAmp = 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]

и