第七色在线视频,2021少妇久久久久久久久久,亚洲欧洲精品成人久久av18,亚洲国产精品特色大片观看完整版,孙宇晨将参加特朗普的晚宴

為了賬號(hào)安全,請(qǐng)及時(shí)綁定郵箱和手機(jī)立即綁定
已解決430363個(gè)問題,去搜搜看,總會(huì)有你想問的

在 Python 3 中使用 Scipy 將多個(gè)洛倫茲擬合到布里淵光譜

在 Python 3 中使用 Scipy 將多個(gè)洛倫茲擬合到布里淵光譜

DIEA 2022-11-09 16:24:50
我正在嘗試使用 scipy.optimize.curve_fit 擬合布里淵光譜(有幾個(gè)峰)。我有多個(gè)具有多個(gè)峰的光譜,我正在嘗試將它們與洛倫茲函數(shù)(每個(gè)峰一個(gè)洛倫茲)擬合。我正在嘗試自動(dòng)化批量分析的過程(即,使用 scipy 的峰值查找算法來獲取峰值位置、峰值寬度和峰值高度,并將它們用作擬合的初始猜測)。我現(xiàn)在正在研究一個(gè)光譜,看看總體思路是否有效,然后我會(huì)將其擴(kuò)展為自動(dòng)并使用我擁有的所有光譜。到目前為止,我已經(jīng)這樣做了:import numpy as npimport matplotlib.pyplot as pltfrom scipy.signal import find_peaksfrom 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 peaksI = properties ['peak_heights'] #amplitudefwhm = (properties['widths']) #full width half maximum #define sum of lorentziansdef 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 arrayguess=np.array(guess)#fitpopt, pcov = curve_fit(func, x_data, y_data, p0=guess, method = 'lm',  maxfev=1000000)print(popt)fit = func(x_data, *popt)#plottingfig= 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 )和像素是自變量(任意在代碼中創(chuàng)建的值數(shù)組)。這就是我得到的:頻譜擬合的結(jié)果。知道為什么最后一個(gè)三胞胎沒有按預(yù)期安裝嗎?我檢查并通過 scipy.signal.find_peaks 函數(shù)正確找到了所有峰值 - 因此我認(rèn)為問題出在 scipy.optimize.curve_fit 上,因?yàn)槲冶仨氃黾?maxfev 的數(shù)量才能“工作”。關(guān)于如何以更聰明的方式解決這個(gè)問題的任何想法?
查看完整描述

1 回答

?
慕田峪4524236

TA貢獻(xiàn)1875條經(jīng)驗(yàn) 獲得超5個(gè)贊

通過一些小的修改,OP 的代碼運(yùn)行得很好?,F(xiàn)在看起來像:


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]

http://img1.sycdn.imooc.com//636b63ef00011c0709240463.jpg

查看完整回答
反對(duì) 回復(fù) 2022-11-09
  • 1 回答
  • 0 關(guān)注
  • 258 瀏覽
慕課專欄
更多

添加回答

舉報(bào)

0/150
提交
取消
微信客服

購課補(bǔ)貼
聯(lián)系客服咨詢優(yōu)惠詳情

幫助反饋 APP下載

慕課網(wǎng)APP
您的移動(dòng)學(xué)習(xí)伙伴

公眾號(hào)

掃描二維碼
關(guān)注慕課網(wǎng)微信公眾號(hào)