以 notebook 形式打开

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)\)。由于反函数的封闭形式解很少可用，通常需要依赖反函数的近似值（例如 scipy.special.ndtri, scipy.special.stdtrit)。一般来说，与 UNU.RAN 中的反演方法相比，特殊函数的实现相当慢。

• 拒绝法：拒绝法，通常称为接受-拒绝法，由 John von Neumann 于 1951 年提出[1]。它包括计算 PDF 的上限（也称为帽子函数），并使用反演法从此上限生成一个随机变量，例如 Y。然后可以在 0 到 Y 处上限值之间绘制一个均匀随机数。如果此数字小于 Y 处的 PDF，则返回样本，否则拒绝它。参见 scipy.stats.sampling.TransformedDensityRejection。

• 均匀比率法：这是一种接受-拒绝法，它使用最小边界矩形来构造帽子函数。参见 scipy.stats.sampling.RatioUniforms。

• 离散分布的反演：与连续情况相比，区别在于 \(F\) 现在是一个阶跃函数。为了在计算机中实现这一点，使用了一种搜索算法，其中最简单的是顺序搜索。从 U(0, 1) 生成一个均匀随机数，并对概率求和，直到累积概率超过均匀随机数。发生这种情况的索引是所需的随机变量，并返回。

有关这些算法的更多详细信息，请参阅 UNU.RAN 用户手册的附录。

在生成分布的随机变量时，有两个因素对于确定生成器的速度非常重要：设置步骤和实际采样。根据情况，不同的生成器可能是最佳的。例如，如果需要从具有固定形状参数的给定分布中重复抽取大量样本，则如果采样速度快，则可以接受缓慢的设置。这称为固定参数情况。如果目标是为不同的形状参数（变化参数情况）生成分布样本，则需要为每个参数重复进行昂贵的设置会导致性能非常差。在这种情况下，快速设置对于实现良好的性能至关重要。下表显示了不同方法的设置和采样速度的概述。

连续分布的方法	必需输入	可选输入	设置速度	采样速度
TransformedDensityRejection	pdf，dpdf	无	慢	快
scipy.stats.sampling.NumericalInverseHermite	cdf	pdf，dpdf	（非常）慢	（非常）快
scipy.stats.sampling.NumericalInversePolynomial	pdf	cdf	（非常）慢	（非常）快
scipy.stats.sampling.SimpleRatioUniforms	pdf	无	快	慢

其中

• pdf：概率密度函数

• dpdf：pdf 的导数

• cdf：累积分布函数

要以最少的精力将数值反演方法 NumericalInversePolynomial 应用于 SciPy 中的大量连续分布，请查看 scipy.stats.sampling.FastGeneratorInversion。

离散分布的方法	必需输入	可选输入	设置速度	采样速度
scipy.stats.sampling.DiscreteAliasUrn	pv	pmf	慢	非常快
scipy.stats.sampling.DiscreteGuideTable	pv	pmf	慢	非常快

其中

• pv：概率向量

• pmf：概率质量函数

有关 UNU.RAN 中实现的生成器的更多详细信息，请参阅 [2] 和 [3]。

接口的基本概念

每个生成器都需要先设置好，然后才能开始从中采样。这可以通过实例化该类的对象来完成。大多数生成器都以分布对象作为输入，其中包含所需方法的实现，例如 PDF、CDF 等。除了分布对象之外，还可以传递用于设置生成器的参数。也可以使用 domain 参数截断分布。所有生成器都需要一串均匀随机数，这些随机数被转换为给定分布的随机变量。这通过传递带有 NumPy BitGenerator 的 random_state 参数作为均匀随机数生成器来实现。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)

如示例所示，我们首先初始化一个分布对象，其中包含生成器所需方法的实现。在我们的例子中，我们使用 TransformedDensityRejection (TDR) 方法，该方法需要 PDF 及其相对于 x（即变量）的导数。

注意

请注意，分布的方法（即 pdf、dpdf 等）不需要向量化。它们应该接受并返回浮点数。

注意

也可以将 SciPy 分布作为参数传递。但是，请注意，该对象并不总是具有某些生成器（如 TDR 方法的 PDF 导数）所需的所有信息。依赖 SciPy 分布也可能会降低性能，因为诸如 pdf 和 cdf 之类的方法会进行矢量化处理。在这两种情况下，都可以实现一个自定义分布对象，该对象包含所有必需的方法，并且不会像上面的示例中那样进行矢量化处理。如果想对 SciPy 中定义的分布应用数值反演方法，请参阅 scipy.stats.sampling.FastGeneratorInversion。

在上面的示例中，我们设置了 TransformedDensityRejection 方法的对象，以便从标准正态分布中采样。现在，我们可以通过调用 rvs 方法开始从我们的分布中采样。

rng.rvs()

-1.0494228051028878

rng.rvs((5, 3))

array([[-0.79011734, -0.97364192,  1.30965235],
       [ 1.77777466, -0.73648314,  0.20012238],
       [ 0.28720556, -0.26408585, -0.01241964],
       [-1.43792314,  0.60251365,  1.44168353],
       [-0.67256502, -1.73323039, -1.30021512]])

我们还可以通过可视化样本的直方图来检查样本是否是从正确的分布中提取的。

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 来处理众所周知的分布（例如，正态分布、beta 分布）和其他分布的变换（例如，正态逆高斯分布 scipy.stats.norminvgauss 和对数正态分布 scipy.stats.lognorm）。如果没有实现特定的方法，scipy.stats.rv_continuous 将默认使用 CDF 的数值反演方法，这种方法非常慢。由于 UNU.RAN 以不同于 SciPy 或 NumPy 的方式转换均匀随机数，即使对于相同的均匀随机数流，生成的随机变量流也是不同的。例如，即使对于相同的 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.34700433,  0.78860855,  0.25656489],
       [-0.12975904,  0.20558451,  0.98210181],
       [-0.28086258,  0.75469075, -0.32780586],
       [ 0.1158423 ,  0.09094358, -0.30694274],
       [-0.95995934,  0.80675682,  0.19385288]])

无效和错误的参数将由 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 抛出的警告也遵循相同的格式。

scipy.stats.sampling 中的生成器

• 转换密度拒绝 (TDR)
• 离散别名瓮 (DAU)
• 基于多项式插值的 CDF 反演 (PINV)
• 离散引导表 (DGT)
• 基于埃尔米特插值的 CDF 反演 (HINV)
• 简单均匀比率 (SROU)

参考资料

[1]

Von Neumann, John. “13. various techniques used in connection with random digits.” Appl. Math Ser 12.36-38 (1951): 3. https://mcnp.lanl.gov/pdf_files/InBook_Computing_1961_Neumann_JohnVonNeumannCollectedWorks_VariousTechniquesUsedinConnectionwithRandomDigits.pdf

[2]

UNU.RAN 用户手册，https://statmath.wu.ac.at/unuran/doc/unuran.html

[3]

Leydold, Josef, Wolfgang Hörmann 和 Halis Sak。“An R Interface to the UNU.RAN Library for Universal Random Variate Generators.”，https://cran.r-project.org.cn/web/packages/Runuran/vignettes/Runuran.pdf