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()

- Harpal

---

在这种情况下，这将是相当困难的，因为两个峰非常靠近 - 较小的“高斯”没有明确的峰值。通常会（我认为）识别所有感兴趣的峰，然后迭代每个峰，遮罩所有其他峰并适应于每个峰。总拟合结果是所有这些拟合的总和。看起来你需要做的是确定大峰和它的范围，然后在拟合较小峰之前从数据中遮蔽它。 - Chris

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()

- Usagi

---

那么我假设对于n个高斯函数，我需要将它们相加并从数据中减去它们？ - Harpal

@Harpal - 是的。你可以修改代码以使用n个曲线。我只会确保以一种不同的方式编写算法，使得没有两个曲线具有相同的平均值。 - Usagi

该行代码应为 y_est = norm(x, plsq[0][0], plsq[0][2]) + norm(x, plsq[0][0] + plsq[0][1], plsq[0][3])，在您的示例中并不明显，因为其中一个均值为零。已进行编辑。除此之外，解决方案非常好 :) - Kyle

Answer 2

你可以使用来自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

- bougui

---

这将适应数据的直方图，而不是数据本身。 - Rob

Answer 3

coeffs 0和4是退化的-数据中绝对没有任何东西可以区分它们。应该使用单个零级参数代替两个（即从代码中删除一个）。这可能是阻止您匹配的原因（请忽略此处说这不可能的评论-该数据中显然至少有两个峰，您肯定可以适应它）。

（可能不清楚为什么我建议这样做，但正在发生的是coeffs 0和4可以相互抵消。它们都可以为零，也可以为100和另一个为-100-无论哪种方式，匹配效果都一样好。这使“拟合程序”感到困惑，因为它花费时间来尝试确定它们应该是什么，而没有单个正确答案，因为无论其中一个值是什么，另一个都可以是其负数，匹配结果都相同）。

实际上，从图中看来，可能根本不需要零点。我会尝试放弃它们并查看匹配的外观如何。

还有，没有必要在最小二乘法中拟合coeffs 1和5（或零点）。相反，因为该模型在线性方面计算了它们的值，所以每次循环都可以计算它们的值。这将使事情更快，但并非关键。我只是注意到您说自己的数学不太好，所以可能会忽略这个建议。

- andrew cooke

---

尽管有些刻薄，但我认为这实际上是可信的。如果您可以一次性将整个模型适配好，那将有无数的优势。已点赞。 - nes1983