python - SciPy 中使用截断法线的混合模型拟合(双峰？)。 python 3

Question

我试图以此为例，但似乎无法使其适应我的数据集，因为我需要截断的法线: https://stackoverflow.com/questions/35990467/fit-two-gaussians-to-a-histogram-from-one-set-of-data-python#=

我有一个数据集，它绝对是 2 个截断法线的混合。域中的最小值为 0，最大值为 1。我想创建一个可以拟合的对象来优化参数并获取从该分布中抽取的数字序列的可能性。一种选择可能是仅使用 KDE 模型并使用 pdf 来获取可能性。但是，我想要两个分布的确切平均值和标准差。我想我可以将数据分成两半，然后分别对 2 个法线进行建模，但我也想学习如何在 SciPy 中使用 optimize。我刚刚开始尝试这种类型的统计分析，所以如果这看起来很天真，我深表歉意。

我不知道如何以这种方式获得一个可以积分为 1 并且域限制在 0 和 1 之间的 pdf。

import requests
from ast import literal_eval
from scipy import optimize, stats
import matplotlib.pyplot as plt
import seaborn as sns
import numpy as np


# Actual Data
u = np.asarray(literal_eval(requests.get("https://pastebin.com/raw/hP5VJ9vr").text))
# u.size ==> 6000
u.min(), u.max()
# (1.3628525454666037e-08, 0.99973136607553781)

# Distribution
with plt.style.context("seaborn-white"):
    fig, ax = plt.subplots()
    sns.kdeplot(u, color="black", ax=ax)
    ax.axvline(0, linestyle=":", color="red")
    ax.axvline(1, linestyle=":", color="red")
kde = stats.gaussian_kde(u)

# KDE Model
def truncated_gaussian_lower(x,mu,sigma,A):
    return np.clip(A*np.exp(-(x-mu)**2/2/sigma**2), a_min=0, a_max=None)
def truncated_gaussian_upper(x,mu,sigma,A):
    return np.clip(A*np.exp(-(x-mu)**2/2/sigma**2), a_min=None, a_max=1)
def mixture_model(x,mu1,sigma1,A1,mu2,sigma2,A2):
    return truncated_gaussian_lower(x,mu1,sigma1,A1) + truncated_gaussian_upper(x,mu2,sigma2,A2)
kde = stats.gaussian_kde(u)

# Estimates: mu sigma A
estimates= [0.1, 1, 3, 
            0.9, 1, 1]
params,cov= optimize.curve_fit(mixture_model,u,kde.pdf(u),estimates )

# ---------------------------------------------------------------------------
# RuntimeError                              Traceback (most recent call last)
# <ipython-input-265-b2efb2ca0e0a> in <module>()
#      32 estimates= [0.1, 1, 3, 
#      33             0.9, 1, 1]
# ---> 34 params,cov= optimize.curve_fit(mixture_model,u,kde.pdf(u),estimates )

# /Users/mu/anaconda/lib/python3.6/site-packages/scipy/optimize/minpack.py in curve_fit(f, xdata, ydata, p0, sigma, absolute_sigma, check_finite, bounds, method, jac, **kwargs)
#     738         cost = np.sum(infodict['fvec'] ** 2)
#     739         if ier not in [1, 2, 3, 4]:
# --> 740             raise RuntimeError("Optimal parameters not found: " + errmsg)
#     741     else:
#     742         # Rename maxfev (leastsq) to max_nfev (least_squares), if specified.

# RuntimeError: Optimal parameters not found: Number of calls to function has reached maxfev = 1400.

回应@Uvar 下面非常有用的解释。我正在尝试测试 0 - 1 的积分，看看它是否等于 1，但我得到的是 0.3。我认为我错过了逻辑上的关键一步:

# KDE Model
def truncated_gaussian(x,mu,sigma,A):
    return A*np.exp(-(x-mu)**2/2/sigma**2)

def mixture_model(x,mu1,sigma1,A1,mu2,sigma2,A2):
    if type(x) == np.ndarray:
        norm_probas = truncated_gaussian(x,mu1,sigma1,A1) + truncated_gaussian(x,mu2,sigma2,A2)
        mask_lower = x < 0
        mask_upper = x > 1
        mask_floor = (mask_lower.astype(int) + mask_upper.astype(int)) > 1
        norm_probas[mask_floor] = 0
        return norm_probas
    else:
        if (x < 0) or (x > 1):
            return 0
        return truncated_gaussian_lower(x,mu1,sigma1,A1) + truncated_gaussian_upper(x,mu2,sigma2,A2)

kde = stats.gaussian_kde(u, bw_method=2e-2)

# # Estimates: mu sigma A
estimates= [0.1, 1, 3, 
            0.9, 1, 1]
params,cov= optimize.curve_fit(mixture_model,u,kde.pdf(u)/integrate.quad(kde, 0 , 1)[0],estimates ,maxfev=5000)
# params
# array([  9.89751700e-01,   1.92831695e-02,   7.84324114e+00,
#          3.73623345e-03,   1.07754038e-02,   3.79238972e+01])

# Test the integral from 0 - 1
x = np.linspace(0,1,1000)
with plt.style.context("seaborn-white"):
    fig, ax = plt.subplots()
    ax.plot(x, kde(x), color="black", label="Data")
    ax.plot(x, mixture_model(x, *params), color="red", label="Model")
    ax.legend()
