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

Question

我试图在双峰分布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)

Answer 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]