Python中大量数据的高斯核密度估计（KDE）

Question

我有1000个大数，随机分布在37231到56661的范围内。

我试图使用stats.gaussian_kde，但是似乎出了些问题（也许是因为我对统计学知识了解不够？）。

以下是代码：

from scipy import stats.gaussian_kde
import matplotlib.pyplot as plt

# 'data' is a 1D array that contains the initial numbers 37231 to 56661
xmin = min(data)
xmax = max(data)   

# get evenly distributed numbers for X axis.
x = linspace(xmin, xmax, 1000)   # get 1000 points on x axis
nPoints = len(x)

# get actual kernel density.
density = gaussian_kde(data)
y = density(x)

# print the output data
for i in range(nPoints):
    print "%s   %s" % (x[i], y[i])

plt.plot(x, density(x))
plt.show()

在打印输出中，第一列中得到了x值，在第二列中得到了零值。图表显示一条平直的线。

我无法找到解决方案。 我尝试了非常广泛的X值，但结果相同。

问题出在哪里？我做错了什么？ 大数是否可能是原因？

- Proteos

---

请注意顶部的格式错误；您可以选择所有代码，然后点击“{}”按钮，在每行前添加必要的四个空格。 - sarnold

@sarnold， 对不起，你指的是哪个错误？我实际上使用了那个{}按钮，在我的Mac上格式看起来很好。（我是一个新手，在此提前为错误道歉） - Proteos

@Proteos：看第一行，以“from scipy import…”开头。它没有标记为代码。 - DSM

哈！这下我可吃了闭门羹，竟然在没看源代码的情况下给出了建议；在源代码之前必须留出一行空白。虽然有点傻，但你是对的，代码全都在那里... - sarnold

Answer 1

我认为发生的情况是您的数据数组由整数组成，这会导致问题：

>>> import numpy, scipy.stats
>>> 
>>> data = numpy.random.randint(37231, 56661,size=10)
>>> xmin, xmax = min(data), max(data)
>>> x = numpy.linspace(xmin, xmax, 10)
>>> 
>>> density = scipy.stats.gaussian_kde(data)
>>> density.dataset
array([[52605, 45451, 46029, 40379, 48885, 41262, 39248, 38247, 55987,
        44019]])
>>> density(x)
array([0, 0, 0, 0, 0, 0, 0, 0, 0, 0])

但是如果我们使用浮动：

>>> density = scipy.stats.gaussian_kde(data*1.0)
>>> density.dataset
array([[ 52605.,  45451.,  46029.,  40379.,  48885.,  41262.,  39248.,
         38247.,  55987.,  44019.]])
>>> density(x)
array([  4.42201513e-05,   5.51130237e-05,   5.94470211e-05,
         5.78485526e-05,   5.21379448e-05,   4.43176188e-05,
         3.66725694e-05,   3.06297511e-05,   2.56191024e-05,
         2.01305127e-05])

- DSM

---

哦，多么天真的错误！ 我以为我缺少一些简单的东西，但是那么简单吗？ :)另一方面，gaussian_kde()函数应该负责转换为浮点数；至少给出需要浮点数的警告。 你不同意吗？好了，问题解决了！ 非常感谢你！ - Proteos

@Proteos：我同意，这看起来像是一个bug。 - DSM

不错的发现。我可以看到争论可以朝任何一方倾斜：结果返回的小数位数与您的输入数据一样多，而且 1e-05 在许多情况下足够接近 0，特别是当输入远离 0 时。尽管如此，这绝对是令人惊讶的。 - sarnold

这个问题已经在scipy中得到了修复，并将在下一个版本（scipy 0.10）中转换为浮点数。 - Josef

Answer 2

我已经编写了一个函数来实现这一点。您可以将带宽作为函数的参数进行调整。换句话说，较小的数值会使结果更尖锐，较大的数值则会使结果更平滑。默认值为0.3。

此函数在IPython笔记本--pylab=inline中可行。

箱数已经被优化和编码，因此将根据您数据中变量的数量而变化。

import scipy.stats as stats
import matplotlib.pyplot as plt
import numpy as np

def hist_with_kde(data, bandwidth = 0.3):
    #set number of bins using Freedman and Diaconis
    q1 = np.percentile(data,25)
    q3 = np.percentile(data,75)


    n = len(data)**(.1/.3)
    rng = max(data) - min(data)
    iqr = 2*(q3-q1)
    bins = int((n*rng)/iqr)

    x = np.linspace(min(data),max(data),200)

    kde = stats.gaussian_kde(data)
    kde.covariance_factor = lambda : bandwidth
    kde._compute_covariance()

    plt.plot(x,kde(x),'r') # distribution function
    plt.hist(data,bins=bins,normed=True) # histogram

data = np.random.randn(500)
hist_with_kde(data,0.25)

- John