SciPy 中的通用非均匀随机数采样#
SciPy 提供了许多通用非均匀随机数生成器的接口，用于从各种单变量连续和离散分布中采样随机变量。使用名为 UNU.RAN 的快速 C 库的实现来提高速度和性能。有关这些方法的深入解释，请参阅 UNU.RAN 的文档。本教程和所有生成器的文档都大量引用了它。
介绍#
随机变量生成是研究的一个小领域，它处理从各种分布中生成随机变量的算法。通常假设有一个均匀随机数生成器。这是一个程序，它产生一系列独立同分布的连续 U(0,1) 随机变量（即区间 (0,1) 上的均匀随机变量）。当然，现实世界的计算机永远无法生成理想的随机数，也无法生成任意精度的数字，但最先进的均匀随机数生成器可以接近这一目标。因此，随机变量生成处理将这种 U(0,1) 随机数序列转换为非均匀随机变量的问题。这些方法是通用的，以黑盒方式工作。
一些执行此操作的方法是
逆方法：当累积分布函数的逆 \(F^{-1}\) 已知时，随机变量生成很容易。我们只需生成一个均匀分布的 U(0,1) 随机数 U 并返回 \(X = F^{-1}(U)\)。由于很少有逆的封闭形式解可用，因此通常需要依赖逆的近似值（例如，
ndtri、stdtrit）。通常，与 UNU.RAN 中的逆方法相比，特殊函数的实现相当缓慢。拒绝方法：拒绝方法，通常称为接受-拒绝方法，由 John von Neumann 在 1951 年提出 [1]。它涉及计算 PDF 的上限（也称为帽函数），并使用逆方法从该上限生成一个随机变量，例如 Y。然后可以在 0 到 Y 处上限的值之间绘制一个均匀随机数。如果该数字小于 Y 处的 PDF，则返回样本，否则拒绝它。请参见
TransformedDensityRejection。均匀比方法：这是一种接受-拒绝方法，它使用最小包围矩形来构建帽函数。请参见
scipy.stats.rvs_ratio_uniforms。离散分布的逆：与连续情况的区别在于 \(F\) 现在是一个阶跃函数。为了在计算机中实现这一点，使用了一种搜索算法，其中最简单的是顺序搜索。从 U(0, 1) 生成一个均匀随机数，并对概率求和，直到累积概率超过均匀随机数。发生这种情况的索引是所需的随机变量，并被返回。
有关这些算法的更多详细信息，请参见 UNU.RAN 用户手册的附录。
在生成分布的随机变量时，有两个因素决定生成器的速度：设置步骤和实际采样。根据情况，不同的生成器可能是最佳的。例如，如果需要从具有固定形状参数的给定分布中重复绘制大型样本，则如果采样速度快，则慢速设置是可以接受的。这称为固定参数情况。如果目标是为不同的形状参数（可变参数情况）生成分布的样本，则需要为每个参数重复的昂贵设置会导致非常糟糕的性能。在这种情况下，快速的设置对于实现良好的性能至关重要。下表概述了不同方法的设置和采样速度。
连续分布的方法 |
必需输入 |
可选输入 |
设置速度 |
采样速度 |
|---|---|---|---|---|
pdf、dpdf |
无 |
慢 |
快 |
|
cdf |
pdf、dpdf |
（非常）慢 |
（非常）快 |
|
cdf |
（非常）慢 |
（非常）快 |
||
无 |
快 |
慢 |
其中
pdf：概率密度函数
dpdf：pdf 的导数
cdf：累积分布函数
要将数值反演方法 NumericalInversePolynomial 应用于 SciPy 中的大量连续分布，并以最小的努力，请查看 scipy.stats.sampling.FastGeneratorInversion。
离散分布的方法 |
必需输入 |
可选输入 |
设置速度 |
采样速度 |
|---|---|---|---|---|
pv |
pmf |
慢 |
非常快 |
|
pv |
pmf |
慢 |
非常快 |
其中
pv：概率向量
pmf：概率质量函数
接口的基本概念#
每个生成器都需要在开始从它采样之前进行设置。这可以通过实例化该类的对象来完成。大多数生成器以分布对象作为输入，该对象包含所需方法（如 PDF、CDF 等）的实现。除了分布对象之外，还可以传递用于设置生成器的参数。还可以使用 domain 参数截断分布。所有生成器都需要一个均匀随机数流，这些流将被转换为给定分布的随机变量。这是通过传递 random_state 参数，并使用 NumPy BitGenerator 作为均匀随机数生成器来完成的。 random_state 可以是整数、numpy.random.Generator 或 numpy.random.RandomState。
警告
不鼓励使用 NumPy < 1.19.0，因为它没有用于生成均匀随机数的快速 Cython API，对于实际使用来说可能太慢。
所有生成器都有一个通用的 rvs 方法，可用于从给定分布中绘制样本。
以下示例显示了此接口
>>> from scipy.stats.sampling import TransformedDensityRejection
>>> from math import exp
>>>
>>> class StandardNormal:
... def pdf(self, x: float) -> float:
... # note that the normalization constant isn't required
... return exp(-0.5 * x*x)
... def dpdf(self, x: float) -> float:
... return -x * exp(-0.5 * x*x)
...
>>> dist = StandardNormal()
>>>
>>> import numpy as np
>>> urng = np.random.default_rng()
>>> rng = TransformedDensityRejection(dist, random_state=urng)
如示例所示，我们首先初始化一个分布对象，该对象包含生成器所需方法的实现。在我们的例子中，我们使用 TransformedDensityRejection (TDR) 方法，该方法需要一个 PDF 及其相对于 x（即变量）的导数。
注意
请注意，分布的方法（即 pdf、dpdf 等）不需要是矢量化的。它们应该接受和返回浮点数。
注意
还可以将 SciPy 分布作为参数传递。但是，请注意，该对象并不总是包含某些生成器所需的所有信息，例如 TDR 方法的 PDF 导数。依赖 SciPy 分布也可能会降低性能，因为方法（如 pdf 和 cdf）的矢量化。在这两种情况下，都可以实现一个自定义分布对象，该对象包含所有必需的方法，并且未被矢量化，如上面的示例所示。如果想将数值反演方法应用于 SciPy 中定义的分布，请也查看 scipy.stats.sampling.FastGeneratorInversion。
在上面的示例中，我们设置了 TransformedDensityRejection 方法的对象，以从标准正态分布中采样。现在，我们可以通过调用 rvs 方法来开始从我们的分布中采样
>>> rng.rvs()
-1.526829048388144
>>> rng.rvs((5, 3))
array([[ 2.06206883, 0.15205036, 1.11587367],
[-0.30775562, 0.29879802, -0.61858268],
[-1.01049115, 0.78853694, -0.23060766],
[-0.60954752, 0.29071797, -0.57167182],
[ 0.9331694 , -0.95605208, 1.72195199]])
我们也可以通过可视化样本的直方图来检查样本是否来自正确的分布。
>>> import matplotlib.pyplot as plt
>>> from scipy.stats import norm
>>> from scipy.stats.sampling import TransformedDensityRejection
>>> from math import exp
>>>
>>> class StandardNormal:
... def pdf(self, x: float) -> float:
... # note that the normalization constant isn't required
... return exp(-0.5 * x*x)
... def dpdf(self, x: float) -> float:
... return -x * exp(-0.5 * x*x)
...
>>>
>>> dist = StandardNormal()
>>> urng = np.random.default_rng()
>>> rng = TransformedDensityRejection(dist, random_state=urng)
>>> rvs = rng.rvs(size=1000)
>>> x = np.linspace(rvs.min()-0.1, rvs.max()+0.1, num=1000)
>>> fx = norm.pdf(x)
>>> plt.plot(x, fx, 'r-', lw=2, label='true distribution')
>>> plt.hist(rvs, bins=20, density=True, alpha=0.8, label='random variates')
>>> plt.xlabel('x')
>>> plt.ylabel('PDF(x)')
>>> plt.title('Transformed Density Rejection Samples')
>>> plt.legend()
>>> plt.show()
注意
请注意 scipy.stats 中的分布的 rvs 方法和这些生成器提供的之间的区别。 UNU.RAN 生成器在某种意义上必须被认为是独立的，这意味着它们通常会产生与 scipy.stats 中等效分布产生的随机数流不同的随机数流，无论种子如何。 scipy.stats.rv_continuous 中 rvs 的实现通常依赖于 NumPy 模块 numpy.random 来实现众所周知的分布（例如，对于正态分布，beta 分布）和其他分布的变换（例如，正态逆高斯 scipy.stats.norminvgauss 和对数正态 scipy.stats.lognorm 分布）。如果没有实现特定的方法，scipy.stats.rv_continuous 默认使用 CDF 的数值反演方法，这非常慢。由于 UNU.RAN 的均匀随机数变换方式不同于 SciPy 或 NumPy，因此即使对于相同的均匀随机数流，产生的 RV 流也不同。例如，SciPy 的 scipy.stats.norm 和 UNU.RAN 的 TransformedDensityRejection 的随机数流即使对于相同的 random_state 也不相同。
>>> from scipy.stats import norm
>>> from scipy.stats.sampling import TransformedDensityRejection
>>> from copy import copy
>>> dist = StandardNormal()
>>> urng1 = np.random.default_rng()
>>> urng1_copy = copy(urng1)
>>> rng = TransformedDensityRejection(dist, random_state=urng1)
>>> rng.rvs()
-1.526829048388144
>>> norm.rvs(random_state=urng1_copy)
1.3194816698862635
我们可以传递 domain 参数来截断分布。
>>> rng = TransformedDensityRejection(dist, domain=(-1, 1), random_state=urng)
>>> rng.rvs((5, 3))
array([[-0.99865691, 0.38104014, 0.31633526], # may vary
[ 0.88433909, -0.45181849, 0.78574461],
[ 0.3337244 , 0.12924307, 0.40499404],
[-0.51865761, 0.43252222, -0.6514866 ],
[-0.82666174, 0.71525582, 0.49006743]])
无效参数和错误参数由 SciPy 或 UNU.RAN 处理。后者会抛出 UNURANError，该错误遵循通用格式。
UNURANError: [objid: <object id>] <error code>: <reason> => <type of error>
其中
<object id>是 UNU.RAN 给出的对象的 ID。<error code>是一个错误代码，代表一种错误类型。<reason>是错误发生的原因。<type of error>是对错误类型的简短描述。
<reason> 显示了导致错误的原因。这本身应该包含足够的信息来帮助调试错误。此外，<error id> 和 <type of error> 可用于调查 UNU.RAN 中的不同错误类别。所有错误代码及其描述的完整列表可以在 UNU.RAN 用户手册的第 8.4 节 中找到。
下面显示了 UNU.RAN 生成的错误示例。
UNURANError: [objid: TDR.003] 50 : PDF(x) < 0.! => (generator) (possible) invalid data
这表明 UNU.RAN 无法初始化 ID 为 TDR.003 的对象，因为 PDF < 0，即为负。这属于“生成器可能无效数据”类型，错误代码为 50。
UNU.RAN 抛出的警告也遵循相同的格式。