如何从numpy数组中确定概率分布函数？

Question

我在网上搜索了一下，没想到这个问题居然没有人回答。

我有一个包含10000个测量值的Numpy数组。我用Matplotlib画了一个直方图，从视觉上看这些值似乎是正态分布的：

不过，我想验证一下这个结果。我找到了一种正态性检验的方法，具体在scipy.stats.mstats.normaltest中实现，但结果却显示不是这样。我得到了这个输出：

(masked_array(data = [1472.8855375088663],
         mask = [False],
   fill_value = 1e+20)
, masked_array(data = [ 0.],
         mask = False,
   fill_value = 1e+20)

)

这意味着这个数据集是正态分布的可能性为0。我重新进行了实验并再次测试，结果还是一样，最好的情况下p值是3.0e-290。

我用以下代码测试了这个函数，似乎达到了我想要的效果：

import numpy
import scipy.stats as stats

mu, sigma = 0, 0.1
s = numpy.random.normal(mu, sigma, 10000)

print stats.normaltest(s)

(1.0491016699730547, 0.59182113002186942)

如果我理解和使用这个函数没错，那就意味着这些值并不是正态分布的。（说实话，我也不知道为什么输出会有差异，也就是细节少了。）

我原本很确定这是一个正态分布（虽然我对统计学的了解很基础），但我不知道还有什么其他的可能性。我该如何检查这个概率分布函数是什么呢？

编辑：

我的这个包含10000个值的Numpy数组是这样生成的（我知道这不是填充Numpy数组的最佳方式），然后运行了normaltest：

values = numpy.empty(shape=10000, 1))
for i in range(0, 10000):
    values[i] = measurement(...) # The function returns a float

print normaltest(values)

编辑 2：

我刚意识到输出之间的差异是因为我不小心用了两个不同的函数（scipy.stats.normaltest()和scipy.stats.mstats.normaltest()），但这并没有影响，因为无论用哪个函数，相关的输出部分都是一样的。

编辑 3：

根据askewchan的建议对直方图进行拟合：

plt.plot(bin_edges, scipy.stats.norm.pdf(bin_edges, loc=values.mean(), scale=values.std()))

结果是这样的：

编辑 4：

根据用户user333700的建议对直方图进行拟合：

scipy.stats.t.fit(data)

结果是这样的：

Answer 1

测试一个大样本是否符合某种分布通常比较复杂，因为如果样本和分布有任何偏差，测试就会把这些偏差当作异常值，从而拒绝这个分布。

所以我一般会使用QQ图来做这个。QQ图是一种图形工具，X轴表示数据的分位数，Y轴表示拟合分布的分位数。通过这种图形分析，可以选择对特定研究重要的分布部分：比如中心部分的分散情况，或者上下尾部。

为了做到这一点，我会使用DrawQQplot这个函数。

import openturns as ot
import numpy as np
sample = ot.Sample(s, 1)
tested_distribution = ot.NormalFactory().build(sample)
QQ_plot = ot.VisualTest.DrawQQplot(sample, tested_distribution)

这会生成以下图形。

QQ图验证了数据点是否在测试线附近。在当前情况下，拟合效果很好，不过我们注意到数据的极端分位数拟合得不太好（这也在意料之中，因为这些事件的概率密度很低）。

为了看看常见的情况，我尝试了BetaFactory，显然这是个错误的选择！

tested_distribution = ot.BetaFactory().build(sample)
QQ_plot = ot.VisualTest.DrawQQplot(sample, tested_distribution)

这会生成：

现在QQ图很明显：在中心区域的拟合是可以接受的，但对于低于-0.2或高于0.2的分位数就不能接受了。注意，Beta分布及其四个参数足够灵活，可以很好地拟合[0.2, 0.2]区间的数据。

如果样本量很大，我更倾向于使用核平滑而不是直方图。因为核平滑更准确，也就是更接近真实的、未知的概率密度函数（在AMISE误差方面，核平滑可以达到1/n^{4/5}，而直方图只能达到1/n^{2/3}），而且它是一个连续分布（你的分布看起来是连续的）。如果样本真的很大，还可以启用分箱，这样可以减少CPU的负担。

Answer 2

假设你正确地使用了测试，我猜测你可能有一个小的偏差，跟正常分布不太一样。而且因为你的样本量非常大，即使是小的偏差也会导致你拒绝“样本来自正常分布”的这个假设。

一种方法是通过绘制一个有很多区间的normed直方图来直观地检查你的数据，同时可以绘制概率密度函数（pdf），使用loc=data.mean()和scale=data.std()来设置。

还有其他测试正常性的方式，比如statsmodels提供的Anderson-Darling和Lillifors（Kolmogorov-Smirnov）测试，适用于当我们估计分布参数时。

不过，我认为结果不会有太大差别，因为样本量很大。

主要的问题是，你想测试你的样本是否“完全”来自正常分布，还是只是想知道你的样本是否来自一个与正常分布非常接近的分布，在实际应用中接近即可。

关于最后一点再详细说一下：

http://jpktd.blogspot.ca/2012/10/tost-statistically-significant.html http://www.graphpad.com/guides/prism/6/statistics/index.htm?testing_for_equivalence2.htm

随着样本量的增加，假设测试的能力会增强，这意味着即使是很小的差异，测试也能拒绝“相等”的假设。如果我们保持显著性水平不变，最终我们会拒绝一些其实并不重要的小差异。

另一种假设测试是我们想证明我们的样本接近某个特定的假设，比如两个样本的均值几乎相同。问题是我们需要定义什么是我们的等效区域。

在拟合优度测试中，我们需要选择一个距离度量，并定义样本与假设分布之间的距离阈值。我还没有找到任何解释，能帮助我们直观地选择这个距离阈值。

stats.normaltest是基于偏度和峰度与正常分布的偏差。

Anderson-Darling是基于累积分布函数（cdf）之间加权平方差的积分。

Kolmogorov-Smirnov是基于累积分布函数（cdf）之间的最大绝对差。

对于分箱数据的卡方检验则是基于加权的平方箱概率的总和。

等等。

我只尝试过对分箱或离散数据进行等效性测试，使用的阈值是来自一些参考案例的，这个阈值还是相对任意的。

在医学等效性测试中，有一些预定义的标准来说明何时可以认为两种治疗是等效的，或者在单侧版本中，何时可以认为一种治疗是劣于或优于另一种。