Scipy文档尝试翻译之二：scipy.stats scipy统计

统计（scipy.stats）

简介

在本教程中，我们讨论scipy.stats的许多——当然不是全部——功能 。我们的目的是向用户提供此软件包的实用知识。有关更多详细信息，请参阅参考手册。

注意：本文档正在完善中，若有不完善之处请与作者联系。

离散统计分布
连续统计分布

随机变量

scipy.stats有两种通用分布类，分别用于封装连续随机分布和离散随机分布。

使用这些类已经实现了80多个连续随机分布（RV）和10个离散随机分布。除此之外，终端用户可以轻松添加新的例程和发行版。（如果创建一个，请贡献它。）

所有统计功能都位于子程序包中 scipy.stats，使用命令info(stats)可以获得这些功能的清单。可用的随机分布列表也可以从stats子软件包的文档中获取。

在下面的讨论中，我们主要关注连续性随机变量。当然了，以下的几乎所有内容都适用于离散变量，但是我们在这里指出了一些区别：离散分布的特殊之处

在下面的代码示例中，我们假设scipy.stats程序包已经按照以下方式导入：

>>> from scipy import stats

并且在某些情况下，我们假设将所需的单个分布（比如以下代码中的正态分布norm）以以下形式导入：

>>> from scipy.stats import norm

为了使Python 2和Python 3保持一致，我们还将确保print是以函数的形式导入的。
<译者注：python2需要执行以下代码以避免语法错误；如果使用的是python3则无需这一步。>

>>> from __future__ import print_function

获取帮助

首先，所有分布都带有帮助函数。为了仅获取一些基本信息，我们将打印相关文档字符串：print(stats.norm.__doc__)。

要找到支持，比如，分布的上下限，请调用：

>>> print('bounds of distribution lower: %s, upper: %s' % (norm.a, norm.b))
bounds of distribution lower: -inf, upper: inf

我们可以使用dir(norm)列出这一分布的所有方法和属性 。其中，有些方法是私有的，尽管它们没有以私有的形式命名（亦即名称不以下划线开头，比如veccdf）。这些私有的方法和变量只能用于内部计算，如果尝试外部调用，则会发出警告。

为了获得真正的主要方法，在这里列出了冻结分布的方法。冻结分布的含义将在下面解释。

>>> rv = norm()
>>> dir(rv)  # reformatted
['__class__', '__delattr__', '__dict__', '__dir__', '__doc__', '__eq__',
 '__format__', '__ge__', '__getattribute__', '__gt__', '__hash__',
 '__init__', '__le__', '__lt__', '__module__', '__ne__', '__new__',
 '__reduce__', '__reduce_ex__', '__repr__', '__setattr__', '__sizeof__',
 '__str__', '__subclasshook__', '__weakref__', 'a', 'args', 'b', 'cdf',
 'dist', 'entropy', 'expect', 'interval', 'isf', 'kwds', 'logcdf',
 'logpdf', 'logpmf', 'logsf', 'mean', 'median', 'moment', 'pdf', 'pmf',
 'ppf', 'random_state', 'rvs', 'sf', 'stats', 'std', 'var']

最后，我们可以通过自省（调用函数dir()）获得可用的分布的列表：

>>> dist_continu = [d for d in dir(stats) if
...                 isinstance(getattr(stats, d), stats.rv_continuous)]
>>> dist_discrete = [d for d in dir(stats) if
...                  isinstance(getattr(stats, d), stats.rv_discrete)]
>>> print('number of continuous distributions: %d' % len(dist_continu))
number of continuous distributions: 100
>>> print('number of discrete distributions:   %d' % len(dist_discrete))
number of discrete distributions:   15

求分布的特征函数

连续随机变量

连续随机变量的主要公共方法是：

rvs：随机变量

pdf：概率密度函数

cdf：累积分布函数

sf：生存函数（1-CDF）<如果是可靠性行业，就是t时刻的可靠度R(t)>

ppf：百分比点函数（CDF的反函数）

isf：逆生存函数（SF的反函数）

stats：返回均值，方差，（费舍尔）偏度或（费舍尔）峰度

moment：分布的非中心矩

让我们以正态随机变量为例。

>>> norm.cdf(0)
0.5

要计算cdf多个点的值，我们可以传递一个列表或一个numpy数组。

>>> norm.cdf([-1., 0, 1])
array([ 0.15865525,  0.5,  0.84134475])
>>> import numpy as np
>>> norm.cdf(np.array([-1., 0, 1]))
array([ 0.15865525,  0.5,  0.84134475])

因此，分布的一些基本方法（例如pdf，cdf等）都有矢量化处理。

一些其他通用的、好用的方法也是支持的：

>>> norm.mean(), norm.std(), norm.var()
(0.0, 1.0, 1.0)
>>> norm.stats(moments="mv")
(array(0.0), array(1.0))

要找到分布的中位数，我们可以使用百分比点函数ppf，它是cdf的反函数：

>>> norm.ppf(0.5)
0.0

要生成随机变量序列(类型为np.ndarray)，请使用size关键字参数：

>>> norm.rvs(size=3)
array([-0.35687759,  1.34347647, -0.11710531])   # random

请注意，随机数的生成依赖于numpy.random 包中的生成器。因此如果你重复输入上面的代码，会得到不同的结果。为了实现可重复性，您可以显式地设置全局随机种子：

>>> np.random.seed(1234)

但是，不建议依赖全局状态。更好的方法是使用随机种子中的random_state参数，该参数接受numpy.random.mtrand.RandomStateclass 的实例或一个整数，然后将其用作内部RandomState对象的种子，这样的话结果就是可重现的：

>>> norm.rvs(size=5, random_state=1234)
array([ 0.47143516, -1.19097569,  1.43270697, -0.3126519 , -0.72058873])

当然，不要认为norm.rvs(5)会产生长度为5的一个array：

>>> norm.rvs(5)
5.471435163732493

在这里，5作为无关键字参数，被解释为第一个可能的关键字参数loc。这是所有连续分布采用的一对关键字参数中的第一个。
如果对这个问题感到困惑，不用担心，我们下一部分的主题将详细讲解。

