Python中的密度估计答案

【问题标题】：Density estimation in PythonPython中的密度估计
【发布时间】：2017-05-25 07:59:32
【问题描述】：

我是 Python 新手，我必须从 2D 样本中估计密度。我的第一个想法是一个简单的直方图估计器，我实现如下：

num = 10**4
sp = 0.01

subsetf1 = np.random.uniform(0,1,size=(num,2)) # I created this set to test the estimator

def f_est(x,y,h=sp, subset=subsetf1):
    indicator = np.zeros(num)

    for i in range(num):
        if (x <= subset[i][0] <= (x + h)) and (y <= subset[i][1] <= (y + h)):
            indicator[i] = 1
        else :  indicator[i] = 0


    return sum(indicator)/(num*h**2)

#f_est should yield values closely to 1 if  0 <x,y <1 , because  subsetf1 contains uniformly distributed values on [0,1).

问题是f_est 经常产生大于 1 的值，这让我相信我的代码是错误的，但我不知道错误可能在哪里。我也尝试过内核密度估计，但如果我尝试类似：

from scipy import stats
xmin = partsetf1[0].min()
xmax = partsetf1[0].max()
ymin = partsetf1[1].min()
ymax = partsetf1[1].max()
X, Y = np.mgrid[xmin:xmax:100j, ymin:ymax:100j]
positions = np.vstack([X.ravel(), Y.ravel()])
values = np.vstack([partsetf1[0], partsetf1[1]])
gkde=stats.gaussian_kde(values)
f = np.reshape(gkde(positions).T, X.shape)

f 会产生不正确的奇怪值。如果有人能告诉我代码中的错误在哪里，以及如何在 Python 中使用 kde，那就太好了，因为我没有找到关于这个主题的好的教程。

【问题讨论】：

您是否将结果与 scipy 的密度估计等内置函数进行了比较？ docs.scipy.org/doc/scipy/reference/generated/… 或者更确切地说，我应该问，如果有的话，您将结果与什么进行了比较？
嗨@Jon。我还尝试了 scipy 的密度估计，因此代码也在我的帖子中。问题是 scipys 密度估计会产生完全错误的结果，我不知道为什么
对不起，我没看到。 Scipy 的替代方案是 statsmodels。你试过这个statsmodels.sourceforge.net/notebooks/generated/…

标签： python statistics kernel-density

【解决方案1】：

核密度估计可以产生大于 1 的值，因为它们返回的不是概率，而是概率密度。在您的情况下，概率密度给出了每单位面积的概率，并且可以在局部远大于 1。

好消息是，由于您的结果与 scipy 或 statsmodels 的结果相匹配，听起来您的代码是正确的。

【讨论】：