通用非均匀随机数抽样在SciPy中的应用¶
SciPy提供了到许多通用非均匀随机数生成器的接口，以从各种单变量连续和离散分布中采样随机变量。名为的快速C库的实现 UNU.RAN 用于速度和性能。请看 UNU.RAN's documentation 以深入解释这些方法。它在编写本教程和所有生成器的文档时被大量引用。
引言¶
随机变量生成是一个很小的研究领域，它涉及从各种分布生成随机变量的算法。通常假设有统一的随机数生成器可用。这是一个产生独立且同分布的连续U(0，1)随机变量(即区间(0，1)上的均匀随机变量)序列的程序。当然，现实世界的计算机永远不能产生理想的随机数，它们也不能产生任意精度的数字，但最先进的均匀随机数生成器接近这一目标。因此，随机变量生成处理将这样的U(0，1)随机数序列转换为非均匀随机变量的问题。这些方法是通用的，并且以黑盒方式工作。
执行此操作的一些方法包括：
反演方法：当反演时 \(F^{{-1}}\) 累积分布函数的大小是已知的，那么随机变量的产生就很容易。我们只生成一个U(0，1)均匀分布的随机数U并返回 \(X = F^{{-1}}(U)\) 。由于逆解的闭合形式解很少可用，人们通常需要依赖逆解的近似(例如，
ndtri，stdtrit)。一般来说，与UNU.RAN的反演方法相比，特殊功能的实现相当缓慢。拒绝法：拒绝法，通常称为接受-拒绝法，由约翰·冯·诺伊曼于1951年提出。 1. 它包括计算PDF的上限(也称为HAT函数)，并使用反演方法从该上限生成随机变量，比如Y。然后可以在0到Y的上界的值之间画一个统一的随机数。如果这个数字小于Y处的PDF，则返回样本，否则拒绝该样本。看见
TransformedDensityRejection。一致性比法：这是一种接受拒绝法，它使用最小外接矩形来构造帽子函数。看见
rvs_ratio_uniforms。离散分布的反演：与连续情况不同的是 \(F\) 现在是阶跃函数。为了在计算机中实现这一点，使用了搜索算法，其中最简单的算法是 顺序搜索 。从U(0，1)生成均匀随机数，并对概率求和，直到累积概率超过均匀随机数。发生这种情况的索引是所需的随机变量，并被返回。
有关这些算法的更多详细信息，请参阅 appendix of the UNU.RAN user manual 。
在生成分布的随机变量时，决定生成器速度的两个因素很重要：设置步骤和实际抽样。根据不同的情况，不同的发电机可以是最优的。例如，如果需要重复使用固定形状参数从给定分布中提取大样本，则如果采样速度较快，则可以接受较慢的设置。这称为固定参数情况。如果目标是为不同形状参数(可变参数情况)生成分布的样本，则需要为每个参数重复的昂贵设置将导致非常差的性能。在这种情况下，快速设置对于实现良好的性能至关重要。下表概述了不同方法的设置和采样速度。
连续分布的几种方法 |
所需输入 |
可选输入 |
设置速度 |
采样速度 |
|---|---|---|---|---|
PDF，dpdf |
无 |
慢的 |
快地 |
|
CDF |
无 |
(非常)慢 |
(非常)快 |
|
PDF格式 |
CDF |
(非常)慢 |
(非常)快 |
|
PDF格式 |
多种多样 |
慢/快 |
中等 |
哪里
PDF：概率密度函数
dpdf：pdf的派生
CDF：累积分布函数
离散分布的方法 |
所需输入 |
可选输入 |
设置速度 |
采样速度 |
|---|---|---|---|---|
光伏发电 |
PMF |
慢的 |
非常快 |
哪里
PV：概率向量
PMF：概率质量函数
接口的基本概念¶
每台发电机都需要安装好，然后才能开始采样。这可以通过实例化该类的对象来完成。大多数生成器以包含PDF、CDF等所需方法的实现的分布对象作为输入。除了分布对象，还可以传递用于设置生成器的参数。还可以使用 domain 参数。所有的生成器都需要一个均匀的随机数流，这些随机数被转换成给定分布的随机变量。这是通过将一个 random_state 参数，并将NumPy BitGenerator用作统一随机数生成器。 random_state 可以是整数， np.random.Generator ，或 np.random.RandomState 。
警告
不鼓励使用NumPy<1.19.0，因为它没有快速的Cython API来生成统一的随机数，并且对于实际使用来说可能太慢了。
所有的发电机都有一个共同的 rvs 方法，该方法可用于从给定分布中提取样本。
此界面的示例如下所示：
>>> from scipy.stats 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)
如示例所示，我们首先初始化一个分发对象，该对象包含生成器所需方法的实现。在我们的示例中，我们使用 TransformedDensityRejection (TDR)方法，需要PDF及其导数w.r.t. x (即变量)。
注解
请注意，分发的方法(即 pdf ， dpdf 等)不需要矢量化。他们应该接受并返还浮动资金。
注解
还可以将SciPy分发作为参数传递，但由于验证和昂贵的NumPy操作，传递速度可能会很慢。此外，它并不总是包含某些生成器(如TDR方法的PDF导数)所需的全部信息。此外，SciPy中的大多数发行版都提供了 rvs 方法，该方法可以替代使用。
在上面的示例中，我们设置了 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, 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()
注解
请注意两者之间的区别 rvs 中出现的分布的方法 scipy.stats 还有这些发电机提供的那个。UNU.RAN生成器必须被认为是独立的，因为它们通常会生成不同于 scipy.stats 对于任何种子。该计划的实施 rvs 在……里面 rv_continuous 通常依赖于NumPy模块 np.random 对于众所周知的分布(例如，对于正态分布、贝塔分布)和其他分布的变换(例如，正态逆高斯 norminvgauss 而对数正态分布 lognorm 分布)。如果没有实现特定的方法， rv_continuous 默认采用CDF的数值反演方法，速度非常慢。由于UNU.RAN变换均匀随机数的方式不同于SciPy或NumPy，因此即使对于相同的均匀随机数流，产生的RV流也是不同的。例如，SciPy的随机数流 ~norm 和UNU.RAN的 TransformedDensityRejection 即使是为了同样的目的也不会是一样的 random_state ：
>>> from scipy.stats 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.99865691, 0.38104014, 0.31633526],
[ 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> shows what caused the error. This, by itself, should contain enough information to help debug the error. In addition, <error id> and <type of error> can be used to investigate different classes of error in UNU.RAN. A complete list of all the error codes and their descriptions can be found in the Section 8.4 of the UNU.RAN user manual 。
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¶
参考文献¶
- 1
冯·诺伊曼，约翰。“13.与随机数字有关的各种技术。”应用程序。数学系列12.36-38(1951)：3.
- 2
- 3
题名/责任者：Reydold，Josef，Wolfgang Hörmann，和Halis Sak。“用于通用随机变量生成器的联合国类库的R接口。”，https://cran.r-project.org/web/packages/Runuran/vignettes/Runuran.pdf