SciPy 中使用截断法线的混合模型拟合（双峰？）。 Python 3
Mixture model fitting (Bimodal?) in SciPy using truncated normals. Python 3
我试图以此为例，但由于我需要截断的法线，因此似乎无法使其适应我的数据集：
我有一个数据集，它绝对是 2 个截断法线的混合体。域中的最小值为 0，最大值为 1。我想创建一个适合优化参数并获得从该分布中提取数字序列的可能性的对象。一种选择可能是只使用 KDE 模型并使用 pdf 来获得可能性。但是，我想要 2 个分布的确切均值和标准差。我想我可以，将数据分成两半，然后分别对 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
您似乎错误地指定了验配程序。
您试图在限制半边界的同时适应 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..
为了检查曲线下的面积，我使用了重新定义 mixture_model..:[=28=]
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 仍然有很大差异这是我无法重新创建的东西，即使在复制您的代码片段时也是如此。您是否可能不小心在某处重新定义了 functions/variables？
我试图以此为例，但由于我需要截断的法线，因此似乎无法使其适应我的数据集：
我有一个数据集，它绝对是 2 个截断法线的混合体。域中的最小值为 0，最大值为 1。我想创建一个适合优化参数并获得从该分布中提取数字序列的可能性的对象。一种选择可能是只使用 KDE 模型并使用 pdf 来获得可能性。但是，我想要 2 个分布的确切均值和标准差。我想我可以，将数据分成两半，然后分别对 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
您似乎错误地指定了验配程序。
您试图在限制半边界的同时适应 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..
为了检查曲线下的面积，我使用了重新定义 mixture_model..:[=28=]
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 仍然有很大差异这是我无法重新创建的东西，即使在复制您的代码片段时也是如此。您是否可能不小心在某处重新定义了 functions/variables？