Python：使用非线性最小二乘法的双曲线高斯拟合

Question

我的数学知识是有限的，这就是为什么我可能卡住了。我有一个谱图，我试图拟合两个高斯峰。我可以适应最大的高峰，但我无法适应最小的高峰。我知道我需要为两个峰值求和高斯函数，但我不知道我出错的地方。我的电流输出的图像显示：

蓝线是我的数据和绿线是我目前的契合。有个肩膀可以在我的数据主峰左侧我目前努力配合，使用下面的代码：

import matplotlib.pyplot as pt 
import numpy as np 
from scipy.optimize import leastsq 
from pylab import * 

time = [] 
counts = [] 


for i in open('/some/folder/to/file.txt', 'r'): 
    segs = i.split() 
    time.append(float(segs[0])) 
    counts.append(segs[1]) 

time_array = arange(len(time), dtype=float) 
counts_array = arange(len(counts)) 
time_array[0:] = time 
counts_array[0:] = counts 


def model(time_array0, coeffs0): 
    a = coeffs0[0] + coeffs0[1] * np.exp(- ((time_array0-coeffs0[2])/coeffs0[3])**2) 
    b = coeffs0[4] + coeffs0[5] * np.exp(- ((time_array0-coeffs0[6])/coeffs0[7])**2) 
    c = a+b 
    return c 


def residuals(coeffs, counts_array, time_array): 
    return counts_array - model(time_array, coeffs) 

# 0 = baseline, 1 = amplitude, 2 = centre, 3 = width 
peak1 = np.array([0,6337,16.2,4.47,0,2300,13.5,2], dtype=float) 
#peak2 = np.array([0,2300,13.5,2], dtype=float) 

x, flag = leastsq(residuals, peak1, args=(counts_array, time_array)) 
#z, flag = leastsq(residuals, peak2, args=(counts_array, time_array)) 

plt.plot(time_array, counts_array) 
plt.plot(time_array, model(time_array, x), color = 'g') 
#plt.plot(time_array, model(time_array, z), color = 'r') 
plt.show() 

Answer 1

此代码为我工作提供，你只拟合函数是一个两个高斯分布的组合。

我只是做了一个残差函数，它添加了两个高斯函数，然后从真实数据中减去它们。

我传递给Numpy最小二乘函数的参数（p）包括：第一个高斯函数的平均值（m），与第一和第二高斯函数的平均值的差值（dm，即水平位移） ，第一个标准偏差（sd1）和第二个标准偏差（sd2）。

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

###################################### 
# Setting up test data 
def norm(x, mean, sd): 
    norm = [] 
    for i in range(x.size): 
    norm += [1.0/(sd*np.sqrt(2*np.pi))*np.exp(-(x[i] - mean)**2/(2*sd**2))] 
    return np.array(norm) 

mean1, mean2 = 0, -2 
std1, std2 = 0.5, 1 

x = np.linspace(-20, 20, 500) 
y_real = norm(x, mean1, std1) + norm(x, mean2, std2) 

###################################### 
# Solving 
m, dm, sd1, sd2 = [5, 10, 1, 1] 
p = [m, dm, sd1, sd2] # Initial guesses for leastsq 
y_init = norm(x, m, sd1) + norm(x, m + dm, sd2) # For final comparison plot 

def res(p, y, x): 
    m, dm, sd1, sd2 = p 
    m1 = m 
    m2 = m1 + dm 
    y_fit = norm(x, m1, sd1) + norm(x, m2, sd2) 
    err = y - y_fit 
    return err 

plsq = leastsq(res, p, args = (y_real, x)) 

y_est = norm(x, plsq[0][0], plsq[0][2]) + norm(x, plsq[0][0] + plsq[0][1], plsq[0][3]) 

plt.plot(x, y_real, label='Real Data') 
plt.plot(x, y_init, 'r.', label='Starting Guess') 
plt.plot(x, y_est, 'g.', label='Fitted') 
plt.legend() 
plt.show() 

Answer 2

coeffs 0和4退化 - 是绝对没有的，可以决定它们之间的数据。你应该使用一个零水平参数而不是两个（即从你的代码中删除其中的一个）。这可能是阻止你的合适（忽略这里的评论，说这是不可能的 - 这些数据中至少有两个高峰，你当然应该能够适应这一点）。

（可能不太清楚为什么我提出这个建议，但是发生的事情是系数0和4可以相互抵消，它们都可以是零，或者一个可以是100，另一个可以是100这种“适应”就是一样的好，这使得适配程序“混淆”了，它花费了时间试图弄清楚他们应该做什么，什么时候没有单一的正确答案，因为无论价值是什么，其他都可能是负面的，并且合适的将是相同的）。实际上，从情节来看，它可能根本不需要零水平。我会试着放弃这两种，并看看适合的外观。

此外，不需要在最小平方中拟合coeffs 1和5（或零点）。相反，因为模型是线性的，你可以在每个循环中计算它们的值。这会使事情变得更快，但并不重要。我只是注意到你说你的数学不太好，所以可能忽略这个。

Answer 3

可以使用高斯混合模型从scikit-learn：

from sklearn import mixture 
import matplotlib.pyplot 
import matplotlib.mlab 
import numpy as np 
clf = mixture.GMM(n_components=2, covariance_type='full') 
clf.fit(yourdata) 
m1, m2 = clf.means_ 
w1, w2 = clf.weights_ 
c1, c2 = clf.covars_ 
histdist = matplotlib.pyplot.hist(yourdata, 100, normed=True) 
plotgauss1 = lambda x: plot(x,w1*matplotlib.mlab.normpdf(x,m1,np.sqrt(c1))[0], linewidth=3) 
plotgauss2 = lambda x: plot(x,w2*matplotlib.mlab.normpdf(x,m2,np.sqrt(c2))[0], linewidth=3) 
plotgauss1(histdist[1]) 
plotgauss2(histdist[1]) 

您也可以使用下面的功能，以适应您想NCOMP参数高斯数量：

from sklearn import mixture 
%pylab 

def fit_mixture(data, ncomp=2, doplot=False): 
    clf = mixture.GMM(n_components=ncomp, covariance_type='full') 
    clf.fit(data) 
    ml = clf.means_ 
    wl = clf.weights_ 
    cl = clf.covars_ 
    ms = [m[0] for m in ml] 
    cs = [numpy.sqrt(c[0][0]) for c in cl] 
    ws = [w for w in wl] 
    if doplot == True: 
     histo = hist(data, 200, normed=True) 
     for w, m, c in zip(ws, ms, cs): 
      plot(histo[1],w*matplotlib.mlab.normpdf(histo[1],m,np.sqrt(c)), linewidth=3) 
    return ms, cs, ws