Python中大数的高斯核密度估计（KDE）答案

【问题标题】：Gaussian Kernel Density Estimation (KDE) of large numbers in PythonPython中大数的高斯核密度估计（KDE）
【发布时间】：2012-04-06 13:13:16
【问题描述】：

我有 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()

在打印输出中，我在第 1 列中得到 x 值，在第 2 列中得到零。该图显示一条平线。

我根本找不到解决方案。我尝试了非常广泛的 X-es，结果相同。

有什么问题？我究竟做错了什么？大数字可能是原因吗？

【问题讨论】：

注意顶部附近的格式错误；您可以选择所有代码并点击{} 按钮在每行之前添加必要的四个空格。
@sarnold，对不起，您指的是哪个错误？我实际上使用了那个 {} 按钮，至少在我的 Mac 上格式看起来不错。（我是这里的新手，但我提前为这个错误道歉）
@Proteos：看第一行，从“from scipy import..”开始。它没有标记为代码。
哈！这将教会我在不查看源代码的情况下提供建议；您必须在源代码之前留一个空行。这有点傻，但你是对的，代码都在那里......

标签： python statistics matplotlib scipy

【解决方案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])

【讨论】：

哦，多么天真的错误！我以为我错过了一些简单的东西，但就这么简单吗？ :) 另一方面， gaussian_kde() 函数应该负责转换为浮点数；至少给出它需要浮动的警告。你不同意吗？好了，案子解决了！非常感谢！
@Proteos：我同意，这对我来说似乎是个错误。
不错的收获。我可以看到这个论点是两种方式：返回的结果与您的输入数据一样多，小数点后的数字，1e-05 足够接近0 用于许多用途，特别是当输入远离@ 987654325@。不过，这绝对令人惊讶。
这已在 scipy 中修复，并将在下一个版本（scipy 0.10）中转换为浮点数

【解决方案2】：

我已经做了一个函数来做到这一点。您可以改变带宽作为函数的参数。也就是说，较小的数字 = 更尖，更大的数字 = 更平滑。默认值为 0.3。

它适用于IPython notebook --pylab=inline

bin 的数量经过优化和编码，因此会因数据中的变量数量而异。

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)

【讨论】：