拟合双峰偏态高斯

Question

我有一些数据正在尝试用双峰偏态高斯拟合。我从标准双峰高斯开始，数据太倾斜了。我还希望能够提取底层高斯。

对于双峰，我有（由于循环和绘图，它比需要的多一点）：

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)

params=np.empty(6,dtype=object)

for i in range(len(files)):
    x=data_x[i]
    y=data_y[i]
    
    max_count=np.max(data_y[i][2:])
     
    peaks, _ = find_peaks(data_y[i], height=max_count-5, distance=4)
    
    
    p_init=[data_x[i][peaks[0]],1,max_count,data_x[i][peaks[1]],1,max_count,]
    params,cov=curve_fit(bimodal,x[0:len(y)],y,p_init)
    
    params=np.vstack([params,params])
    x=x.astype(int)
    plt.bar(x[0:len(y)],y)
    
    plt.plot(x,bimodal(x,*params),'r')
    plt.title(naming_string[i]+' -> index = '+ str(i))
    plt.show()

我拥有的数据如下所示：

它确实需要倾斜拟合，左侧曲线具有正倾斜，右侧曲线具有负倾斜。

我想做的只是修改输入函数来做到这一点。然而，所有偏态正态分布都使用 scipy 绕过输入函数。这不允许我将其修改为双峰（我可以告诉）。

我尝试过的函数是scipy函数：

def bimodal_skew_normal(x, amp1, mean1, std1, skew1, amp2, mean2, std2, skew2):
    pdf1 = amp1 * skewnorm.pdf(x, skew1, loc=mean1, scale=std1)
    pdf2 = amp2 * skewnorm.pdf(x, skew2, loc=mean2, scale=std2)
    return pdf1 + pdf2

在对初始参数进行一些调整后，这可能会起作用，但是我不确定如何从中获得底层高斯。脱离拟合的幅度和均值不对应于未偏斜的高斯。

Answer 1

如果没有任何数据集，我们必须首先创建一个合成数据集。 假设您的数据确实遵循双偏态正态分布，那么我们就这样做吧。

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

def peak(x, A, a, loc, scale):
    return A * stats.skewnorm.pdf(x, a=a, loc=loc, scale=scale)

def model(x, A1, a1, loc1, scale1, A2, a2, loc2, scale2):
    return peak(x, A1, a1, loc1, scale1) + peak(x, A2, a2, loc2, scale2)

p0 = (500, 1.2, 140, 10, 600, 1.5, 170, 15)

xexp = np.linspace(120, 220, 70)
yexp = model(xexp, *p0)

现在我们可以用模型拟合这个分布：

popt, pcov = optimize.curve_fit(
    model, xexp, yexp,
    p0=(1000, 1, 150, 1, 1000, 1, 180, 1),
    bounds = [
        (1, 0.1, 140, 0.1, 1, 0.1, 170, 0.1),
        (1e4, 2.0, 160, 100, 1e4, 2.0, 190, 100)
    ]
)
# array([500. ,   1.2, 140. ,  10. , 600. ,   1.5, 170. ,  15. ])

初始猜测和界限很重要。

最终结果，看起来像：

xlin = np.linspace(xexp.min(), xexp.max(), 250)
yhat = model(xlin, *popt)

fig, axe = plt.subplots()
axe.step(xexp, yexp, where="mid")
axe.plot(xlin, yhat)
axe.grid()