使用scipy高斯核密度估计计算CDF反函数

Question

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 是否存在这样的函数？或者是否可以通过已经实现的方法简单构造出这样的函数？

- aph

---

如果您的最终目标是重新采样，为什么不使用“resample”方法呢？ - Paul Panzer

那么根查找器呢？会不会太慢？ - Paul Panzer

我从未使用过这个“kde”东西，但“integrate_box_1d”方法听起来对我来说几乎像是累积分布函数，也许你甚至可以将“-inf”作为边界？而且您可以使用根查找器反转cdf-显然不是最快的方法。 - Paul Panzer

是的，这是一个非常好的建议。使用scipy.integrate.quad进行PDF的数值积分比“integrate_box_1d”慢了约50倍。因此，这实际上会非常快。 如果您将其编写为建议的答案，那么很可能会被接受。否则，在明确清晰地解释之后，我会将其编写出来。无论哪种方式，感谢您的提示！ - aph

看起来我被抢先了。不管怎样，很高兴能够帮到你。 - Paul Panzer

Answer 1

方法integrate_box_1d可用于计算CDF，但它不是向量化的；您需要循环处理点。如果内存不是问题，将其源代码（本质上只是对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)。如果目标是在特定点计算反函数，请考虑使用插值样条的反函数，对CDF的计算值进行插值。

- user6655984

---

我发现我必须稍微修改pde_cdf这一行：pde_cdf = ndtr(np.subtract.outer(x, n)/stdev).mean(axis=1)。您会在您指向的源代码中看到除以标准差的操作。从数学上讲，我认为这是必需的。如果你足够幸运能够看到像正态分布这样的东西，这样做还算可以。但是如果没有“stdev”，它在查看任何不是“正常”的东西时真的会出问题。这里之所以有效，是因为高斯函数的宽度为1.0。 - Gordon

Answer 2

其他答案已经回答了这个问题，但我花了一些时间来理清思路。以下是最终解决方案的完整示例：

import numpy as np 
from scipy import interpolate
from scipy.special import ndtr
import matplotlib.pyplot as plt
from scipy.stats import norm, gaussian_kde

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

# grid for plotting
npts_sample = int(1e3)
x = np.linspace(-3, 3, npts_sample)

# evaluate pdfs
kde_pdf = kde.evaluate(x)
norm_pdf = norm.pdf(x)

# cdf and inv cdf are available directly from scipy
norm_cdf = norm.cdf(x)
norm_inv = norm.ppf(x)

# estimate cdf
cdf = tuple(ndtr(np.ravel(item - kde.dataset) / kde.factor).mean()
            for item in x)

# estimate inv cdf
inversefunction = interpolate.interp1d(cdf, x, kind='cubic', bounds_error=False)

fig, ax = plt.subplots(1, 3, figsize=(6, 3))
ax[0].plot(x, norm_pdf, c='k')
ax[0].plot(x, kde_pdf, c='r', ls='--')
ax[0].set_title('PDF')
ax[1].plot(x, norm_cdf, c='k')
ax[1].plot(x, cdf, c='r', ls='--')
ax[1].set_title('CDF')
ax[2].plot(x, norm_inv, c='k')
ax[2].plot(x, inversefunction(x), c='r', ls='--')
ax[2].set_title("Inverse CDF")

- kilojoules

Answer 3

你可以使用一些Python技巧来快速和节省内存地估计CDF（基于这个答案）：

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

它的速度与这个答案一样快，但其空间复杂度是线性的(len(kde.dataset))，而不是二次的(实际上，是len(kde.dataset) * len(x))。

接下来要做的就是使用逆近似方法，例如从statsmodels中获取。

- Dmitry