分布的移动和缩放

所有连续分布均以loc和scale作为关键字参数来调整分布的位置和比例，例如，对于标准正态分布，位置是平均值，比例是标准偏差。

>>> norm.stats(loc=3, scale=4, moments="mv")
(array(3.0), array(16.0))

在许多情况下，X 可通过以下变换获得随机变量的标准化分布：
$$\frac{x-loc}{scale}$$
式中loc为位置参数，scale为比例参数。默认值为loc = 0和scale = 1。

巧妙地使用loc和scale可以通过多种方式帮助修改标准分布。为了进一步说明缩放比例， cdf具有均值$\frac{1}{\lambda }$的指数分布随机变量可以使用下式定义：
$$F(x)=1-e^{-\lambda x}$$

通过应用上面的缩放规则，可以看出通过采用scale = 1/lambda，我们就可以得到适当的缩放比例。

>>> from scipy.stats import expon
>>> expon.mean(scale=3.)
3.0

注意<翻译不全>
要获得理想的分布形状，分布可能不仅仅需要简单应用loc和/或 scale。例如，给定长度R的矢量，服从正态分布$N（0,{\sigma }^2）$，则二维矢量长度的分布 受独立N（0， ）每个成分的偏差是米，scale = ）。第一个参数是一个形状参数，需要与x一起调整。

均匀分布也很有趣：

>>> from scipy.stats import uniform
>>> uniform.cdf([0, 1, 2, 3, 4, 5], loc=1, scale=4)
array([ 0.  ,  0.  ,  0.25,  0.5 ,  0.75,  1.  ])

最后，回想一下上一段，我们还剩下norm.rvs(5)的含义问题。事实证明，调用这样的分布，第一个参数（即5）被传递以设置loc参数。让我们来看看：

>>> np.mean(norm.rvs(5, size=500))
5.0098355106969992

因此，为了解释最后一部分示例的输出： 由于default，norm.rvs(5)产生了一个服从正态分布随机变量。它的loc=5，这是由于默认状态下，size=1。连续生成500个之后取均值也是在5附近的。

我们建议您通过将值作为关键字参数，而不是通过参数传递来设置loc和scale参数。如下所述，当使用冻结分布的技术调用给定随机变量分布的不止一种方法时，可以尽量避免重复。

分布的形状参数

虽然一般的连续随机变量可以使用loc和scale参数进行移动和缩放，但某些分布还需要其他形状参数。例如，具有密度的伽玛分布
$$\gamma (x,a)=\frac{\lambda (\lambda x{)}^{\alpha -1}}{\Gamma (a)}{e}^{-\lambda x}$$
需要形状参数a。我们发现，设置$\lambda$可以依靠设置形状参数$\frac{1}{\lambda }$来实现。
观察该设置 可以通过将scale关键字 设置为获得。

让我们检查伽玛分布的形状参数的数量和名称。（我们从上面知道应该为1。）
<此处翻译存疑>

>>>
>>> from scipy.stats import gamma
>>> gamma.numargs
1
>>> gamma.shapes
'a'

现在，我们将shape变量的值设置为1以获得指数分布，这样我们就可以轻松地比较是否获得期望的结果。

>>> gamma(1, scale=2.).stats(moments="mv")
(array(2.0), array(4.0))

请注意，我们还可以将形状参数a指定为关键字：

>>> gamma(a=1, scale=2.).stats(moments="mv")
(array(2.0), array(4.0))

冻结分布
反复传递loc和scale关键字可能会变得很麻烦。冻结随机变量的概念用于解决此类问题。

>>> rv = gamma(1, scale=2.)

通过使用经过指定参数的随机变量rv，我们不再需要包含比例或形状参数。因此，可以通过以下两种方式之一使用分发：通过将所有分发参数传递给每个方法调用（如我们之前所做的那样），或者通过冻结分发实例的参数。让我们检查一下：

>>> rv.mean(), rv.std()
(2.0, 2.0)

的确，这就是我们应该得到的。

分布的广播规则（numpy）

基本方法pdf等满足通常的numpy广播规则。例如，我们可以针对不同的概率和自由度计算t分布的上尾的临界值。

>>> stats.t.isf([0.1, 0.05, 0.01], [[10], [11]])
array([[ 1.37218364,  1.81246112,  2.76376946],
       [ 1.36343032,  1.79588482,  2.71807918]])

在此，第一行包含10个自由度的临界值，第二行包含11个自由度（dof）的临界值。因此，广播规则给出isf两次调用的相同结果：

>>> stats.t.isf([0.1, 0.05, 0.01], 10)
array([ 1.37218364,  1.81246112,  2.76376946])
>>> stats.t.isf([0.1, 0.05, 0.01], 11)
array([ 1.36343032,  1.79588482,  2.71807918])

如果具有概率的数组$[0.1,0.05,0.01]$和自由度数组$[10,11,12]$具有相同的数组形状，则使用按元素匹配。
例如，我们可以通过调用获得10自由度的10％尾巴，11自由度的5％尾巴和12自由度的1％尾巴。[0.1, 0.05, 0.01][10, 11, 12]

>>> stats.t.isf([0.1, 0.05, 0.01], [10, 11, 12])
array([ 1.37218364,  1.79588482,  2.68099799])

离散分布的特定点

离散分布与连续分布具有基本相同的基本方法。但是pdf应当由由概率质量函数pmf代替，并且没有可用的估算方法（例如拟合），并且scale不是有效的关键字参数。位置参数keyword loc仍可用于移动分布。

CDF的计算需要额外注意。在连续分布的情况下，在大多数标准情况下，累积分布函数在范围（a，b）中严格单调递增，因此具有唯一的逆。但是，离散分布的cdf是阶跃函数，因此逆cdf（即百分比点函数）需要不同的定义：

ppf(q) = min{x : cdf(x) >= q, x integer}
有关更多信息，请参见此处的文档。

我们可以以超几何分布为例

