高斯拟合 python 中的直方图数据：Trust Region v/s Levenberg Marquardt答案

【问题标题】：Gaussian fit to a histogram data in python: Trust Region v/s Levenberg Marquardt高斯拟合 python 中的直方图数据：Trust Region v/s Levenberg Marquardt
【发布时间】：2014-07-08 15:37:29
【问题描述】：

我的直方图清楚地显示了两个峰值。但是，当用双高斯曲线拟合它时，它只显示一个峰值。遵循stackoverflow中显示的几乎所有答案。但未能得到正确的结果。我的老师以前在 Fortran 中做过，他得到了两个峰值。我在一次试验中使用了 python 的scipy.optimize 的leastsq。我也应该提供我的数据吗？这是我的代码。

binss = (max(x) - min(x))/0.05 #0.05 is my bin width
n, bins, patches = plt.hist(x, binss, color = 'grey') #gives the histogram

x_a = []
for item in range(len(bins)-1):
    b = (bins[item]+bins[item+1])/2
    x_a.append(b)

x_avg = np.array(x_a)
y_real = n

def gauss(x, A, mu, sigma):
    gaus = []
    for item in range(len(x)):
        gaus.append(A*e**(-(x[item]-mu)**2./(2.*sigma**2)))
    return np.array(gaus)
A1, A2, m1, m2, sd1, sd2 = [25, 30, 0.3, 0.6, -0.9, -0.9]

#Initial guesses for leastsq
p = [A1, A2, m1, m2, sd1, sd2]
y_init = gauss(x_avg, A1, m1, sd1) + gauss(x_avg, A2, m2, sd2)    #initially guessed y

def residual(p, x, y):
    A1, A2, m1, m2, sd1, sd2 = p
    y_fit = gauss(x, A1, m1, sd1) + gauss(x, A2, m2, sd2)
    err = y - y_fit
    return err

sf = leastsq(residual, p, args = (x_avg , y_real))

y_fitted1 = gauss(x_avg, sf[0][0], sf[0][2], sf[0][4])
y_fitted2 = gauss(x_avg, sf[0][1], sf[0][3], sf[0][5])

y_fitted = y_fitted1 + y_fitted2

plt.plot(x_avg, y_init, 'b', label='Starting Guess')
plt.plot(x_avg, y_fitted, color = 'red', label = 'Fitted Data')
plt.plot(x_avg, y_fitted1, color= 'black', label = 'Fitted1 Data')
plt.plot(x_avg, y_fitted2, color = 'green', label = 'Fitted2 Data')

即使我得到的数字也不平滑。 x_avg 中只有 54 分，请帮忙。甚至不能在这里发布数字。

在 MATLAB 上绘图时，得到了正确的结果。原因： MATLAB 使用 Trust Region 算法而不是 Levenberg-Marquardt 算法这不适合绑定约束。

只有当它显示为 3 的总和时才会出现正确的结果单个高斯，而不是 2。

我如何决定使用哪种算法以及何时使用？

【问题讨论】：

标签： python python-2.7 curve-fitting gaussian

【解决方案1】：

您的问题似乎与mixtures of Gaussian 也称为Gaussian mixture model 有关。有几种实现方式。 sklearn 值得考虑。

import numpy as np
from sklearn import mixture
import matplotlib.pyplot as plt

comp0 = np.random.randn(1000) - 5 # samples of the 1st component
comp1 = np.random.randn(1000) + 5 # samples of the 2nd component

x = np.hstack((comp0, comp1)) # merge them

gmm = mixture.GMM(n_components=2) # gmm for two components
gmm.fit(x) # train it!

linspace = np.linspace(-10, 10, 1000)

fig, ax1 = plt.subplots()
ax2 = ax1.twinx()

ax1.hist(x, 100) # draw samples
ax2.plot(linspace, np.exp(gmm.score_samples(linspace)[0]), 'r') # draw GMM
plt.show()

输出是

【讨论】：

我可以给你沿 x_axis 的数据吗？我仍然得到one_peak_gaussian。感谢您分享您的知识。
当然，分享你的数据！
下载整个捆绑包。 wakari.io/sharing/bundle/MycrofD/data
就像上面给出的答案一样，有两个直方图。如何为最右侧的直方图显示具有 teo 峰值的高斯曲线。它应该有两个峰而不是一个峰，这不是很清楚吗？
在您的数据上，单个高斯似乎优于两个高斯。

【解决方案2】：

我添加了另一个高斯项。所以p 总共取了 9 个参数。因此

p = [A1, A2, A3, m1, m2, m3, sd1, sd2, sd3]

然后另一个术语y_fitted3 被定义并添加到y_fitted。然后它给出了一个完美拟合的两个峰值的正确数字，除了曲线根本不平滑的事实！然后在stackoverflow中搜索导致我使用spline。即

from scipy.interpolate import spline

然后在最后，

x_new = np.linspace(x_avg.min(),x_avg.max(),30000)
ysmooth = spline(x_avg, y_fitted, x_new)
plt.plot(x_new, ysmooth)

然后就到了。检查维基百科，它说python中的L-M也使用T-R。所以再次尝试leastsq 给出了结果。但是，我仍然不清楚MATLAB 中显示的差异。额外的投入将不胜感激！谢谢。

【讨论】：