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

aph*_*aph 7 python numpy scientific-computing scipy

gaussian_kde函数scipy.stats具有如下功能: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)
Run Code Online (Sandbox Code Playgroud)

正态分布 PDF 的 KDE 近似演示

是否有类似的简单方法来计算逆 CDF?该norm函数有一个非常方便的isf函数,可以做到这一点:

cdf_value = np.sort(np.random.rand(npts_sample))
cdf_inv = norm.isf(1 - cdf_value)
Run Code Online (Sandbox Code Playgroud)

正态分布 CDF 的所需 KDE 近似演示

是否存在这样的函数kde_gaussian?或者从已经实现的方法构造这样的函数是否很简单?

小智 4

该方法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)
Run Code Online (Sandbox Code Playgroud)

反函数的图为plot(pde_cdf, x)。如果目标是计算特定点处的反函数,请考虑使用插值样条函数的反函数,对 CDF 的计算值进行插值。

  • 我发现我必须稍微修改 pde_cdf 行:`pde_cdf = ndtr(np.subtract.outer(x, n)/stdev).mean(axis=1)`。您将在指向的源代码中看到除以标准差。从数学上来说,我认为无论如何这是必需的。如果你足够幸运,能够观察正态分布之类的东西,那么这种情况就可以了。但是如果没有“stdev”,如果您正在查看任何不“正常”的东西,它真的会崩溃。它在这里起作用是因为你的高斯宽度是 1.0。 (2认同)