>>> from scipy.stats import hypergeom
>>> [M, n, N] = [20, 7, 12]

如果我们在某些整数点使用cdf，然后在这些cdf值上评估ppf，则可以返回初始整数，例如

>>> x = np.arange(4)*2
>>> x
array([0, 2, 4, 6])
>>> prb = hypergeom.cdf(x, M, n, N)
>>> prb
array([  1.03199174e-04,   5.21155831e-02,   6.08359133e-01,
         9.89783282e-01])
>>> hypergeom.ppf(prb, M, n, N)
array([ 0.,  2.,  4.,  6.])

但是，如果我们使用的值不属于cdf步函数，那么我们将获得下一个更大的整数：？？？<翻译存疑。原文：If we use values that are not at the kinks of the cdf step function, we get the next higher integer back:>

>>> hypergeom.ppf(prb + 1e-8, M, n, N)
array([ 1.,  3.,  5.,  7.])
>>> hypergeom.ppf(prb - 1e-8, M, n, N)
array([ 0.,  2.,  4.,  6.])

拟合分布

非冻结分布的主要其他方法与分布参数的估计有关：

拟合：分布参数（包括位置）的最大似然估计
和尺度

fit_loc_scale：给出形状参数时位置和比例的估计

nnlf：负对数似然函数

期望：根据pdf或pmf计算函数的期望
<文档中缺少这方面的例子，在后面有相应的调用（t检验附近），译者注>

性能问题和注意事项

就速度而言，各个方法的性能会因分布和方法而有很大差异。方法的结果可以通过以下两种方式之一获得：通过显式计算或通过独立于特定分布的通用算法。

一方面，显式计算要求通过分析公式或in scipy.special或 numpy.randomfor中的特殊函数直接为给定分布指定方法rvs。这些通常是相对较快的计算。

另一方面，如果分布没有指定任何显式计算，则使用通用方法。要定义分布，只需要pdf或cdf之一；所有其他方法都可以使用数值积分和求根法导出。但是，这些间接方法可能非常慢。例如，以一种非常间接的方式创建随机变量，如代码rgh = stats.gausshyper.rvs(0.5, 2, 2, 2, size=100)，在我的计算机上花费100个随机变量大约需要19秒，而来自标准正态或t分布的一百万个随机变量仅需一秒钟多一点。

剩下的问题

scipy.stats最近对分布进行了更正和改进，并获得了大量测试套件；但是，仍然存在一些问题：

分布已在某些参数范围内进行了测试；但是，在某些转角范围内，可能会保留一些不正确的结果。

拟合中的最大似然估计不适用于所有分布的默认启动参数，并且用户需要提供良好的启动参数。同样，对于某些使用最大似然估计器的分布，固有地可能不是最佳选择。

建立自定义的具体分布

下一个示例显示了如何构建自己的分布。其他示例显示了分布的用法和一些统计检验。

通过子类化rv_continuous产生连续分布

产生连续分布非常简单。

>>> from scipy import stats
>>> class deterministic_gen(stats.rv_continuous):
...     def _cdf(self, x):
...         return np.where(x < 0, 0., 1.)
...     def _stats(self):
...         return 0., 0., 0., 0.
>>> deterministic = deterministic_gen(name="deterministic")
>>> deterministic.cdf(np.arange(-3, 3, 0.5))
array([ 0.,  0.,  0.,  0.,  0.,  0.,  1.,  1.,  1.,  1.,  1.,  1.])

有趣的是，这个示例中，pdf计算函数无需定义，便可以自动计算：

>>> deterministic.pdf(np.arange(-3, 3, 0.5))
array([  0.00000000e+00,   0.00000000e+00,   0.00000000e+00,
         0.00000000e+00,   0.00000000e+00,   0.00000000e+00,
         5.83333333e+04,   4.16333634e-12,   4.16333634e-12,
         4.16333634e-12,   4.16333634e-12,   4.16333634e-12])

请注意性能问题和注意事项中提到的性能问题。未指定的通用方法的计算可能会变得非常缓慢，因为仅调用了通用方法，从本质上讲，它们无法使用有关分布的任何特定信息。
因此，以下示例可以作为一个警告：

>>> from scipy.integrate import quad
>>> quad(deterministic.pdf, -1e-1, 1e-1)
(4.163336342344337e-13, 0.0)

但这是不正确的：此pdf上的积分应为1。让我们减小积分间隔：

>>> quad(deterministic.pdf, -1e-3, 1e-3)  # warning removed
(1.000076872229173, 0.0010625571718182458)

这看起来更好。但是，问题的根源是pdf方法未能在分布类中给出其独特的计算方法，因此效率较低。
<译者注：若未定义，那么pdf是按照通用方法算出的，也就是对cdf求导。这种方法显然比自定义一个pdf方法按照pdf的公式求值要慢得多。又由于上面的分布函数有一个阶跃，在0处是间断的，所以区间划分太小的话就无法求出导数>

子类rv_discrete生成离散分布

在下文中，我们将stats.rv_discrete用于生成离散分布，该离散分布对于以整数为中心的间隔具有截断正态的概率。

基本信息

在rv_discrete的文档字符串中help(stats.rv_discrete)，

“您可以通过传递rv_discrete初始化方法（通过values =关键字）来构造任意离散的rv，其中 P { X = x k } = p k P \{X = x_k\} = p_k P{X=xk​}=pk​的序列元组$（x_k，p_k）$仅描述X的那些值$x_k$以非零概率$p_k$发生。”

除此之外，此方法还需要满足一些其他要求：

关键字名称为必填项。
分布xk的支撑点必须是整数。<?>
需要指定有效位数（小数）。

实际上，如果不满足最后两个要求，则可能会引发异常，或者所产生的数字可能不正确。

一个例子

让我们开干吧！首先按照下面的办法做：

