SciPy 中的通用非均匀随机数采样#
SciPy 提供了一系列通用非均匀随机数生成器的接口，用于从各种一元连续和离散分布中抽取随机变量。为了保证速度和性能，系统使用了名为 UNU.RAN 的快速 C 库实现。如需深入了解这些方法，请参阅 UNU.RAN 的文档。本教程及所有生成器的文档在编写过程中都大量参考了该文档。
简介#
随机变量生成是一个研究如何从各种分布中生成随机变量算法的小型研究领域。通常假设已知一个均匀随机数生成器。这是一个能够产生一系列独立同分布的连续 U(0,1) 随机变量（即 (0,1) 区间上的均匀随机变量）的程序。当然，现实世界的计算机永远无法生成理想的随机数，也无法产生任意精度的数字，但最先进的均匀随机数生成器已接近这一目标。因此，随机变量生成处理的是如何将此类 U(0,1) 随机数序列转换为非均匀随机变量的问题。这些方法具有通用性，并以黑盒方式运行。
实现此目标的一些方法包括：
逆转法（The Inversion method）：当累积分布函数的逆函数 \(F^{-1}\) 已知时，随机变量生成就变得很简单。我们只需生成一个均匀分布的随机数 U ~ U(0,1)，并返回 \(X = F^{-1}(U)\)。由于很少能获得逆函数的闭式解，通常需要依赖逆函数的近似值（例如
scipy.special.ndtri,scipy.special.stdtrit）。通常情况下，与 UNU.RAN 中的逆转法相比，特殊函数的实现速度相当慢。拒绝法（The Rejection Method）：拒绝法通常被称为接受-拒绝法，由约翰·冯·诺依曼于 1951 年提出[1]。它涉及计算 PDF 的上界（也称为帽子函数），并使用逆转法从该上界生成一个随机变量，假设为 Y。然后，在 0 到上界在 Y 处的值之间抽取一个均匀随机数。如果该数字小于 Y 处的 PDF，则返回该样本，否则予以拒绝。参见
scipy.stats.sampling.TransformedDensityRejection。均匀比率法（The Ratio-of-Uniforms Method）：这是一种接受-拒绝法，它使用最小边界矩形来构建帽子函数。参见
scipy.stats.sampling.RatioUniforms。离散分布的逆转法：与连续情况的区别在于 \(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 参数对分布进行截断。所有生成器都需要一个均匀随机数流，这些随机数会被转换为给定分布的随机变量。这是通过传递带有 NumPy BitGenerator 作为均匀随机数生成器的 random_state 参数来完成 combat 的。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
import numpy as np
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)
如示例所示，我们首先初始化一个分布对象，该对象包含生成器所需方法的实现。在我们的例子中，我们使用变换密度拒绝 (TDR) 方法，该方法需要 PDF 及其相对于 x（即变量）的导数。
注意
请注意，分布的方法（即 pdf、dpdf 等）不需要向量化。它们应该接受并返回浮点数。
注意
也可以将 SciPy 分布作为参数传递。但是，请注意，该对象并不总是包含某些生成器（如 TDR 方法所需的 PDF 导数）所需的所有信息。由于 pdf 和 cdf 等方法的向量化，依赖 SciPy 分布也可能会降低性能。在这两种情况下，都可以实现一个自定义分布对象，其中包含所有必需的方法且不进行向量化，如上例所示。如果想对 SciPy 中定义的分布应用数值逆转法，请参考 scipy.stats.sampling.FastGeneratorInversion。
在上面的示例中，我们设置了一个 TransformedDensityRejection 方法的对象，用于从标准正态分布中采样。现在，我们可以通过调用 rvs 方法开始从我们的分布中采样
rng.rvs()
0.7911236383287276
rng.rvs((5, 3))
array([[ 0.94717699, -0.18902938, -0.08389448],
[-0.90655761, 0.09918658, 0.13561142],
[-2.01232016, -0.20579 , 0.23439292],
[-1.22889898, -0.18955583, 1.12464389],
[ 0.52648626, -0.47670397, 0.95772632]])
我们还可以通过可视化样本的直方图来检查样本是否是从正确的分布中抽取的
import numpy as np
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 来处理已知分布（例如，正态分布、贝塔分布）和其它分布的变换（例如，正态逆高斯分布 scipy.stats.norminvgauss 和对数正态分布 scipy.stats.lognorm）。如果没有实现特定方法，scipy.stats.rv_continuous 默认使用非常慢的 CDF 数值逆转法。由于 UNU.RAN 转换均匀随机数的方式与 SciPy 或 NumPy 不同，因此即使对于相同的均匀随机数流，生成的 RV（随机变量）流也是不同的。例如，即使对于相同的 random_state，SciPy 的 scipy.stats.norm 和 UNU.RAN 的 TransformedDensityRejection 的随机数流也不会相同
from scipy.stats.sampling import norm, 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.14146894, -0.37218024, 0.04759399],
[-0.72646744, 0.96456598, 0.14687657],
[-0.59245914, -0.08711739, 0.45031751],
[ 0.27099605, -0.00765304, -0.7415469 ],
[-0.52377992, 0.24028669, 0.9135163 ]])
无效和错误的参数由 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 抛出的警告也遵循相同的格式。