使用 Python 的低峰宽的高斯拟合不准确答案

【问题标题】：gaussian fitting inaccurate for lower peak width using Python使用 Python 的低峰宽的高斯拟合不准确
【发布时间】：2021-11-29 14:27:24
【问题描述】：

我正在尝试将高斯拟合到嘈杂的吸收光谱。但是，它似乎不适用于所有情况。当我尝试将峰宽减小到例如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)")

【问题讨论】：

标签： python optimization curve-fitting gaussian

【解决方案1】：

嗯，你说的是光谱，所以我猜想有不止一个峰。在这种情况下，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()

【讨论】：

有趣。感谢分享 Jacquelin 的方法。他的方法在我的系统上也快 10-20 倍：%timeit jacquelin_method(xl, yl)245 µs ± 8.15 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each) 与 %timeit curvefit_method(xl, yl)5 ms ± 82.4 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)。顺便说一句，最新版本的 SciPy 中不再存在 cumtrapz。您可能可以在最新版本中使用scipy.integrate.quad。
@Roald 我在 Python3 上有 1.5.4 版，cumtrapz() 存在。在任何情况下，quad() 需要一个可调用的而不是数据列表，所以这是行不通的。如果你真的没有，累积列表集成应该不难写。你可以做一个样条拟合，或者其他一些分段的东西来获得一个可调用的并将其提供给四边形。不过，这不是我最喜欢的。
对于我拥有的 1.6.1，cumtrapz() 也存在。我只是在寻找文档，只能找到scipy v=0.18.1 的页面。没有最新版本（1.7.1）的文档，所以我认为它已被删除。但也许它仍然存在。
@Roald 我已经看到了，我想尝试点更新，这将我从 1.5.3 更改为 1.5.4。 cumtrapz 出席。虽然不知道宣布的 1.7.1。
@Roald 我猜它还在。在“使用样本集成”部分的 1.7.1 文档中，“...梯形和辛普森...”已链接。第一个然后显示为下一页 scipy.integrate.cumulative_trapezoid。它还在那里；在 ODE 的部分...不知道为什么，搜索功能不起作用。

【解决方案2】：

我认为问题在于使用curve_fit 或任何基于梯度的算法。在处理嘈杂的数据时，他们往往会陷入局部最小值。有时你很幸运，有时你不是。如果增加噪声幅度，您将获得与改变峰宽相同的（不幸的）效果。

您应该提供一些好的初始猜测或将拟合方法更改为完全不同的东西。您需要使用一些hyperparameter optimization 方法搜索函数参数的空间。

也许一些简单的级联（精度越来越高）网格搜索就足够了？

【讨论】：