我可以拟合质心和峰高随机分布的多个高斯谱吗？

Question

这里是新手，但我已尝试在发布之前进行尽职调查。对任何无意的失礼深表歉意。

我正在以电压与时间序列的形式从示波器采集数据。时间箱的宽度为 0.8 纳秒。我运行多个“数据捕获”周期。单个捕获将具有固定数量的样本，以及 5 到 15 个高斯峰，确切的峰数未知。高斯峰具有相对受限的 FWHM（在 2 到 3 纳秒之间）、变化的峰高和随机到达时间（即质心位置不是周期性的）。

我一直在使用 Python 对这些数据进行高斯拟合，并且使用 scipy.optimise 库和 astropy 库取得了一些成功。下面包含使用 scipy.optimise 的代码。我可以拟合多个高斯，但我的代码中的一个关键步骤是提供峰数的“猜测”，并为每个峰估计质心位置、峰高和峰宽。有没有办法概括这段代码而不必提供“猜测”？如果我放松“猜测”中的条件，那么拟合就会失去质量。我知道这些峰将是具有良好约束宽度的高斯峰，但我想推广代码以适应任何给定数据捕获中的峰质心和峰高。

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

#Get data from file
with open('test3.txt') as f:
    w, h = [float(x) for x in next(f).split()]
    print(w, h)
    array = [[float(x) for x in line.split()] for line in f]

#Separate    
x, z = zip(*array)
#Change sign since fitting routine seems to
#prefer positive numbers
y=[ -p for p in z]

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 the peak positions, heights, and widths
guess = [16, 5, 2, 75, 5, 2, 105, 5, 2, 139, 5, 2, 225, 5, 2, 315, 5, 2, 330, 5, 2]

#Fit and print parameters to screen
popt, pcov = curve_fit(func, x, y, p0=guess)
print(popt)
fit = func(x, *popt)

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

结果如下所示： Plot of Data and Fits

数据文件在这里：https://spaces.hightail.com/receive/5MY7Vc7r9R

这类似于How can I fit multiple Gaussian curved to mass spectrometry data in Python? 和fit multiple gaussians to the data in python，但是这两个依赖于拟合周期性数据集或具有已知峰位置、宽度和高度的数据集。他们让我走到这一步很有用，但我现在被困住了。有什么想法或建议我可以跟进？

谢谢， 阿迪

Answer 1

我的想法是我们将曲线的值与其平均值进行比较。
multiplier 变量表示该值必须大于平均值多少倍才能让我们了解这是峰值之一。超过该值的峰值的第一个点被认为是逼近该峰值平均值的起点。
我还将列表替换为 x 和 y 的数组。

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

#Get data from file
with open('test3.txt') as f:
    w, h = [float(x) for x in next(f).split()]
    array = [[float(x) for x in line.split()] for line in f]

#Separate    
x, z = zip(*array)
x = np.array(x)
y = np.array([ -p for p in z])
#Change sign since fitting routine seems to
#prefer positive numbers


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 the peak positions, heights, and widths
# guess = [16, 5, 2, 75, 5, 2, 105, 5, 2, 139, 5, 2, 225, 5, 2, 315, 5, 2, 330, 5, 2]

def getPeaks(x, y, multiplier):
    x_peaks = []
    isPeak = False
    for i, j in zip(x, y):
        if j > y.mean() * multiplier:
            if not isPeak:
                isPeak = True
                x_peaks.append(i)
        else:
            isPeak = False
    return x_peaks

multiplier = 3
x_peaks = getPeaks(x, y, multiplier)

guess = []
for i in range(len(x_peaks)):
    guess.append(x_peaks[i])
    guess.append(5)
    guess.append(2)
    

#Fit and print parameters to screen
popt, pcov = curve_fit(func, x, y, p0=guess)
print(popt)
fit = func(x, *popt)

#Plot stuff
plt.plot(x, y)
plt.plot(x, fit , 'r--')
# plt.plot(popt[::3], np.ones_like(popt[::3]) * multiplier, 'ko')
plt.show()

Answer 2

正如 cmets 中提到的，每个迭代算法估计都需要从一些超参数开始。在您描述的问题中，您有初始高斯参数和高斯数。 在估计高斯分布时，EM 算法被证明是收敛的。我建议将它与随机初始高斯参数和网格搜索一起使用，以寻找分布数量的最佳解决方案。从 5 到 15 开始，计算每个解的距离并取最小距离解。 (https://en.m.wikipedia.org/wiki/Expectation%E2%80%93maximization_algorithm)