>>> npoints = 20   # 整数数据点的数目为20个。number of integer support points of the distribution minus 1
>>> npointsh = npoints // 2
>>> npointsf = float(npoints)
>>> nbound = 4   # 正态分布截尾的范围 bounds for the truncated normal
>>> normbound = (1+1/npointsf) * nbound   # 正态分布的实际分布范围 actual bounds of truncated normal
>>> grid = np.arange(-npointsh, npointsh+2, 1)   # integer grid
>>> gridlimitsnorm = (grid-0.5) / npointsh * nbound   # bin limits for the truncnorm
>>> gridlimits = grid - 0.5   # 在之后的分析中用到
>>> grid = grid[:-1]
>>> probs = np.diff(stats.truncnorm.cdf(gridlimitsnorm, -normbound, normbound))
>>> gridint = grid

最后，我们可以继承rv_discrete：

>>> normdiscrete = stats.rv_discrete(values=(gridint,
...              np.round(probs, decimals=7)), name='normdiscrete')
现在我们已经定义了分布，我们可以访问所有常见的离散分布方法。

>>> print('mean = %6.4f, variance = %6.4f, skew = %6.4f, kurtosis = %6.4f' %
...       normdiscrete.stats(moments='mvsk'))
mean = -0.0000, variance = 6.3302, skew = 0.0000, kurtosis = -0.0076

>>> nd_std = np.sqrt(normdiscrete.stats(moments='v'))

测试以上的结果

让我们生成一个随机样本，并将观察到的频率与概率进行比较。

>>> n_sample = 500
>>> np.random.seed(87655678)   # fix the seed for replicability
>>> rvs = normdiscrete.rvs(size=n_sample)
>>> f, l = np.histogram(rvs, bins=gridlimits)
>>> sfreq = np.vstack([gridint, f, probs*n_sample]).T
>>> print(sfreq)
[[-1.00000000e+01  0.00000000e+00  2.95019349e-02]
 [-9.00000000e+00  0.00000000e+00  1.32294142e-01]
 [-8.00000000e+00  0.00000000e+00  5.06497902e-01]
 [-7.00000000e+00  2.00000000e+00  1.65568919e+00]
 [-6.00000000e+00  1.00000000e+00  4.62125309e+00]
 [-5.00000000e+00  9.00000000e+00  1.10137298e+01]
 [-4.00000000e+00  2.60000000e+01  2.24137683e+01]
 [-3.00000000e+00  3.70000000e+01  3.89503370e+01]
 [-2.00000000e+00  5.10000000e+01  5.78004747e+01]
 [-1.00000000e+00  7.10000000e+01  7.32455414e+01]
 [ 0.00000000e+00  7.40000000e+01  7.92618251e+01]
 [ 1.00000000e+00  8.90000000e+01  7.32455414e+01]
 [ 2.00000000e+00  5.50000000e+01  5.78004747e+01]
 [ 3.00000000e+00  5.00000000e+01  3.89503370e+01]
 [ 4.00000000e+00  1.70000000e+01  2.24137683e+01]
 [ 5.00000000e+00  1.10000000e+01  1.10137298e+01]
 [ 6.00000000e+00  4.00000000e+00  4.62125309e+00]
 [ 7.00000000e+00  3.00000000e+00  1.65568919e+00]
 [ 8.00000000e+00  0.00000000e+00  5.06497902e-01]
 [ 9.00000000e+00  0.00000000e+00  1.32294142e-01]
 [ 1.00000000e+01  0.00000000e+00  2.95019349e-02]]
