Python加载数据并进行多高斯拟合

Question

我一直在寻找一种方法来对我的数据进行多重高斯拟合。到目前为止，我发现的大多数示例都使用正态分布来生成随机数。但我有兴趣查看我的数据图并检查是否有 1-3 个峰值。

我可以对一个峰执行此操作，但我不知道如何对更多峰执行此操作。

例如，我有这个数据:http://www.filedropper.com/data_11

我尝试过使用 lmfit，当然还有 scipy，但没有很好的结果。

感谢您的帮助!

Answer 1

简单地制作单个高斯总和的参数化模型函数。为您的初始猜测选择一个合适的值(这是非常关键的一步)，然后让 scipy.optimize 稍微调整这些数字。

下面是你可能会怎么做:

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

data = np.genfromtxt('data.txt')
def gaussian(x, height, center, width, offset):
    return height*np.exp(-(x - center)**2/(2*width**2)) + offset
def three_gaussians(x, h1, c1, w1, h2, c2, w2, h3, c3, w3, offset):
    return (gaussian(x, h1, c1, w1, offset=0) +
        gaussian(x, h2, c2, w2, offset=0) +
        gaussian(x, h3, c3, w3, offset=0) + offset)

def two_gaussians(x, h1, c1, w1, h2, c2, w2, offset):
    return three_gaussians(x, h1, c1, w1, h2, c2, w2, 0,0,1, offset)

errfunc3 = lambda p, x, y: (three_gaussians(x, *p) - y)**2
errfunc2 = lambda p, x, y: (two_gaussians(x, *p) - y)**2

guess3 = [0.49, 0.55, 0.01, 0.6, 0.61, 0.01, 1, 0.64, 0.01, 0]  # I guess there are 3 peaks, 2 are clear, but between them there seems to be another one, based on the change in slope smoothness there
guess2 = [0.49, 0.55, 0.01, 1, 0.64, 0.01, 0]  # I removed the peak I'm not too sure about
optim3, success = optimize.leastsq(errfunc3, guess3[:], args=(data[:,0], data[:,1]))
optim2, success = optimize.leastsq(errfunc2, guess2[:], args=(data[:,0], data[:,1]))
optim3

plt.plot(data[:,0], data[:,1], lw=5, c='g', label='measurement')
plt.plot(data[:,0], three_gaussians(data[:,0], *optim3),
    lw=3, c='b', label='fit of 3 Gaussians')
plt.plot(data[:,0], two_gaussians(data[:,0], *optim2),
    lw=1, c='r', ls='--', label='fit of 2 Gaussians')
plt.legend(loc='best')
plt.savefig('result.png')

如您所见，这两种拟合(视觉上)几乎没有区别。因此，您无法确定源中是否存在 3 个高斯函数或仅存在 2 个。但是，如果您必须进行猜测，请检查最小残差:

err3 = np.sqrt(errfunc3(optim3, data[:,0], data[:,1])).sum()
err2 = np.sqrt(errfunc2(optim2, data[:,0], data[:,1])).sum()
print('Residual error when fitting 3 Gaussians: {}\n'
    'Residual error when fitting 2 Gaussians: {}'.format(err3, err2))
# Residual error when fitting 3 Gaussians: 3.52000910965
# Residual error when fitting 2 Gaussians: 3.82054499044

在这种情况下，3 个高斯给出了更好的结果，但我也使我的初始猜测相当准确。