scipy.stats.kde 和 scipy.stats.kstest

Question

我该如何以符合规范的方式使用 scipy.stats.kde.gaussian_kde 和 scipy.stats.kstest 呢？

比如，下面这段代码：

from numpy import inf
import scipy.stat
my_pdf = scipy.stats.kde.gaussian_kde(sample)
scipy.stats.kstest(sample, lambda x: my_pdf.integrate_box_1d(-inf, x))

得到了这样的结果： (0.5396735893479544, 0.0)

这个结果是不正确的，因为一个样本显然是属于基于这个样本构建的分布的。

Answer 1

首先，如果你想测试两个样本是否可能来自同一个分布，最合适的测试方法是双样本KS检验，这个方法在scipy.stats.ks_2samp中实现，它直接比较经验累积分布函数（CDF）。而KDE是密度估计，它会对CDF进行平滑处理，因此这其实是多余的工作，反而会让你的估计在统计上变得更糟。

你遇到这个问题的原因是你的CDF参数的签名不太对。kstest会调用cdf(vals) (源代码)，其中vals是排序后的样本，用来获取每个样本的CDF值。在你的代码中，这最终调用了my_pdf.integrate_box_1d(-np.inf, samps)，但integrate_box_1d希望两个参数都是单个数值（标量）。所以签名不对，如果你用大多数数组来试这个，会出现ValueError错误：

>>> my_pdf.integrate_box_1d(-np.inf, samp[:10])
---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
<ipython-input-38-81d0253a33bf> in <module>()
----> 1 my_pdf.integrate_box_1d(-np.inf, samp[:10])

/Library/Python/2.7/site-packages/scipy-0.12.0.dev_ddd617d_20120725-py2.7-macosx-10.8-x86_64.egg/scipy/stats/kde.pyc in integrate_box_1d(self, low, high)
    311 
    312         normalized_low = ravel((low - self.dataset) / stdev)
--> 313         normalized_high = ravel((high - self.dataset) / stdev)
    314 
    315         value = np.mean(special.ndtr(normalized_high) - \

ValueError: operands could not be broadcast together with shapes (10) (1,1000) 

但不幸的是，当第二个参数是samp时，它可以正常广播，因为这两个数组的形状是一样的，然后一切就乱套了。可以推测integrate_box_1d应该检查它的参数形状，但这里有一种正确的方法：

>>> my_cdf = lambda ary: np.array([my_pdf.integrate_box_1d(-np.inf, x) for x in ary])
>>> scipy.stats.kstest(sample, my_cdf)
(0.015597917205996903, 0.96809912578616597)

如果你愿意的话，也可以使用np.vectorize。

（不过再说一次，你可能实际上更想用ks_2samp。）