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

如果 seed 是 None（或 np.random），则使用 numpy.random.RandomState 单例。如果 seed 是 int，则使用新的 RandomState 实例，并以 seed 为种子。如果 seed 已经是 Generator 或 RandomState 实例，则使用该实例。

方法

rvs([size])

随机变量的采样

注释

给定一个单变量概率密度函数 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) 并返回 V/U + c，如果 (U, V) 也在 A 中，这可以直接验证。

如果将 pdf 替换为任意常数 k > 0 的 k * pdf，算法不会改变。因此，通常使用与概率密度函数成比例的函数来工作是很方便的，通过去掉不必要的归一化因子。

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

如果边界矩形没有正确指定（即，如果它不包含 A），算法会从一个不同于 pdf 给定的分布中采样。因此，建议执行如 kstest 这样的测试作为检查。

参考文献

[1] (1,2)

L. Devroye, “Non-Uniform Random Variate Generation”, Springer-Verlag, 1986.

[2]

W. Hoermann and J. Leydold, “Generating generalized inverse Gaussian random variates”, Statistics and Computing, 24(4), p. 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
On this page