在Python中将多个高斯函数拟合到数据

Question

在Python中将多个高斯函数拟合到数据

11

我在想是否有一种简单的方法来实现高斯/洛伦兹函数拟合10个峰并提取FWHM，同时确定FWHM在x值上的位置。复杂的方法是分离峰并拟合数据并提取FWHM。

数据为[https://drive.google.com/file/d/0B6sUnnbyNGuOT2RZb2UwYXU4dlE/view?usp=sharing]。

非常感谢任何建议。谢谢。

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

data = np.loadtxt('data.txt', delimiter=',')
x, y = data

plt.plot(x,y)
plt.show()

def func(x, *params):
    y = np.zeros_like(x)
    print len(params)
    for i in range(0, len(params), 3):
        ctr = params[i]
        amp = params[i+1]
        wid = params[i+2]
        y = y + amp * np.exp( -((x - ctr)/wid)**2)



guess = [0, 60000, 80, 1000, 60000, 80]
for i in range(12):
    guess += [60+80*i, 46000, 25]


popt, pcov = curve_fit(func, x, y, p0=guess)
print popt
fit = func(x, *popt)

plt.plot(x, y)
plt.plot(x, fit , 'r-')
plt.show()



Traceback (most recent call last):
File "C:\Users\test.py", line 33, in <module>
popt, pcov = curve_fit(func, x, y, p0=guess)
File "C:\Python27\lib\site-packages\scipy\optimize\minpack.py", line 533, in curve_fit
res = leastsq(func, p0, args=args, full_output=1, **kw)
File "C:\Python27\lib\site-packages\scipy\optimize\minpack.py", line 368, in leastsq
shape, dtype = _check_func('leastsq', 'func', func, x0, args, n)
File "C:\Python27\lib\site-packages\scipy\optimize\minpack.py", line 19, in _check_func
res = atleast_1d(thefunc(*((x0[:numinputs],) + args)))
File "C:\Python27\lib\site-packages\scipy\optimize\minpack.py", line 444, in    _ general_function
return function(xdata, *params) - ydata
TypeError: unsupported operand type(s) for -: 'NoneType' and 'float'

- Rocky

1

@LokeshA.R. fwhm 的通常含义是“半峰全宽”。它是衡量光谱峰宽度的一种方便方法。 - John1024

仅供参考，数据已不再可用。 - MortenSickel

3个回答

2

@john1024的回答很好，但需要手动过程来生成初始猜测。这里有一种简单的方法来自动化起始猜测。将john1024代码中相关的3行替换为以下内容：

    import scipy.signal
    i_pk = scipy.signal.find_peaks_cwt(y, widths=range(3,len(x)//Npks))
    DX = (np.max(x)-np.min(x))/float(Npks) # starting guess for component width
    guess = np.ravel([[x[i], y[i], DX] for i in i_pk]) # starting guess for (x, amp, width) for each component

- adriaan peens-hough

-1

在像这样的问题中，我认为始终建议绘制残差（数据-模型）。您还需要查看拟合的ChiSq。

- Doug Wood

1

这条信息应该是一条注释，而不是一个答案 ;) - Gowachin

网页内容由stack overflow 提供, 点击上面的

可以查看英文原文，
原文链接

- John1024 · Accepted Answer

这需要进行非线性拟合。一个很好的工具是scipy的curve_fit函数。

要使用curve_fit，我们需要一个模型函数，称之为func，它将x和我们（猜测的）参数作为参数，并返回相应的y值。我们使用高斯函数的总和作为我们的模型：

from scipy.optimize import curve_fit
import numpy as np

def func(x, *params):
    y = np.zeros_like(x)
    for i in range(0, len(params), 3):
        ctr = params[i]
        amp = params[i+1]
        wid = params[i+2]
        y = y + amp * np.exp( -((x - ctr)/wid)**2)
    return y

现在，让我们为我们的参数创建一个初始猜测。这个猜测从在x=0和x=1,000处具有振峰，振幅为60,000，e-折叠宽度为80开始。然后，我们添加候选峰在x=60, 140, 220, ...处，振幅为46,000，宽度为25：

guess = [0, 60000, 80, 1000, 60000, 80]
for i in range(12):
    guess += [60+80*i, 46000, 25]

现在，我们已经准备好进行拟合：

popt, pcov = curve_fit(func, x, y, p0=guess)
fit = func(x, *popt)

为了看到我们的表现如何，让我们绘制实际的y值(实线黑色曲线)和fit(虚线红色曲线)，并与x对比:

enter image description here

正如您所见，拟合效果相当不错。

完整的工作代码

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

data = np.loadtxt('data.txt', delimiter=',')
x, y = data

plt.plot(x,y)
plt.show()

def func(x, *params):
    y = np.zeros_like(x)
    for i in range(0, len(params), 3):
        ctr = params[i]
        amp = params[i+1]
        wid = params[i+2]
        y = y + amp * np.exp( -((x - ctr)/wid)**2)
    return y

guess = [0, 60000, 80, 1000, 60000, 80]
for i in range(12):
    guess += [60+80*i, 46000, 25]   

popt, pcov = curve_fit(func, x, y, p0=guess)
print popt
fit = func(x, *popt)

plt.plot(x, y)
plt.plot(x, fit , 'r-')
plt.show()