scipy.stats.sampling.

RatioUniforms#

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

使用比率均匀方法从概率密度函数生成随机样本。

参数:
pdf可调用对象

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

umax浮点数

u 方向上边界矩形的上界。

vmin浮点数

v 方向上边界矩形的下界。

vmax浮点数

v 方向上边界矩形的上界。

c浮点数,可选。

比率均匀方法的偏移参数,请参阅注释。默认为 0。

random_state{None, int, numpy.random.Generator,

如果 seed 为 None(或 np.random),则使用 numpy.random.RandomState 单例。如果 seed 为整数,则使用新的 RandomState 实例,并使用 seed 进行初始化。如果 seed 已经是 GeneratorRandomState 实例,则使用该实例。

注释

给定一个单变量概率密度函数 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) 有界(即次二次尾部),则这些值是有限的。可以在 R 上均匀生成 (U, V),如果 (U, V) 也在 A 中,则返回 V/U + c,这可以直接验证。

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

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

如果边界矩形未正确指定(即,如果不包含 A),则算法会从与 pdf 给出的分布不同的分布中采样。因此,建议执行诸如 kstest 之类的测试作为检查。

参考文献

[1] (1,2)

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

[2]

W. Hoermann 和 J. Leydold,“生成广义逆高斯随机变量”,《统计和计算》,24(4),第 547–557 页,2014 年。

[3]

A.J. Kinderman 和 J.F. Monahan,“使用均匀偏差比率的随机变量计算机生成”,《ACM 数学软件汇刊》,3(3),第 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

指数分布提供了另一个可以显式确定边界矩形的示例。

>>> 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

方法

rvs([size])

随机变量采样