# Integrating from 0 to 1
integrate.quad(lambda x: mixture_model(x, *params), 0,1)[0]
# 0.3026863969781809

Answer 1

您似乎错误地指定了拟合程序。您正在尝试在限制半边界的同时适应kde.pdf(u)。

foo = kde.pdf(u)

min(foo)
Out[329]: 0.22903365654960098

max(foo)
Out[330]: 4.0119283429320332

如您所见，u 的概率密度函数不限于 [0,1]。因此，只需删除裁剪操作即可实现精确拟合。

def truncated_gaussian_lower(x,mu,sigma,A):
    return A*np.exp((-(x-mu)**2)/(2*sigma**2))
def truncated_gaussian_upper(x,mu,sigma,A):
    return A * np.exp((-(x-mu)**2)/(2*sigma**2))
def mixture_model(x,mu1,sigma1,A1,mu2,sigma2,A2):
    return truncated_gaussian_lower(x,mu1,sigma1,A1) + truncated_gaussian_upper(x,mu2,sigma2,A2)

estimates= [0.15, 1, 3, 
            0.95, 1, 1]
params,cov= optimize.curve_fit(f=mixture_model, xdata=u, ydata=kde.pdf(u), p0=estimates)

params
Out[327]: 
array([ 0.00672248,  0.07462657,  4.01188383,  0.98006841,  0.07654998,
        1.30569665])

y3 = mixture_model(u, params[0], params[1], params[2], params[3], params[4], params[5])

plt.plot(kde.pdf(u)+0.1) #add offset for visual inspection purpose

plt.plot(y3)

所以，现在假设我改变了我的计划:

plt.figure(); plt.plot(u,y3,'.')

因为，确实:

np.allclose(y3, kde(u), atol=1e-2)
>>True

您可以将混合模型稍微编辑为域外的 0 [0, 1]:

def mixture_model(x,mu1,sigma1,A1,mu2,sigma2,A2):
    if (x < 0) or (x > 1):
        return 0
    return truncated_gaussian_lower(x,mu1,sigma1,A1) + truncated_gaussian_upper(x,mu2,sigma2,A2)

但是，这样做将失去通过 x 数组立即计算函数的选项。因此，为了便于讨论，我现在将其保留。

无论如何，我们希望我们的积分在域 [0, 1] 中总和为 1，并且有一种方法可以做到这一点(请随意使用 中的带宽估计器stats.gaussian_kde 以及..)是将概率密度估计除以其在域上的积分。请注意，在此实现中，optimize.curve_fit 仅需要 1400 次迭代，因此初始参数估计很重要。

from scipy import integrate
sum_prob = integrate.quad(kde, 0 , 1)[0]
y = kde(u)/sum_prob
# Estimates: mu sigma A
estimates= [0.15, 1, 5, 
            0.95, 0.5, 3]
params,cov= optimize.curve_fit(f=mixture_model, xdata=u, ydata=y, p0=estimates)
>>array([  6.72247814e-03,   7.46265651e-02,   7.23699661e+00,
     9.80068414e-01,   7.65499825e-02,   2.35533297e+00])

y3 = mixture_model(np.arange(0,1,0.001), params[0], params[1], params[2], 
    params[3], params[4], params[5])

with plt.style.context("seaborn-white"):
    fig, ax = plt.subplots()
    sns.kdeplot(u, color="black", ax=ax)
    ax.axvline(0, linestyle=":", color="red")
    ax.axvline(1, linestyle=":", color="red")
    plt.plot(np.arange(0,1,0.001), y3) #The red line is now your custom pdf with area-under-curve = 0.998 in the domain..

为了检查曲线下的面积，我使用了重新定义 mix_model 的这个 hacky 解决方案:

def mixture_model(x):
    mu1=params[0]; sigma1=params[1]; A1=params[2]; mu2=params[3]; sigma2=params[4]; A2=params[5]
    return truncated_gaussian_lower(x,mu1,sigma1,A1) + truncated_gaussian_upper(x,mu2,sigma2,A2)

from scipy import integrate
integrated_value, error = integrate.quad(mixture_model, 0, 1) #0 lower bound, 1 upper bound
>>(0.9978588016186962, 5.222293368393178e-14)

或者用第二种方法进行积分:

import sympy
x = sympy.symbols('x', real=True, nonnegative=True)
foo = sympy.integrate(params[2]*sympy.exp((-(x-params[0])**2)/(2*params[1]**2))+params[5]*sympy.exp((-(x-params[3])**2)/(2*params[4]**2)),(x,0,1), manual=True)
foo.doit()

>>0.562981541724715*sqrt(pi) #this evaluates to 0.9978588016186956

实际上按照您编辑的问题中所述的方式进行操作:

def mixture_model(x,mu1,sigma1,A1,mu2,sigma2,A2):
    return truncated_gaussian_lower(x,mu1,sigma1,A1) + truncated_gaussian_upper(x,mu2,sigma2,A2)
integrate.quad(lambda x: mixture_model(x, *params), 0,1)[0]
>>0.9978588016186962

如果我将带宽设置为您的级别 (2e-2)，确实评估会下降到 0.92，这比我们之前的 0.998 更糟糕，但这仍然与您报告的 0.3 有很大不同即使复制您的代码片段，我也无法重新创建一些东西。您是否可能不小心在某处重新定义了函数/变量？