问题描述

我正在尝试将高斯与噪声吸收光谱相匹配。然而，它似乎并不适用于所有情况。当我尝试将峰值宽度减小到例如PEAK_WIDTH=10时，下面的代码没有产生很好的匹配，只有一行。同样，如果我将峰值的位置再向右移动x_Peak_loc=160，它也不起作用。我如何才能更好地适应这些情况呢？谢谢!代码如下：

import numpy as np
from scipy.optimize import curve_fit
import matplotlib.pyplot as mpl
import matplotlib.pyplot as plt
import scipy.integrate as integrate
def func(x, a, x0, sigma):
    #return a*(1/(np.sqrt(2*np.pi*sigma**2 ))) *np.exp(-(x-x0)**2/(2*sigma**2))
    return a*np.exp(-(x-x0)**2/(2*sigma**2))
amplitude=-10
peak_width=30
x_peak_loc=70

# Generating clean data
x = np.linspace(0, 200, 1000)
y = func(x, amplitude,x_peak_loc, peak_width)
# Adding noise to the data
mn = 0
N=0.2
std=np.sqrt(N)
noise2=np.random.normal(mn,std,size=len(x))
yn = y + noise2
fig = mpl.figure(1)
ax = fig.add_subplot(111)
ax.plot(x, y, c='k', label='analytic function')
ax.scatter(x, yn, s=5, label='fake noisy data')
fig.savefig('model_and_noise.png')
popt, pcov = curve_fit(func, x, yn)
print (popt)

ym = func(x, popt[0], popt[1], popt[2])
ax.plot(x, ym, c='r', label='Best fit')
ax.legend()
fig.savefig('model_fit.png')
plt.legend(loc='upper left')
plt.xlabel("v")
plt.ylabel("f(v)")

推荐答案

你说光谱，所以我猜有不止一个峰。在这种情况下，scipy.signal.find_peaks()可能非常有用。在可以分离峰值的情况下，绝对不需要机器学习过度，因为它可以通过Jean Jacquelin

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

from scipy.integrate import cumtrapz

def func( x, a, x0, sigma ):
    return a * np.exp( -( x - x0 )**2 / ( 2 * sigma**2 ) )
    
amplitude = -10
peak_width = 30
x_peak_loc = 160

# Generating clean data
xl = np.linspace( 0, 200, 1000 )
y0 = func( xl, amplitude, x_peak_loc, peak_width )

mn = 0
N = 0.2
std = np.sqrt( N )
noise2 = np.random.normal( mn, std, size=len( xl ) )
yl = y0 + noise2

"""
Most simple implementation of the linear fit of mu and sigma
"""
nul = yl - yl[0]

Sk = cumtrapz( yl, x=xl, initial=0 )
Tk = cumtrapz( xl * yl, x=xl, initial=0 )

MX = [
        [ np.dot( Sk, Sk ), np.dot( Sk, Tk ) ],
        [ np.dot( Sk, Tk ), np.dot( Tk, Tk ) ]
    ]
Vek = [ np.dot( nul, Sk ), np.dot( nul, Tk ) ]

res = np.dot( np.linalg.inv( MX ), Vek )
sig = np.sqrt( -1 / res[1] )
mu = res[0] * sig**2 
print( sig )
print( mu )
"""
Most simple linear fit of amplitude, probably not required but...hey.
"""
fk = func( xl, 1, mu, sig )
amp = np.dot( yl, fk ) / np.dot( fk, fk )
print( amp )
print( "probably quite close, already")

popt, pcov = curve_fit( func, xl, yl, p0=( sig, mu, amp ) ) 
print( popt )

fig = plt.figure( 1 )
ax = fig.add_subplot( 1, 1, 1 )
ax.plot( xl, y0, c='k', label="analytic function" )
ax.scatter( xl, yl, s=5, label="fake noisy data" )
ym = func( xl, *popt )
ax.plot( xl, ym, c='r', label="Best fit" )
ax.legend()
fig.savefig( "model_fit.png" )
plt.legend( loc="lower left" )
plt.xlabel( "v" )
plt.ylabel( "f(v)" )
plt.show()

这篇关于使用Python对较低的峰宽进行不精确的高斯拟合的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持编程学习网！