Python SciPy sampling.RatioUniforms用法及代码示例

本文简要介绍 python 语言中 scipy.stats.sampling.RatioUniforms 的用法。

用法:

class  scipy.stats.sampling.RatioUniforms(pdf, *, umax, vmin, vmax, c=0, random_state=None)#

使用 ratio-of-uniforms 方法从概率密度函数生成随机样本。

参数 ：：

pdf： 可调用的

具有签名 pdf(x) 的函数，与分布的概率密度函数成比例。

umax： 浮点数

u-direction 中边界矩形的上限。

vmin： 浮点数

v-direction 中边界矩形的下界。

vmax： 浮点数

v-direction 中边界矩形的上限。

c： 浮点数，可选。

ratio-of-uniforms 方法的移位参数，请参见注释。默认值为 0。

random_state： {无，整数， numpy.random.Generator ，

numpy.random.RandomState}, optional

如果种子是无(或np.random)， 这numpy.random.RandomState使用单例。如果种子是一个 int，一个新的RandomState使用实例，播种种子.如果种子已经是一个Generator或者RandomState实例然后使用该实例。

注意：

给定一个单变量概率密度函数pdf和一个常数c, 定义集合A = {(u, v) : 0 < u <= sqrt(pdf(v/u + c))}.如果(U, V)是一个均匀分布的随机向量A， 然后V/U + c遵循一个分布pdf.

上述结果(见[1],[2]) 可用于仅使用 PDF 对随机变量进行采样，即不需要 CDF 的反转。典型选择c是零或模式pdf.套装A是矩形的子集R = [0, umax] x [vmin, vmax]其中

• umax = sup sqrt(pdf(x))

• vmin = inf (x - c) sqrt(pdf(x))

• vmax = sup (x - c) sqrt(pdf(x))

特别是，这些值是有限的，如果pdf是有界的并且x**2 * pdf(x)是有界的(即次二次尾部)。一个可以生成(U, V)一致地R并返回V/U + c如果(U, V)也在A可以直接验证。

如果对于任何常数 k > 0，将 pdf 替换为 k * pdf，则该算法不会改变。因此，通过删除不必要的归一化因子，使用与概率密度函数成比例的函数通常很方便。

直观地说，如果A填满大部分封闭矩形，使得概率很高(U, V)在于A每当它位于R因为否则所需的迭代次数会变得太大。更准确地说，请注意要绘制的预期迭代次数(U, V)均匀分布在R这样(U, V)也在A由比率给出area(R) / area(A) = 2 * umax * (vmax - vmin) / area(pdf)，其中面积(pdf)是积分pdf(如果使用概率密度函数，则等于一，但如果使用与密度成比例的函数，则可以取其他值)。等式成立，因为面积A等于0.5 * area(pdf)(定理 7.1 中[1])。如果采样在 50000 次迭代后未能生成单个随机变量(即没有一次抽签在A)，引发异常。

如果未正确指定边界矩形(即，如果它不包含A)，算法从不同于下式给出的分布中采样pdf.因此，建议执行测试，例如scipy.stats.kstest作为支票。

参考：

[1] (1,2)

L. Devroye，“非均匀随机变量生成”，Springer-Verlag，1986 年。

[2]

W. Hoermann 和 J. Leydold，“生成广义逆高斯随机变量”，统计与计算，24(4)，p。 547-557，2014 年。

[3]

A.J. Kinderman 和 J.F. Monahan，“使用均匀偏差比率的随机变量的计算机生成”，ACM Transactions on Mathematical Software，3(3)，p。 257-260，1977。

例子：

>>> import numpy as np
>>> from scipy import stats

>>> from scipy.stats.sampling import RatioUniforms
>>> rng = np.random.default_rng()

模拟正态分布的随机变量。在这种情况下，很容易显式地计算边界矩形。为简单起见，我们删除了密度的归一化因子。

>>> f = lambda x: np.exp(-x**2 / 2)
>>> v = np.sqrt(f(np.sqrt(2))) * np.sqrt(2)
>>> umax = np.sqrt(f(0))
>>> gen = RatioUniforms(f, umax=umax, vmin=-v, vmax=v, random_state=rng)
>>> r = gen.rvs(size=2500)

K-S 检验确认随机变量确实是正态分布的(正态性在 5% 显著性水平上不被拒绝)：

>>> stats.kstest(r, 'norm')[1]
0.250634764150542

index 分布提供了另一个可以明确确定边界矩形的示例。

>>> gen = RatioUniforms(lambda x: np.exp(-x), umax=1, vmin=0,
...                     vmax=2*np.exp(-1), random_state=rng)
>>> r = gen.rvs(1000)
>>> stats.kstest(r, 'expon')[1]
0.21121052054580314