在表达较少的双峰数据上拟合两个高斯

Posted 2023-03-12

【中文标题】在表达较少的双峰数据上拟合两个高斯

【英文标题】：Fitting two Gaussians on less expressed bimodal data

【发布时间】：2018-04-17 00:58:23

【问题描述】：

我试图在双峰分布data 上拟合两个高斯分布，但大多数优化器总是根据如下开始猜测给我错误的结果

我还从scikit-learn 尝试了GMM，但没有多大帮助。我想知道我可能做错了什么以及更好的方法，以便我们可以测试和拟合双峰数据。使用curve_fit和data的示例代码之一如下

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

def gauss(x,mu,sigma,A):
    return A*np.exp(-(x-mu)**2/2/sigma**2)

def bimodal(x,mu1,sigma1,A1,mu2,sigma2,A2):
    return gauss(x,mu1,sigma1,A1)+gauss(x,mu2,sigma2,A2)

def rmse(p0):
    mu1,sigma1,A1,mu2,sigma2,A2 =p0
    y_sim = bimodal(x,mu1,sigma1,A1,mu2,sigma2,A2)
    rms = np.sqrt((y-y_sim)**2/len(y))

data = pd.read_csv('data.csv')
x, y = data.index, data['24hr'].values

expected=(400,720,500,700,774,150)

params,cov=curve_fit(bimodal,x,y,expected, maxfev=100000)
sigma=np.sqrt(np.diag(cov))
plt.plot(x,bimodal(x,*params),color='red',lw=3,label='model')
plt.plot(x,y,label='data')
plt.legend()
print(params,'\n',sigma)

【问题讨论】：

你为什么使用高斯函数？大山峰高度倾斜，看起来几乎是三角形。

我使用了高斯分布，因为后面阶段的数据会扩散到松散的双峰，看起来像高斯分布。我想保持拟合参数一致，以便在后期进行数据间比较。

我没有尝试运行你的代码，但是从结果来看，你可能没有做错任何事情。假设其中一种模式与右侧的小峰很好地对齐。然后另一个高斯模式必须以某种方式适应大的、倾斜的、非高斯峰值，如果没有相当高的 RMS 误差，它就无法做到这一点。尽管不能很好地拟合小峰值，但将两种模式组合在一起以降低匹配大峰值所需的误差时，总体误差可能会小得多。

感谢您指出这个方向。那么可能我需要更改双峰分布函数中的分布类型。我会调查的。

我有一个使用 scipy 的 scipy.optimize.differential_evolution 遗传算法来确定初始参数的示例，用于将双洛伦兹峰方程拟合到碳纳米管的拉曼光谱数据bitbucket.org/zunzuncode/RamanSpectroscopyFit - 替换方程和数据用你自己的，你应该完成。

【参考方案1】：

你可以试试skewed Gaussian。使用参数alpha->0，这变成了一个正常的高斯，允许进行比较：

import matplotlib.pyplot as plt
import numpy as np
from scipy.special import erf
from scipy.optimize import minimize,leastsq, curve_fit


def gauss(x):
    return np.exp( -0.5 * x**2 / np.sqrt( 2 * np.pi ) )


def Phi(x):
    return ( 0.5 * ( 1. + erf(x/np.sqrt(2) ) ) )


def skewed(x, x0, s, a):
    return 2./s * gauss( ( x - x0 ) / s ) * Phi( a * ( x - x0 ) / s)


def my_double_peak(x, A0, x0, s0, a, A1, x1, s1):
    return A0 * skewed( x, x0, s0, a ) + A1 / s1 * gauss( ( x - x1 ) / s1 )

data = np.loadtxt("data.csv", skiprows=1, delimiter=',')
xData = range(len(data))

fitResult, ier = curve_fit( my_double_peak, xData, data[:,1], p0=(45e3, 400., 60,4. ,15e3, 700., 30 )  ) 
print fitResult
bestFit = [my_double_peak(x, *fitResult ) for x in range(len(data)) ]


fig1=plt.figure(1)
ax= fig1.add_subplot( 1, 1, 1 )
ax.plot( data[:,1] )
ax.plot( bestFit )

plt.show()

提供：

>>> [  6.77971459e+04   3.48661227e+02   8.60938473e+01   
       8.43422033e+00   3.86660495e+03   7.22528635e+02   
       2.49055201e+01]

【讨论】：

您使用什么方法来确定代码中 p0 初始参数的值？

@JamesPhillips 目前我只是使用有根据的猜测。对于自动化方法，我需要了解更多有关数据的信息，例如，可以假设未倾斜的峰值始终在右侧，或者倾斜始终具有相同的符号等。这不是万无一失的，但它相当健壮。例如。 p0=(1700, 400., 1., 1. , 1700 / 3., 2 * 400., 1. ) 有效。因此，您只需要知道近似最大值并假设第二个高斯更小并且在右侧近似 2 * x0，设置 sigmas 和偏度 1。其他简单的假设可能适用于真实数据。

谢谢 - 根据我所见，您有根据的猜测非常有效。