CQF预备知识：Python相关库 -- 通用非均匀随机数抽样 scipy.stats

文中内容仅限技术学习与代码实践参考，市场存在不确定性，技术分析需谨慎验证，不构成任何投资建议。

通用非均匀随机数抽样

SciPy 提供了许多通用非均匀随机数生成器的接口，用于从广泛的单变量连续和离散分布中抽取随机变量。这些实现基于一个快速的 C 库 UNU.RAN，以提高速度和性能。有关这些方法的深入解释，请参阅 UNU.RAN 的文档。编写本教程和所有生成器的文档时，也大量参考了 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 中的反演方法慢得多。

• 拒绝法：拒绝法（也称为接受-拒绝法）是由约翰·冯·诺伊曼在 1951 ¹年提出的。它涉及计算概率密度函数（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：累积分布函数

为了在 SciPy 中以最小的努力将数值反演方法 NumericalInversePolynomial 应用于大量连续分布，可以参考 scipy.stats.sampling.FastGeneratorInversion。

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

其中：

• pv：概率向量

• pmf：概率质量函数

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

接口的基本概念

在开始抽样之前，每个生成器都需要进行设置。这可以通过实例化该类的对象来完成。大多数生成器需要一个分布对象作为输入，该对象包含实现所需方法（如 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 npclass StandardNormal:def pdf(self, x: float) -> float:# 注意不需要归一化常数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 方法开始从该分布中抽样：

array([[-0.5188201 ,  0.19758648,  1.56394908],[ 1.37478434, -0.25089228, -0.06377169],[ 0.04548633,  0.60338512, -0.18978775],[-0.58243905,  1.5650833 , -0.964274  ],[ 2.41625536, -1.60428538, -0.54605791]])

我们还可以通过绘制样本的直方图来验证样本是否来自正确的分布：

import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import norm
from scipy.stats.sampling import TransformedDensityRejection
from math import expclass StandardNormal:def pdf(self, x: float) -> float:# 注意不需要归一化常数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='真实分布')
plt.hist(rvs, bins=20, density=True, alpha=0.8, label='随机变量')
plt.xlabel('x')
plt.ylabel('PDF(x)')
plt.title('变换密度拒绝法抽样')
plt.legend()
plt.show()

注意

请注意 scipy.stats 中分布的 rvs 方法与这些生成器提供的 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 转换均匀随机数的方式不同，即使对于相同的均匀随机数流，生成的随机变量流也会不同。例如，SciPy 的 scipy.stats.norm 和 UNU.RAN 的 TransformedDensityRejection 即使对于相同的 random_state，其随机数流也不会相同：

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.25668751, -0.29186077, -0.16906123],[ 0.28342675,  0.14713262,  0.44345019],[-0.86292589,  0.40110739, -0.3304392 ],[ 0.18390038,  0.7311804 , -0.98740077],[ 0.67457194, -0.77107968, -0.63685895]])

无效或错误的参数将由 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 User Manual, https://statmath.wu.ac.at/unuran/doc/unuran.html ↩︎

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