统计 (scipy.stats)#
简介#
在本教程中，我们将讨论 scipy.stats 的许多（但肯定不是全部）功能。这里的目的是为用户提供对该包的工作知识。有关更多详细信息，请参阅 参考手册。
注意：此文档正在进行中。
随机变量#
已经实现了两个通用的分布类，用于封装 连续随机变量 和 离散随机变量。使用这些类已经实现了 80 多个连续随机变量 (RV) 和 10 个离散随机变量。除此之外，最终用户可以轻松添加新的例程和分布。（如果您创建了一个，请贡献它。）
所有统计函数都位于子包 scipy.stats 中，可以使用 info(stats) 获取这些函数的相当完整的列表。还可以从 stats 子包的文档字符串中获取可用随机变量的列表。
在下面的讨论中，我们主要关注连续 RV。几乎所有内容也适用于离散变量，但我们在这里指出了一些差异：离散分布的特定点。
在下面的代码示例中，我们假设 scipy.stats 包被导入为
>>> from scipy import stats
在某些情况下，我们假设单个对象被导入为
>>> from scipy.stats import norm
获取帮助#
首先，所有分布都附带帮助函数。要获取一些基本信息，我们可以打印相关的文档字符串：print(stats.norm.__doc__)。
要查找支持，即分布的上限和下限，请调用
>>> print('bounds of distribution lower: %s, upper: %s' % norm.support())
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']
最后，我们可以通过内省获取可用分布的列表
>>> 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: 108
>>> print('number of discrete distributions: %d' % len(dist_discrete))
number of discrete distributions: 20
常用方法#
连续 RV 的主要公共方法是
rvs：随机变量
pdf：概率密度函数
cdf：累积分布函数
sf：生存函数（1-CDF）
ppf：百分位点函数（CDF 的逆函数）
isf：逆生存函数（SF 的逆函数）
stats：返回均值、方差、（费舍尔的）偏度或（费舍尔的）峰度
moment：分布的非中心矩
让我们以一个正态 RV 为例。
>>> 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
要生成一系列随机变量，请使用 size 关键字参数
>>> norm.rvs(size=3)
array([-0.35687759, 1.34347647, -0.11710531]) # random
不要认为 norm.rvs(5) 会生成 5 个变量
>>> norm.rvs(5)
5.471435163732493 # random
这里，5 没有关键字，被解释为第一个可能的关键字参数 loc，它是所有连续分布接受的两个关键字参数中的第一个。这将我们引入了下一小节的主题。
随机数生成#
随机数的绘制依赖于来自 numpy.random 包的生成器。在上面的示例中，随机数的特定流在运行之间不可重现。为了实现可重现性，您可以显式地 *播种* 随机数生成器。在 NumPy 中，生成器是 numpy.random.Generator 的实例。以下是创建生成器的规范方法
>>> from numpy.random import default_rng
>>> rng = default_rng()
固定种子可以这样完成
>>> # do NOT copy this value
>>> rng = default_rng(301439351238479871608357552876690613766)
警告
不要使用此数字或常见的数字，例如 0。使用一小部分种子来实例化更大的状态空间意味着有些初始状态无法到达。如果每个人都使用这些值，就会产生一些偏差。获取种子的好方法是使用 numpy.random.SeedSequence
>>> from numpy.random import SeedSequence
>>> print(SeedSequence().entropy)
301439351238479871608357552876690613766 # random
分布中的 random_state 参数接受 numpy.random.Generator 类的实例，或一个整数，然后用于播种内部的 Generator 对象
>>> norm.rvs(size=5, random_state=rng)
array([ 0.47143516, -1.19097569, 1.43270697, -0.3126519 , -0.72058873]) # random
有关更多信息，请参阅 NumPy 的文档。
要了解有关 SciPy 中实现的随机数采样器的更多信息，请参阅 非均匀随机数采样教程 和 准蒙特卡罗教程
移位和缩放#
所有连续分布都接受 loc 和 scale 作为关键字参数来调整分布的位置和尺度，例如，对于标准正态分布，位置是均值，尺度是标准差。
>>> norm.stats(loc=3, scale=4, moments="mv")
(array(3.0), array(16.0))
在许多情况下，随机变量 X 的标准化分布是通过变换 (X - loc) / scale 获得的。默认值为 loc = 0 和 scale = 1。
巧妙地使用 loc 和 scale 可以帮助以多种方式修改标准分布。为了进一步说明缩放，具有均值为 \(1/\lambda\) 的指数分布 RV 的 cdf 由下式给出
通过应用上面的缩放规则，可以看出，通过取 scale = 1./lambda，我们可以得到正确的尺度。
>>> from scipy.stats import expon
>>> expon.mean(scale=3.)
3.0
注意
接受形状参数的分布可能需要比简单地应用 loc 和/或 scale 更复杂的操作才能达到所需的形状。例如，给定长度为 \(R\) 的常数向量，在每个分量中受到独立的 N(0, \(\sigma^2\)) 偏差扰动后，2-D 向量长度的分布为 rice(\(R/\sigma\), scale= \(\sigma\))。第一个参数是一个形状参数，需要与 \(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 # random
因此，为了解释上一节示例的输出：norm.rvs(5) 生成一个均值为 loc=5 的单一正态分布随机变量，因为默认值为 size=1。
我们建议您明确设置 loc 和 scale 参数，通过将值作为关键字而不是参数传递。当调用给定 RV 的多个方法时，可以通过使用 冻结分布 的技术来最大程度地减少重复，如下所述。
形状参数#
虽然可以使用 loc 和 scale 参数对一般连续随机变量进行平移和缩放，但某些分布需要额外的形状参数。例如，密度为
的伽马分布需要形状参数 \(a\)。观察到通过将 scale 关键字设置为 \(1/\lambda\) 可以获得 \(\lambda\) 的设置。
让我们检查伽马分布的形状参数的数量和名称。（我们从上面知道这应该是 1。）
>>> from scipy.stats import gamma
>>> gamma.numargs
1
>>> gamma.shapes
'a'
现在，我们将形状变量的值设置为 1 以获得指数分布，以便我们轻松比较是否获得了预期的结果。
>>> gamma(1, scale=2.).stats(moments="mv")
(array(2.0), array(4.0))
请注意，我们也可以将形状参数指定为关键字
>>> gamma(a=1, scale=2.).stats(moments="mv")
(array(2.0), array(4.0))
冻结分布#
反复传递 loc 和 scale 关键字可能会变得很麻烦。使用 冻结 RV 的概念来解决此类问题。
>>> rv = gamma(1, scale=2.)
通过使用 rv，我们不再需要包含比例或形状参数。因此，分布可以使用两种方式之一，要么将所有分布参数传递给每个方法调用（就像我们之前所做的那样），要么冻结分布实例的参数。让我们检查一下
>>> rv.mean(), rv.std()
(2.0, 2.0)
这确实是我们应该得到的。
广播#
基本方法 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 个自由度的临界值（d.o.f.）。因此，广播规则给出了与调用 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% 尾部：
>>> stats.t.isf([0.1, 0.05, 0.01], [10, 11, 12])
array([ 1.37218364, 1.79588482, 2.68099799])
离散分布的具体要点#
离散分布与连续分布具有基本相同的基本方法。但是，pdf 被概率质量函数 pmf 替换，没有可用的估计方法，例如拟合，并且 scale 不是有效的关键字参数。位置参数，关键字 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 阶梯函数拐点的值，我们会得到下一个更高的整数
>>> 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：分布参数的最大似然估计，包括位置
和尺度
fit_loc_scale：当给出形状参数时，位置和尺度的估计
nnlf：负对数似然函数
expect：计算函数相对于 pdf 或 pmf 的期望
性能问题和注意事项#
各个方法的性能，在速度方面，因分布和方法而异。方法的结果通过两种方式之一获得：要么通过显式计算，要么通过独立于特定分布的通用算法。
一方面，显式计算要求直接为给定分布指定方法，无论是通过解析公式，还是通过 scipy.special 或 numpy.random 中的特殊函数来实现 rvs。这些通常是相对快速的计算。
另一方面，如果分布没有指定任何显式计算，则使用通用方法。要定义一个分布，只需要 pdf 或 cdf 中的一个；所有其他方法都可以使用数值积分和求根来推导。但是，这些间接方法可能非常慢。例如，rgh = stats.gausshyper.rvs(0.5, 2, 2, 2, size=100) 以非常间接的方式创建随机变量，在我的计算机上，100 个随机变量需要大约 19 秒，而标准正态分布或 t 分布的一百万个随机变量只需要一秒多一点。
剩余问题#
在 scipy.stats 中的分布最近得到了修正和改进，并获得了相当多的测试套件；然而，仍然存在一些问题。
这些分布已经在一定范围的参数上进行了测试；然而，在一些角落范围内，可能仍然存在一些错误的结果。
在 fit 中的最大似然估计不适用于所有分布的默认起始参数，用户需要提供良好的起始参数。此外，对于某些分布，使用最大似然估计可能本质上不是最佳选择。
构建特定分布#
以下示例展示了如何构建自己的分布。进一步的示例展示了分布的使用和一些统计检验。
创建连续分布，即子类化 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 在确定性分布的类定义中没有被指定。
子类化 rv_discrete#
在以下内容中，我们使用 stats.rv_discrete 生成一个离散分布，该分布具有围绕整数为中心的区间的截断正态分布的概率。
一般信息
从 rv_discrete 的文档字符串中，help(stats.rv_discrete)，
“您可以通过将包含两个序列元组 (xk, pk) 传递给 rv_discrete 初始化方法（通过 values= 关键字）来构建一个任意的离散随机变量，其中 P{X=xk} = pk，该元组仅描述了 X (xk) 的那些以非零概率 (pk) 出现的取值。”
除此之外，这种方法还需要满足一些其他要求才能正常工作。
关键字 name 是必需的。
分布的支撑点 xk 必须是整数。
需要指定有效数字（小数）的数量。
事实上，如果最后两个要求没有满足，可能会引发异常或导致结果数字不正确。
示例
让我们来做一下。首先
>>> npoints = 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 # used later in the analysis
>>> 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
>>> 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] # random
[-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]]
接下来，我们可以测试我们的样本是否是由我们的 norm-discrete 分布生成的。这也验证了随机数是否生成正确。
卡方检验要求每个箱体中至少有最小数量的观测值。我们将尾部箱体合并成更大的箱体，以便它们包含足够的观测值。
>>> 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 # random
在这种情况下，p 值很高，因此我们可以非常自信地认为我们的随机样本实际上是由该分布生成的。
分析一个样本#
首先，我们创建一些随机变量。我们设置一个种子，以便在每次运行中获得相同的結果。例如，我们从学生 t 分布中抽取一个样本
>>> x = stats.t.rvs(10, size=1000)
在这里，我们将 t 分布的所需形状参数（在统计学中对应于自由度）设置为 10。使用 size=1000 表示我们的样本包含 1000 个独立抽取的（伪）随机数。由于我们没有指定关键字参数 loc 和 scale，因此它们被设置为默认值零和一。
描述性统计#
x 是一个 numpy 数组，我们可以直接访问所有数组方法，例如：
>>> print(x.min()) # equivalent to np.min(x)
-3.78975572422 # random
>>> print(x.max()) # equivalent to np.max(x)
5.26327732981 # random
>>> print(x.mean()) # equivalent to np.mean(x)
0.0140610663985 # random
>>> print(x.var()) # equivalent to np.var(x))
1.28899386208 # random
样本属性与其理论对应物相比如何？
>>> m, v, s, k = stats.t.stats(10, moments='mvsk')
>>> n, (smin, smax), sm, sv, ss, sk = stats.describe(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 # random
>>> print(sstr % ('sample:', sm, sv, ss, sk))
sample: mean = 0.0141, variance = 1.2903, skew = 0.2165, kurtosis = 1.0556 # random
注意：stats.describe 使用方差的无偏估计量，而 np.var 是有偏估计量。
对于我们的样本，样本统计量与其理论对应物之间存在很小的差异。
t 检验和 KS 检验#
我们可以使用 t 检验来检验样本均值是否与理论预期值在统计学上存在显著差异。
>>> print('t-statistic = %6.3f pvalue = %6.4f' % stats.ttest_1samp(x, m))
t-statistic = 0.391 pvalue = 0.6955 # random
p 值为 0.7，这意味着在例如 10% 的 alpha 错误率下，我们不能拒绝样本均值等于零（标准 t 分布的期望值）的假设。
作为练习，我们可以直接计算 t 检验，而无需使用提供的函数，这应该会给我们相同的答案，事实也确实如此。
>>> 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 # random
Kolmogorov-Smirnov 检验可用于检验样本是否来自标准 t 分布的假设。
>>> print('KS-statistic D = %6.3f pvalue = %6.4f' % stats.kstest(x, 't', (10,)))
KS-statistic D = 0.016 pvalue = 0.9571 # random
同样，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.3918 # random
但是，标准正态分布的方差为 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.2397 # random
注意：Kolmogorov-Smirnov 检验假设我们针对具有给定参数的分布进行检验，因为在最后一种情况下，我们估计了均值和方差，所以这个假设被违反了，并且基于 p 值的检验统计量的分布是不正确的。
分布的尾部#
最后，我们可以检查分布的上尾。我们可以使用百分位点函数 ppf（它是 cdf 函数的逆函数）来获得临界值，或者更直接地说，我们可以使用生存函数的逆函数。
>>> 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 # random
在这三种情况下，我们的样本在顶端尾部比基础分布具有更大的权重。我们可以简要检查一个更大的样本，看看我们是否能得到更接近的匹配。在这种情况下，经验频率非常接近理论概率，但如果我们重复此操作几次，波动仍然相当大。
>>> 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 # random
我们也可以将其与正态分布的尾部进行比较，正态分布的尾部权重较小。
>>> 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
卡方检验可用于检验对于有限数量的箱子，观察到的频率是否与假设分布的概率显著不同。
>>> 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 # random
>>> print('chisquare for normal: chi2 = %6.2f pvalue = %6.4f' % (nch, npval))
chisquare for normal: chi2 = 64.60 pvalue = 0.0000 # random
我们看到标准正态分布被明确拒绝，而标准 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 # random
>>> print('chisquare for normal: chi2 = %6.2f pvalue = %6.4f' % (nch, npval))
chisquare for normal: chi2 = 11.08 pvalue = 0.0858 # random
考虑到估计参数，我们仍然可以拒绝我们的样本来自正态分布的假设（在 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 # random
>>> print('normal kurtosistest teststat = %6.3f pvalue = %6.4f' % stats.kurtosistest(x))
normal kurtosistest teststat = 4.757 pvalue = 0.0000 # random
这两个检验在正态性检验中结合在一起。
>>> print('normaltest teststat = %6.3f pvalue = %6.4f' % stats.normaltest(x))
normaltest teststat = 30.379 pvalue = 0.0000 # random
在这三种检验中，p 值都非常低，我们可以拒绝我们的样本具有正态分布的偏度和峰度的假设。
由于我们样本的偏度和峰度是基于中心矩的，因此如果我们检验标准化样本，我们将得到完全相同的结果。
>>> print('normaltest teststat = %6.3f pvalue = %6.4f' %
... stats.normaltest((x-x.mean())/x.std()))
normaltest teststat = 30.379 pvalue = 0.0000 # random
由于正态性被强烈拒绝，我们可以检查 normaltest 是否对其他情况给出合理的结果。
>>> print('normaltest teststat = %6.3f pvalue = %6.4f' %
... stats.normaltest(stats.t.rvs(10, size=100)))
normaltest teststat = 4.698 pvalue = 0.0955 # random
>>> print('normaltest teststat = %6.3f pvalue = %6.4f' %
... stats.normaltest(stats.norm.rvs(size=1000)))
normaltest teststat = 0.613 pvalue = 0.7361 # random
当检验 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) # random
具有不同均值的样本检验
>>> rvs3 = stats.norm.rvs(loc=8, scale=10, size=500)
>>> stats.ttest_ind(rvs1, rvs3)
Ttest_indResult(statistic=-4.533414290175026, pvalue=6.507128186389019e-06) # random
两个样本的 Kolmogorov-Smirnov 检验 ks_2samp#
对于两个样本都来自相同分布的示例，我们无法拒绝零假设，因为 p 值很高。
>>> stats.ks_2samp(rvs1, rvs2)
KstestResult(statistic=0.026, pvalue=0.9959527565364388) # random
在第二个示例中，具有不同的位置，即均值，我们可以拒绝零假设，因为 p 值低于 1% 。
>>> stats.ks_2samp(rvs1, rvs3)
KstestResult(statistic=0.114, pvalue=0.00299005061044668) # random
核密度估计#
统计学中的一项常见任务是从一组数据样本中估计随机变量的概率密度函数 (PDF)。这项任务称为密度估计。最著名的工具是直方图。直方图是可视化的有用工具（主要是因为每个人都理解它），但它并没有有效地利用可用数据。核密度估计 (KDE) 是用于同一任务的更有效的工具。 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.float64)
>>> kde1 = stats.gaussian_kde(x1)
>>> kde2 = stats.gaussian_kde(x1, bw_method='silverman')
>>> fig = plt.figure()
>>> ax = fig.add_subplot(111)
>>> 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 规则和 Silverman 规则之间几乎没有区别，并且使用有限数量的数据进行的带宽选择可能有点太宽。我们可以定义自己的带宽函数来获得更平滑的结果。
>>> 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
rng = np.random.default_rng()
x1 = rng.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_state=rng) # 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 倍的带宽则不够平滑。然而，在这种情况下，我们真正需要的是非均匀（自适应）带宽。
多元估计#
使用 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()
让我们先看看一些线性数据
>>> rng = np.random.default_rng()
>>> x = np.linspace(-1, 1, num=100)
>>> y = x + 0.3 * rng.random(x.size)
模拟关系可以在下面绘制
>>> mgc_plot(x, y, "Linear", only_viz=True)
现在，我们可以看到下面可视化的检验统计量、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)
从这里可以清楚地看到，MGC 能够确定输入数据矩阵之间的关系，因为 p 值非常低，而 MGC 检验统计量相对较高。MGC 图表明了**强烈的线性关系**。直观地说，这是因为拥有更多邻居将有助于识别 \(x\) 和 \(y\) 之间的线性关系。在这种情况下，最佳尺度**等同于全局尺度**，在图上用红色点标记。
非线性数据集也可以进行同样的操作。以下 \(x\) 和 \(y\) 数组来自非线性模拟
>>> unif = np.array(rng.uniform(0, 5, size=100))
>>> x = unif * np.cos(np.pi * unif)
>>> y = unif * np.sin(np.pi * unif) + 0.4 * rng.random(x.size)
模拟关系可以在下面绘制
>>> mgc_plot(x, y, "Spiral", only_viz=True)
现在，我们可以看到下面可视化的检验统计量、p 值和 MGC 图。最佳尺度在图上显示为红色“x”
>>> stat, pvalue, mgc_dict = multiscale_graphcorr(x, y)
>>> print("MGC test statistic: ", round(stat, 1))
MGC test statistic: 0.2 # random
>>> print("P-value: ", round(pvalue, 1))
P-value: 0.0
>>> mgc_plot(x, y, "Spiral", mgc_dict, only_mgc=True)
从这里可以清楚地看到，MGC 能够再次确定一种关系，因为 p 值非常低，而 MGC 检验统计量相对较高。MGC 地图表明存在 **强烈的非线性关系**。在这种情况下，最佳尺度 **相当于局部尺度**，在地图上用红点标记。
准蒙特卡罗#
在谈论准蒙特卡罗 (QMC) 之前，先简单介绍一下蒙特卡罗 (MC)。MC 方法或 MC 实验是一类广泛的计算算法，它们依赖于重复随机抽样来获得数值结果。其基本概念是利用随机性来解决原则上可能是确定性的问题。它们经常用于物理和数学问题，在难以或不可能使用其他方法时最为有用。MC 方法主要用于三类问题：优化、数值积分和从概率分布中生成抽样。
生成具有特定属性的随机数比听起来更复杂。简单的 MC 方法旨在对点进行采样，使其独立且同分布 (IID)。但生成多组随机点可能会产生截然不同的结果。
在上图的两种情况下，点都是随机生成的，没有任何关于先前绘制点的知识。很明显，空间的某些区域没有被探索 - 这可能会在模拟中造成问题，因为特定的一组点可能会触发完全不同的行为。
MC 的一个巨大优势是它具有已知的收敛特性。让我们看一下 5 维空间中平方和的平均值
其中 \(x_j \sim \mathcal{U}(0,1)\)。它有一个已知的平均值，\(\mu = 5/3+5(5-1)/4\)。使用 MC 抽样，我们可以用数值方法计算该平均值，并且近似误差遵循 \(O(n^{-1/2})\) 的理论速率。
虽然收敛性得到保证，但从业人员往往希望有一个更确定的探索过程。对于正常的 MC，可以使用种子来获得可重复的过程。但固定种子会破坏收敛特性：给定种子可能适用于给定类的问题，而对于另一类问题则会失效。
通常的做法是使用跨越所有参数维度的规则网格来以确定性方式遍历空间，也称为饱和设计。让我们考虑单位超立方体，所有边界范围从 0 到 1。现在，如果点之间的距离为 0.1，则填充单位区间的点数将为 10。在二维超立方体中，相同的间距将需要 100 个点，在三维空间中需要 1,000 个点。随着维数的增加，填充空间所需的实验次数呈指数增长，因为空间的维数增加。这种指数增长被称为“维数灾难”。
>>> import numpy as np
>>> disc = 10
>>> x1 = np.linspace(0, 1, disc)
>>> x2 = np.linspace(0, 1, disc)
>>> x3 = np.linspace(0, 1, disc)
>>> x1, x2, x3 = np.meshgrid(x1, x2, x3)
为了缓解这个问题，已经设计了 QMC 方法。它们是确定性的，对空间有很好的覆盖率，其中一些可以继续并保持良好的特性。与 MC 方法的主要区别在于，这些点不是 IID，而是知道先前点的。因此，一些方法也被称为序列。
此图显示了 2 组 256 个点。左侧的设计是普通 MC，而右侧的设计是使用 *Sobol'* 方法的 QMC 设计。我们清楚地看到，QMC 版本更加均匀。这些点在边界附近更好地采样，并且集群或间隙更少。
评估均匀性的一种方法是使用称为差异的度量。这里 *Sobol'* 点的差异优于粗略的 MC。
回到均值的计算，QMC 方法对误差也有更好的收敛速度。它们可以为该函数实现 \(O(n^{-1})\)，并且在非常平滑的函数上甚至可以实现更好的速度。此图显示 *Sobol'* 方法的速度为 \(O(n^{-1})\)
有关更多数学细节，请参阅 scipy.stats.qmc 的文档。
计算差异#
让我们考虑两组点。从下图可以清楚地看出，左侧的设计比右侧的设计覆盖了更多的空间。这可以使用 discrepancy 度量来量化。差异越低，样本越均匀。
>>> import numpy as np
>>> from scipy.stats import qmc
>>> space_1 = np.array([[1, 3], [2, 6], [3, 2], [4, 5], [5, 1], [6, 4]])
>>> space_2 = np.array([[1, 5], [2, 4], [3, 3], [4, 2], [5, 1], [6, 6]])
>>> l_bounds = [0.5, 0.5]
>>> u_bounds = [6.5, 6.5]
>>> space_1 = qmc.scale(space_1, l_bounds, u_bounds, reverse=True)
>>> space_2 = qmc.scale(space_2, l_bounds, u_bounds, reverse=True)
>>> qmc.discrepancy(space_1)
0.008142039609053464
>>> qmc.discrepancy(space_2)
0.010456854423869011
使用 QMC 引擎#
实现了多个 QMC 采样器/引擎。这里我们看一下两种最常用的 QMC 方法：Sobol 和 Halton 序列。
"""Sobol' and Halton sequences."""
from scipy.stats import qmc
import numpy as np
import matplotlib.pyplot as plt
rng = np.random.default_rng()
n_sample = 256
dim = 2
sample = {}
# Sobol'
engine = qmc.Sobol(d=dim, seed=rng)
sample["Sobol'"] = engine.random(n_sample)
# Halton
engine = qmc.Halton(d=dim, seed=rng)
sample["Halton"] = engine.random(n_sample)
fig, axs = plt.subplots(1, 2, figsize=(8, 4))
for i, kind in enumerate(sample):
axs[i].scatter(sample[kind][:, 0], sample[kind][:, 1])
axs[i].set_aspect('equal')
axs[i].set_xlabel(r'$x_1$')
axs[i].set_ylabel(r'$x_2$')
axs[i].set_title(f'{kind}—$C^2 = ${qmc.discrepancy(sample[kind]):.2}')
plt.tight_layout()
plt.show()
警告
QMC 方法需要特别注意，用户必须阅读文档以避免常见的陷阱。例如，*Sobol'* 需要 2 的幂的点数。此外，稀疏、燃烧或其他点选择可能会破坏序列的特性，并导致一组点，这些点不如 MC 好。
QMC 引擎是状态感知的。这意味着您可以继续序列、跳过一些点或重置它。让我们从 Halton 中获取 5 个点。然后要求第二组 5 个点。
>>> from scipy.stats import qmc
>>> engine = qmc.Halton(d=2)
>>> engine.random(5)
array([[0.22166437, 0.07980522], # random
[0.72166437, 0.93165708],
[0.47166437, 0.41313856],
[0.97166437, 0.19091633],
[0.01853937, 0.74647189]])
>>> engine.random(5)
array([[0.51853937, 0.52424967], # random
[0.26853937, 0.30202745],
[0.76853937, 0.857583 ],
[0.14353937, 0.63536078],
[0.64353937, 0.01807683]])
现在我们重置序列。要求 5 个点会导致相同的最初 5 个点。
>>> engine.reset()
>>> engine.random(5)
array([[0.22166437, 0.07980522], # random
[0.72166437, 0.93165708],
[0.47166437, 0.41313856],
[0.97166437, 0.19091633],
[0.01853937, 0.74647189]])
在这里，我们推进序列以获得相同的第二组 5 个点。
>>> engine.reset()
>>> engine.fast_forward(5)
>>> engine.random(5)
array([[0.51853937, 0.52424967], # random
[0.26853937, 0.30202745],
[0.76853937, 0.857583 ],
[0.14353937, 0.63536078],
[0.64353937, 0.01807683]])
创建 QMC 引擎，即子类化 QMCEngine#
要创建自己的 QMCEngine，必须定义一些方法。以下是一个包装 numpy.random.Generator 的示例。
>>> import numpy as np
>>> from scipy.stats import qmc
>>> class RandomEngine(qmc.QMCEngine):
... def __init__(self, d, seed=None):
... super().__init__(d=d, seed=seed)
... self.rng = np.random.default_rng(self.rng_seed)
...
...
... def _random(self, n=1, *, workers=1):
... return self.rng.random((n, self.d))
...
...
... def reset(self):
... self.rng = np.random.default_rng(self.rng_seed)
... self.num_generated = 0
... return self
...
...
... def fast_forward(self, n):
... self.random(n)
... return self
然后我们像使用其他任何 QMC 引擎一样使用它。
>>> engine = RandomEngine(2)
>>> engine.random(5)
array([[0.22733602, 0.31675834], # random
[0.79736546, 0.67625467],
[0.39110955, 0.33281393],
[0.59830875, 0.18673419],
[0.67275604, 0.94180287]])
>>> engine.reset()
>>> engine.random(5)
array([[0.22733602, 0.31675834], # random
[0.79736546, 0.67625467],
[0.39110955, 0.33281393],
[0.59830875, 0.18673419],
[0.67275604, 0.94180287]])