计算两个叠加的高斯函数的质心

Posted 2023-03-12

【中文标题】计算两个叠加的高斯函数的质心

【英文标题】：Calculating the centroids of two superposed gaussian functions

【发布时间】：2019-03-20 11:53:47

【问题描述】：

我正在尝试找到以下问题的解决方案。我有一组点，它们应该对以不同点为中心的 2 个高斯函数的总和进行建模。我需要找到这两点。到目前为止，我的方法一直是找到整个集合的质心，并在其下方和上方切割日期集。然后我计算每件的质心，这些就是我的中心。然而，这种方法会减少泄漏到数据右半部分的左高斯信息。这使得当高斯靠近在一起时过程会失败。有没有办法更智能地做到这一点？由于计算困难，我希望解决方案不涉及曲线拟合。

【问题讨论】：

我有一个 Python 示例，该示例使用遗传算法将两个洛伦兹峰拟合到碳纳米管的拉曼光谱，以在此处自动提供初始参数估计，用您自己的数据和方程替换它可能会很容易工作：bitbucket.org/zunzuncode/RamanSpectroscopyFit

对于 FFT 中的频率查找，有一种使用高斯窗函数的方法。在对数刻度中，您可以从两个连续的 bin 计算中心。然而，这仅适用于单个峰值，如果接近最大值则更好。你可以从高峰尝试更多。但是如果两个高斯分布强烈混合，很难看出如何避免拟合。不过，我没有看到计算困难。

这是混合高斯模型。找到这种模型的参数的常规方法是所谓的期望最大化算法，它只是您已经拥有的切割和拟合过程的概括。网络搜索“混合高斯期望最大化”应该会找到很多资源。许多编程语言中都有实现此算法的包。

【参考方案1】：

由于 OP 没有显示任何数据，因此不清楚数据的噪声程度。此外，尚不清楚此处如何定义“靠近在一起”。在下文中，我有一个简单的近似值，它适用于低噪声，并假设左侧数据由左高斯主导，而右侧数据由右高斯主导。这对位置、高度，尤其是标准偏差有一些限制。

对于单峰确实有效，但对于混合双峰来说还可以（在上述限制范围内）

#!/usr/bin/python
import matplotlib.pyplot as plt
import numpy as np

def gaussian( x, x0, s, a):
    return a * np.exp( -( x - x0 )**2 / ( 2 * s**2 ) ) / np.sqrt( 2 * np.pi * s**2 )

def get_x0( x1, x2, x3, y1, y2, y3 ):
    l12= np.log( y1 / y2 )
    l13= np.log( y1 / y3 )
    return ( ( x2 + x1 )/2. - ( x3 + x1 )/2. * l12/l13 * ( x3 - x1 ) / ( x2 - x1 ) ) / ( 1 - l12 / l13 * (x3 - x1 ) / ( x2 - x1 ) )


fig = plt.figure( )
ax = fig.add_subplot( 2, 1, 1 )

xL = np.linspace(-8, 8, 150 )
yL = np.fromiter( ( gaussian( x,-2.1, 1.2, 8 ) for x in xL ), np.float )
marker=[10,15,20]

x1 = xL[ marker[0] ]
x2 = xL[ marker[1] ]
x3 = xL[ marker[2] ]
y1 = yL[ marker[0] ]
y2 = yL[ marker[1] ]
y3 = yL[ marker[2] ]

print get_x0( x1, x2, x3, y1, y2, y3 )

ax.plot( xL, yL )
ax.scatter( [ x1, x2, x3 ],[ y1, y2, y3 ])

bx = fig.add_subplot( 2, 1, 2 )
yL = np.fromiter( ( gaussian( x,-2.1, 1.2, 8) + gaussian( x,0.7, 1.4, 6 ) for x in xL ), np.float )
marker=[10,15,20]
x1 = xL[ marker[0] ]
x2 = xL[ marker[1] ]
x3 = xL[ marker[2] ]
y1 = yL[ marker[0] ]
y2 = yL[ marker[1] ]
y3 = yL[ marker[2] ]
bx.scatter( [ x1, x2, x3 ],[ y1, y2, y3 ])
print get_x0( x1, x2, x3, y1, y2, y3 )
marker=[-20,-25,-30]
x1 = xL[ marker[0] ]
x2 = xL[ marker[1] ]
x3 = xL[ marker[2] ]
y1 = yL[ marker[0] ]
y2 = yL[ marker[1] ]
y3 = yL[ marker[2] ]
bx.scatter( [ x1, x2, x3 ],[ y1, y2, y3 ])
print get_x0( x1, x2, x3, y1, y2, y3 )
bx.plot( xL, yL )

plt.show() 

演出：

#Single
-2.0999999999999455
#Double
-2.0951188129317813
0.6998760921436634

与-2.1 和0.7 非常接近

如果出现噪音，可能需要进行一些平均。