使用 scipy 高斯核密度估计计算 CDF 逆答案

【问题标题】：Using scipy gaussian kernel density estimation to calculate CDF inverse使用 scipy 高斯核密度估计计算 CDF 逆
【发布时间】：2017-11-21 16:28:03
【问题描述】：

scipy.stats中的gaussian_kde函数有一个函数evaluate可以返回一个输入点的PDF值。我正在尝试使用gaussian_kde 来估计逆CDF。其动机是生成一些输入数据的蒙特卡罗实现，这些输入数据的统计分布是使用 KDE 数值估计的。是否有绑定到gaussian_kde 的方法可以达到此目的？

下面的示例显示了这在高斯分布的情况下应该如何工作。首先，我展示了如何进行 PDF 计算以设置我想要实现的特定 API：

import numpy as np 
from scipy.stats import norm, gaussian_kde

npts_kde = int(5e3)
n = np.random.normal(loc=0, scale=1, size=npts_kde)
kde = gaussian_kde(n)

npts_sample = int(1e3)
x = np.linspace(-3, 3, npts_sample)
kde_pdf = kde.evaluate(x)
norm_pdf = norm.pdf(x)

是否有类似的简单方法来计算逆 CDF？ norm 函数有一个非常方便的 isf 函数，它正是这样做的：

cdf_value = np.sort(np.random.rand(npts_sample))
cdf_inv = norm.isf(1 - cdf_value)

kde_gaussian 是否存在这样的功能？还是直接从已经实现的方法构造这样的函数？

【问题讨论】：

如果您的最终目标是重新采样，为什么不使用resample 方法？
我的最终目标不是随机重采样，而是相关重采样。如果resample 方法允许我传递CDF 值（就像isf 方法所做的那样），那么我的问题就解决了。但是resample 假设我想使用随机制服通过逆 CDF 生成样本，我不这样做。
那么找根目录怎么样？会不会太慢了？
我从来没有用过这个kde 的东西，但是integrate_box_1d 方法对我来说听起来几乎是cdf，也许你甚至可以把-inf 作为一个边界？还有你可以使用根查找器反转的 cdf - 显然不是最快的。
看来我被抢了。无论如何，很高兴能帮上忙。

标签： python numpy scipy scientific-computing

【解决方案1】：

方法integrate_box_1d可以用来计算CDF，但不是向量化的；你需要循环点。如果内存不是问题，则以向量形式重写其source code（本质上只是对special.ndtr 的调用）可能会加快速度。

from scipy.special import ndtr
stdev = np.sqrt(kde.covariance)[0, 0]
pde_cdf = ndtr(np.subtract.outer(x, n)).mean(axis=1)
plot(x, pde_cdf)

反函数的绘图将是plot(pde_cdf, x)。如果目标是计算特定点的反函数，请考虑使用the inverse of interpolating spline，对 CDF 的计算值进行插值。

【讨论】：

我发现我不得不稍微修改 pde_cdf 行：pde_cdf = ndtr(np.subtract.outer(x, n)/stdev).mean(axis=1)。您将在指向的源代码中看到除以标准差。从数学上讲，我认为无论如何这是必需的。如果你有幸看到正态分布之类的东西，这种方法就可以了。但是如果没有stdev，如果您正在查看任何不“正常”的东西，它真的会崩溃。它在这里工作是因为你的高斯宽度是 1.0。

【解决方案2】：

您可以使用一些 Python 技巧来快速且有效地估计 CDF（基于 this answer）：

    from scipy.special import ndtr
    cdf = tuple(ndtr(np.ravel(item - kde.dataset) / kde.factor).mean()
                for item in x)

它的运行速度与 this answer 一样快，但具有线性 (len(kde.dataset)) 空间复杂度，而不是二次（实际上是 len(kde.dataset) * len(x)）空间复杂度。

接下来您要做的就是使用逆近似，例如，来自statsmodels。

【讨论】：