../_images/normdiscr_plot1.png
../_images/normdiscr_plot2.png
![在这里插入图片描述](https://img-blog.csdnimg.cn/20200605091613544.png?x-oss-process=image/watermark,type_ZmFuZ3poZW5naGVpdGk,shadow_10,text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3dlaXhpbl80MTEwMjY3Mg==,size_16,color_FFFFFF,t_70)
接下来，我们可以测试样本是否由我们的离散正态分布生成。这也可以验证随机数是否正确生成。

卡方检验要求每个仓中的观测值最少。我们将尾箱型合并为较大的箱型，以便它们包含足够的观测值。
```python
>>> f2 = np.hstack([f[:5].sum(), f[5:-5], f[-5:].sum()])
>>> p2 = np.hstack([probs[:5].sum(), probs[5:-5], probs[-5:].sum()])
>>> ch2, pval = stats.chisquare(f2, p2*n_sample)

>>> print('chisquare for normdiscrete: chi2 = %6.3f pvalue = %6.4f' % (ch2, pval))
chisquare for normdiscrete: chi2 = 12.466 pvalue = 0.4090

在这种情况下，p值很高，因此我们可以确信我们的随机样本实际上是由分布生成的。

分析一个样品

首先，我们创建一些随机变量。我们设置一个种子，以便在每次运行中得到相同的结果。例如，我们从Student t分布中抽取一个样本：

>>> np.random.seed(282629734)
>>> x = stats.t.rvs(10, size=1000)

在这里，我们将t分布的所需形状参数（在统计中对应于自由度）设置为10。size = 1000表示我们的样本包含1000个独立绘制的（伪）随机数。由于我们未指定关键字参数loc和scale，因此它们被系统分别设置为默认值0和1。

描述性统计

对于以上步骤生成的x这个numpy数组，我们可以直接访问所有数组方法，例如，

>>> print(x.min())   # 等效于 np.min(x)，求最小值
-3.78975572422
>>> print(x.max())   # 等效于 np.max(x)，求最大值
5.26327732981
>>> print(x.mean())  # 等效于 np.mean(x)，求均值
0.0140610663985
>>> print(x.var())   # 等效于 np.var(x)，求方差
1.28899386208

如果要检验样本性质与理论性质有多少差异，我们可以运行以下代码：

>>> m, v, s, k = stats.t.stats(10, moments='mvsk')# 标准分布的相关值。
>>> n, (smin, smax), sm, sv, ss, sk = stats.describe(x)# 对样本x进行描述统计所得到的值。

>>> sstr = '%-14s mean = %6.4f, variance = %6.4f, skew = %6.4f, kurtosis = %6.4f'
>>> print(sstr % ('distribution:', m, v, s ,k))
distribution:  mean = 0.0000, variance = 1.2500, skew = 0.0000, kurtosis = 1.0000
>>> print(sstr % ('sample:', sm, sv, ss, sk))
sample:        mean = 0.0141, variance = 1.2903, skew = 0.2165, kurtosis = 1.0556

注意：stats.describe使用无偏估计量作为方差，而np.var是有偏估计量。因此上面可以看到x.var()比后面的sv要小一些，盖因前者有偏，后者无偏。

对于我们的样本，样本统计量与理论值相差很小。

T检验和KS检验

我们可以使用t检验来检验样本均值是否在统计学上与理论预期有显着差异。

>>> print('t-statistic = %6.3f pvalue = %6.4f' %  stats.ttest_1samp(x, m))
t-statistic =  0.391 pvalue = 0.6955

p值为0.7，这意味着在alpha误差为10％的情况下，我们不能拒绝样本均值等于零（标准t分布的期望）的假设。

作为练习，我们也可以不使用提供的函数而直接计算ttest，该函数应该给我们相同的答案，所以它做到了：

>>> tt = (sm-m)/np.sqrt(sv/float(n))  # t-statistic for mean
>>> pval = stats.t.sf(np.abs(tt), n-1)*2  # two-sided pvalue = Prob(abs(t)>tt)
>>> print('t-statistic = %6.3f pvalue = %6.4f' % (tt, pval))
t-statistic =  0.391 pvalue = 0.6955

Kolmogorov-Smirnov检验可用于检验以下假设：样本来自标准t分布

>>> print('KS-statistic D = %6.3f pvalue = %6.4f' % stats.kstest(x, 't', (10,)))
KS-statistic D =  0.016 pvalue = 0.9606

同样，p值足够高，我们不能拒绝这样的假设：随机样本实际上是根据t分布进行分布的。在实际的应用中，我们不知道事实上的分布是什么。如果我们针对标准正态分布对样本执行Kolmogorov-Smirnov检验，则在此示例中，p值几乎为40％，我们也不能拒绝假设“样本是由正态分布生成的”。

>>> print('KS-statistic D = %6.3f pvalue = %6.4f' % stats.kstest(x, 'norm'))
KS-statistic D =  0.028 pvalue = 0.3949

但是，标准正态分布的方差为1，而我们的样本的方差为1.29。如果我们将样本标准化并再次针对正态分布进行测试，则p值将再次足够大，以至于我们无法拒绝样本来自正态分布的假设。

>>> d, pval = stats.kstest((x-x.mean())/x.std(), 'norm')
>>> print('KS-statistic D = %6.3f pvalue = %6.4f' % (d, pval))
KS-statistic D =  0.032 pvalue = 0.2402

注意：Kolmogorov-Smirnov检验假定我们针对给定参数的分布进行检验。但是又由于在最后一种的例子中，我们估计均值和方差时违反了这个假设，所以建立在测试分布上的p值是不正确的。<这段翻译存疑。>

分布尾巴

最后，我们可以检查分布的上尾。使用百分比点函数ppf（它是cdf函数的反函数）可以获取临界值，或者更直接地，可以使用生存函数的反函数isf，如下所示：

>>> crit01, crit05, crit10 = stats.t.ppf([1-0.01, 1-0.05, 1-0.10], 10)
>>> print('critical values from ppf at 1%%, 5%% and 10%% %8.4f %8.4f %8.4f' % (crit01, crit05, crit10))
critical values from ppf at 1%, 5% and 10%   2.7638   1.8125   1.3722
>>> print('critical values from isf at 1%%, 5%% and 10%% %8.4f %8.4f %8.4f' % tuple(stats.t.isf([0.01,0.05,0.10],10)))
critical values from isf at 1%, 5% and 10%   2.7638   1.8125   1.3722

>>> freq01 = np.sum(x>crit01) / float(n) * 100
>>> freq05 = np.sum(x>crit05) / float(n) * 100
>>> freq10 = np.sum(x>crit10) / float(n) * 100
>>> print('sample %%-frequency at 1%%, 5%% and 10%% tail %8.4f %8.4f %8.4f' % (freq01, freq05, freq10))
sample %-frequency at 1%, 5% and 10% tail   1.4000   5.8000  10.5000

在这三种情况下，我们的样本在顶部尾部的权重均大于其基础分布。我们可以简单地检查一个更大的样本，看看是否能找到更接近的匹配。在这种情况下，经验频率与理论概率非常接近，但是如果我们重复几次，则波动仍然很大。

>>> freq05l = np.sum(stats.t.rvs(10, size=10000) > crit05) / 10000.0 * 100
>>> print('larger sample %%-frequency at 5%% tail %8.4f' % freq05l)
larger sample %-frequency at 5% tail   4.8000

我们还可以将其与正态分布的尾部进行比较，该正态分布的尾部权重较小：

>>> print('tail prob. of normal at 1%%, 5%% and 10%% %8.4f %8.4f %8.4f' %
...       tuple(stats.norm.sf([crit01, crit05, crit10])*100))
tail prob. of normal at 1%, 5% and 10%   0.2857   3.4957   8.5003

卡方检验chisquare可用于检验在一定数量的统计数据中，观察到的频率是否与假设分布的概率显着不同。

>>> quantiles = [0.0, 0.01, 0.05, 0.1, 1-0.10, 1-0.05, 1-0.01, 1.0]
>>> crit = stats.t.ppf(quantiles, 10)
>>> crit
array([       -inf, -2.76376946, -1.81246112, -1.37218364,  1.37218364,
        1.81246112,  2.76376946,         inf])
>>> n_sample = x.size
>>> freqcount = np.histogram(x, bins=crit)[0]
>>> tprob = np.diff(quantiles)
>>> nprob = np.diff(stats.norm.cdf(crit))
>>> tch, tpval = stats.chisquare(freqcount, tprob*n_sample)
>>> nch, npval = stats.chisquare(freqcount, nprob*n_sample)
>>> print('chisquare for t:      chi2 = %6.2f pvalue = %6.4f' % (tch, tpval))
chisquare for t:      chi2 =  2.30 pvalue = 0.8901
>>> print('chisquare for normal: chi2 = %6.2f pvalue = %6.4f' % (nch, npval))
chisquare for normal: chi2 = 64.60 pvalue = 0.0000

我们看到标准正态分布明显被拒绝（p值几乎接近于0），而标准t分布却不能被拒绝。由于我们样本的方差不同于两个标准分布，因此我们可以再次考虑比例和位置的估计值来重做测试。

分布的拟合方法可用于估计分布的参数，并使用估计的分布的概率重复进行测试。

>>> tdof, tloc, tscale = stats.t.fit(x)
>>> nloc, nscale = stats.norm.fit(x)
>>> tprob = np.diff(stats.t.cdf(crit, tdof, loc=tloc, scale=tscale))
>>> nprob = np.diff(stats.norm.cdf(crit, loc=nloc, scale=nscale))
>>> tch, tpval = stats.chisquare(freqcount, tprob*n_sample)
>>> nch, npval = stats.chisquare(freqcount, nprob*n_sample)
>>> print('chisquare for t:      chi2 = %6.2f pvalue = %6.4f' % (tch, tpval))
chisquare for t:      chi2 =  1.58 pvalue = 0.9542
>>> print('chisquare for normal: chi2 = %6.2f pvalue = %6.4f' % (nch, npval))
chisquare for normal: chi2 = 11.08 pvalue = 0.0858

考虑到估计的参数，我们仍然可以拒绝“我们的样本来自正态分布（5％的水平）”的假设。但是同样地，当p值为0.95时，我们不能拒绝t分布。

针对正态分布的特殊测试

由于正态分布是统计中最常见的分布，因此有几个附加功能可用于测试样本是否来自正态分布。

首先，我们可以测试样本的偏斜和峰度是否与正态分布的偏斜和峰度显着不同：

>>> print('normal skewtest teststat = %6.3f pvalue = %6.4f' % stats.skewtest(x))
normal skewtest teststat =  2.785 pvalue = 0.0054
>>> print('normal kurtosistest teststat = %6.3f pvalue = %6.4f' % stats.kurtosistest(x))
normal kurtosistest teststat =  4.757 pvalue = 0.0000

这两个测试已被合并到正态性检验函数normaltest中

>>> print('normaltest teststat = %6.3f pvalue = %6.4f' % stats.normaltest(x))
normaltest teststat = 30.379 pvalue = 0.0000

在所有三个测试中，p值都非常低，我们可以拒绝以下假设：“我们的样本具有正态分布的偏斜和峰度。”

由于我们样本的偏斜和峰度基于中心矩，因此，如果我们测试标准化样本，我们将获得精确相同的结果：

>>> print('normaltest teststat = %6.3f pvalue = %6.4f' %
...       stats.normaltest((x-x.mean())/x.std()))
normaltest teststat = 30.379 pvalue = 0.0000

可以看见p值（pvalue）是0，说明其正态性被强烈地拒绝。由此，我们可以检查正态检验是否可在其他情况下给出了合理的结果：

>>> print('normaltest teststat = %6.3f pvalue = %6.4f' %
...       stats.normaltest(stats.t.rvs(10, size=100)))
normaltest teststat =  4.698 pvalue = 0.0955
>>> print('normaltest teststat = %6.3f pvalue = %6.4f' %
...              stats.normaltest(stats.norm.rvs(size=1000)))
normaltest teststat =  0.613 pvalue = 0.7361

当测试小样本t分布观测值和大样本正态分布观测值的正态性时，在任何情况下，我们都不能拒绝样本来自正态分布的原假设。在第一种情况下，这是因为测试的能力不足以区分小样本中的t分布随机变量和正态分布的随机变量。

比较两个样本

下面，我们给了两个样本，它们可以来自相同或不同的分布。我们要测试这些样本是否具有相同的统计特性。

比较手段

用相同的方法测试样品：

>>> rvs1 = stats.norm.rvs(loc=5, scale=10, size=500)
>>> rvs2 = stats.norm.rvs(loc=5, scale=10, size=500)
>>> stats.ttest_ind(rvs1, rvs2)
Ttest_indResult(statistic=-0.5489036175088705, pvalue=0.5831943748663959)

用不同的方式测试样本：

>>> rvs3 = stats.norm.rvs(loc=8, scale=10, size=500)
>>> stats.ttest_ind(rvs1, rvs3)
Ttest_indResult(statistic=-4.533414290175026, pvalue=6.507128186389019e-06)

两个样本ks_2samp的Kolmogorov-Smirnov测试

对于该示例，两个样本均来自相同的分布，我们不能拒绝原假设，因为p值很高

>>> stats.ks_2samp(rvs1, rvs2)
Ks_2sampResult(statistic=0.026, pvalue=0.9959527565364388)

在第下一个示例中，由于p值低于1％，因此在不同的位置（即均值）下，我们可以拒绝原假设。

>>> stats.ks_2samp(rvs1, rvs3)
Ks_2sampResult(statistic=0.114, pvalue=0.00299005061044668)

核密度估计<翻译到这里>

统计中的一项常见任务是从一组数据样本中估计随机变量的概率密度函数（PDF）。此任务称为密度估计。

对于这种工作，直方图是最著名的工具。直方图是用于可视化的有用工具（主要是因为每个人都了解它），但是不能非常有效地使用可用数据。而内核密度估计（KDE）则是用于同一任务的更有效的工具。stats内置的gaussian_kde估计可以用来估计单变量的PDF以及多变量数据。如果数据是单峰的，则效果最好。

单因素估计

我们以最少的数据量开始，以了解其gaussian_kde 的工作方式以及带宽选择的不同选项的作用。从PDF采样的数据在图的底部显示为蓝色虚线（这称为地毯图）：

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

>>> x1 = np.array([-7, -5, 1, 4, 5], dtype=np.float)
>>> kde1 = stats.gaussian_kde(x1)
>>> kde2 = stats.gaussian_kde(x1, bw_method='silverman')

>>> fig = plt.figure()
>>> ax = fig.add_subplot(111)
```py
>>> ax.plot(x1, np.zeros(x1.shape), 'b+', ms=20)  # rug plot
>>> x_eval = np.linspace(-10, 10, num=200)
>>> ax.plot(x_eval, kde1(x_eval), 'k-', label="Scott's Rule")
>>> ax.plot(x_eval, kde2(x_eval), 'r-', label="Silverman's Rule")

>>> plt.show()

我们看到，斯科特规则(Scott’s Rule)和西尔弗曼规则(Silverman’s Rule)之间几乎没有什么区别，并且数据量有限的带宽选择可能有点太宽。我们可以定义自己的带宽函数以获得较少被平滑的结果。

>>> def my_kde_bandwidth(obj, fac=1./5):
...     """We use Scott's Rule, multiplied by a constant factor."""
...     return np.power(obj.n, -1./(obj.d+4)) * fac

>>> fig = plt.figure()
>>> ax = fig.add_subplot(111)

>>> ax.plot(x1, np.zeros(x1.shape), 'b+', ms=20)  # rug plot
>>> kde3 = stats.gaussian_kde(x1, bw_method=my_kde_bandwidth)
>>> ax.plot(x_eval, kde3(x_eval), 'g-', label="With smaller BW")

>>> plt.show()

我们看到，如果将带宽设置得很窄，则获得的概率密度函数（PDF）估计值就是每个数据点周围的高斯分布之和。<对吗?>

现在，我们以一个更现实的示例为例，看看两个可用带宽选择规则之间的区别。众所周知，这些规则对于（接近）正态分布非常有效，但即使对于非正态分布非常强烈的单峰分布，它们也相当有效。作为非正态分布，我们采用自由度为5的学生T分布。

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


np.random.seed(12456)
x1 = np.random.normal(size=200)  # random data, normal distribution
xs = np.linspace(x1.min()-1, x1.max()+1, 200)

kde1 = stats.gaussian_kde(x1)
kde2 = stats.gaussian_kde(x1, bw_method='silverman')

fig = plt.figure(figsize=(8, 6))

ax1 = fig.add_subplot(211)
ax1.plot(x1, np.zeros(x1.shape), 'b+', ms=12)  # rug plot
ax1.plot(xs, kde1(xs), 'k-', label="Scott's Rule")
ax1.plot(xs, kde2(xs), 'b-', label="Silverman's Rule")
ax1.plot(xs, stats.norm.pdf(xs), 'r--', label="True PDF")

ax1.set_xlabel('x')
ax1.set_ylabel('Density')
ax1.set_title("Normal (top) and Student's T$_{df=5}$ (bottom) distributions")
ax1.legend(loc=1)

x2 = stats.t.rvs(5, size=200)  # random data, T distribution
xs = np.linspace(x2.min() - 1, x2.max() + 1, 200)

kde3 = stats.gaussian_kde(x2)
kde4 = stats.gaussian_kde(x2, bw_method='silverman')

ax2 = fig.add_subplot(212)
ax2.plot(x2, np.zeros(x2.shape), 'b+', ms=12)  # rug plot
ax2.plot(xs, kde3(xs), 'k-', label="Scott's Rule")
ax2.plot(xs, kde4(xs), 'b-', label="Silverman's Rule")
ax2.plot(xs, stats.t.pdf(xs, 5), 'r--', label="True PDF")

ax2.set_xlabel('x')
ax2.set_ylabel('Density')

plt.show()

现在我们来看看具有一个较宽和一个较窄的高斯特征的双峰分布。由于准确解析每个功能所需的带宽不同，我们预计这将是一个更难估计的密度。

>>> from functools import partial

>>> loc1, scale1, size1 = (-2, 1, 175)
>>> loc2, scale2, size2 = (2, 0.2, 50)
>>> x2 = np.concatenate([np.random.normal(loc=loc1, scale=scale1, size=size1),
...                      np.random.normal(loc=loc2, scale=scale2, size=size2)])

>>> x_eval = np.linspace(x2.min() - 1, x2.max() + 1, 500)

>>> kde = stats.gaussian_kde(x2)
>>> kde2 = stats.gaussian_kde(x2, bw_method='silverman')
>>> kde3 = stats.gaussian_kde(x2, bw_method=partial(my_kde_bandwidth, fac=0.2))
>>> kde4 = stats.gaussian_kde(x2, bw_method=partial(my_kde_bandwidth, fac=0.5))

>>> pdf = stats.norm.pdf
>>> bimodal_pdf = pdf(x_eval, loc=loc1, scale=scale1) * float(size1) / x2.size + \
...               pdf(x_eval, loc=loc2, scale=scale2) * float(size2) / x2.size

>>> fig = plt.figure(figsize=(8, 6))
>>> ax = fig.add_subplot(111)

>>> ax.plot(x2, np.zeros(x2.shape), 'b+', ms=12)
>>> ax.plot(x_eval, kde(x_eval), 'k-', label="Scott's Rule")
>>> ax.plot(x_eval, kde2(x_eval), 'b-', label="Silverman's Rule")
>>> ax.plot(x_eval, kde3(x_eval), 'g-', label="Scott * 0.2")
>>> ax.plot(x_eval, kde4(x_eval), 'c-', label="Scott * 0.5")
>>> ax.plot(x_eval, bimodal_pdf, 'r--', label="Actual PDF")

>>> ax.set_xlim([x_eval.min(), x_eval.max()])
>>> ax.legend(loc=2)
>>> ax.set_xlabel('x')
>>> ax.set_ylabel('Density')
>>> plt.show()

不出所料，由于双峰分布的两个特征的特征尺寸不同，KDE并不像我们想要的那样接近真实的PDF。通过将默认带宽减半（Scott*0.5），我们可以做得更好，同时使用比默认带宽小5倍的带宽还不够平滑。但是，在这种情况下，我们真正需要的是非均匀（自适应）带宽。Scott * 0.5

多因素估计

有了gaussian_kde，我们便可以执行多变量估计，就和单变量估计的方式相似。我们证明了双变量情况。首先，我们创建一个模型，在这个模型中两个变量相关。然后使用这个模型，生成一些随机数据。

>>> def measure(n):
...     """Measurement model, return two coupled measurements."""
...     m1 = np.random.normal(size=n)
...     m2 = np.random.normal(scale=0.5, size=n)
...     return m1+m2, m1-m2

>>> m1, m2 = measure(2000)
>>> xmin = m1.min()
>>> xmax = m1.max()
>>> ymin = m2.min()
>>> ymax = m2.max()

然后我们将KDE应用于数据：

>>> X, Y = np.mgrid[xmin:xmax:100j, ymin:ymax:100j]
>>> positions = np.vstack([X.ravel(), Y.ravel()])
>>> values = np.vstack([m1, m2])
>>> kernel = stats.gaussian_kde(values)
>>> Z = np.reshape(kernel.evaluate(positions).T, X.shape)

最后，我们将估计的双变量分布绘制为颜色图，并在顶部绘制各个数据点。

>>> fig = plt.figure(figsize=(8, 6))
>>> ax = fig.add_subplot(111)

>>> ax.imshow(np.rot90(Z), cmap=plt.cm.gist_earth_r,
...           extent=[xmin, xmax, ymin, ymax])
>>> ax.plot(m1, m2, 'k.', markersize=2)

>>> ax.set_xlim([xmin, xmax])
>>> ax.set_ylim([ymin, ymax])

>>> plt.show()

多尺度图相关（MGC）<翻译到这里，以下未翻译>

使用multiscale_graphcorr，我们可以测试高维和非线性数据的独立性。在开始之前，让我们导入一些有用的软件包：

import numpy as np
import matplotlib.pyplot as plt; plt.style.use(‘classic’)
from scipy.stats import multiscale_graphcorr
让我们使用一个自定义绘图功能来绘制数据关系：

def mgc_plot(x, y, sim_name, mgc_dict=None, only_viz=False,
… only_mgc=False):
… “”“Plot sim and MGC-plot”""
… if not only_mgc:
… # simulation
… plt.figure(figsize=(8, 8))
… ax = plt.gca()
… ax.set_title(sim_name + " Simulation", fontsize=20)
… ax.scatter(x, y)
… ax.set_xlabel(‘X’, fontsize=15)
… ax.set_ylabel(‘Y’, fontsize=15)
… ax.axis(‘equal’)
… ax.tick_params(axis=“x”, labelsize=15)
… ax.tick_params(axis=“y”, labelsize=15)
… plt.show()
… if not only_viz:
… # local correlation map
… plt.figure(figsize=(8,8))
… ax = plt.gca()
… mgc_map = mgc_dict[“mgc_map”]
… # draw heatmap
… ax.set_title(“Local Correlation Map”, fontsize=20)
… im = ax.imshow(mgc_map, cmap=‘YlGnBu’)
… # colorbar
… cbar = ax.figure.colorbar(im, ax=ax)
… cbar.ax.set_ylabel("", rotation=-90, va=“bottom”)
… ax.invert_yaxis()
… # Turn spines off and create white grid.
… for edge, spine in ax.spines.items():
… spine.set_visible(False)
… # optimal scale
… opt_scale = mgc_dict[“opt_scale”]
… ax.scatter(opt_scale[0], opt_scale[1],
… marker=‘X’, s=200, color=‘red’)
… # other formatting
… ax.tick_params(bottom=“off”, left=“off”)
… ax.set_xlabel(’#Neighbors for X’, fontsize=15)
… ax.set_ylabel(’#Neighbors for Y’, fontsize=15)
… ax.tick_params(axis=“x”, labelsize=15)
… ax.tick_params(axis=“y”, labelsize=15)
… ax.set_xlim(0, 100)
… ax.set_ylim(0, 100)
… plt.show()
让我们先来看一些线性数据：

np.random.seed(12345678)
x = np.linspace(-1, 1, num=100)
y = x + 0.3 * np.random.random(x.size)
仿真关系可以画在下面：

mgc_plot(x, y, “Linear”, only_viz=True)
…/_images/mgc_plot1_01_00.png
现在，我们可以在下面看到可视化的测试统计量，p值和MGC图。最佳比例在地图上显示为红色“ x”：

stat, pvalue, mgc_dict = multiscale_graphcorr(x, y)
print("MGC test statistic: ", round(stat, 1))
MGC test statistic: 1.0

print("P-value: ", round(pvalue, 1))
P-value: 0.0

mgc_plot(x, y, “Linear”, mgc_dict, only_mgc=True)
…/_images/mgc_plot2.png
从这里可以清楚地看到，MGC能够确定输入数据矩阵之间的关系，因为p值非常低，而MGC测试统计信息相对较高。MGC-map表示强烈的线性关系。直觉上，这是因为拥有更多的邻居将有助于确定之间的线性关系。 和 。在这种情况下，最佳比例等于全局比例，在地图上用红色斑点标记。

对于非线性数据集也可以这样做。以下 和 数组是从非线性仿真得出的：

np.random.seed(12345678)
unif = np.array(np.random.uniform(0, 5, size=100))
x = unif * np.cos(np.pi * unif)
y = unif * np.sin(np.pi * unif) + 0.4 * np.random.random(x.size)
仿真关系可以画在下面：

mgc_plot(x, y, “Spiral”, only_viz=True)
…/_images/mgc_plot3_01_00.png
现在，我们可以在下面看到可视化的测试统计量，p值和MGC图。最佳比例在地图上显示为红色“ x”：

stat, pvalue, mgc_dict = multiscale_graphcorr(x, y)
print("MGC test statistic: ", round(stat, 1))
MGC test statistic: 0.2

print("P-value: ", round(pvalue, 1))
P-value: 0.0

mgc_plot(x, y, “Spiral”, mgc_dict, only_mgc=True)
…/_images/mgc_plot4.png
从这里可以清楚地看到，由于p值非常低且MGC测试统计值相对较高，因此MGC可以再次确定关系。MGC映射表示强烈的非线性关系。在这种情况下，最佳比例等于局部比例，在地图上用红色斑点标记。