RatioUniforms#
- class scipy.stats.sampling.RatioUniforms(pdf, *, umax, vmin, vmax, c=0, random_state=None)[source]#
使用均匀比方法从概率密度函数生成随机样本。
- 参数:
- pdfcallable
一个具有签名 pdf(x) 的函数,它与分布的概率密度函数成比例。
- umaxfloat
u 方向上边界矩形的上限。
- vminfloat
v 方向上边界矩形的下限。
- vmaxfloat
v 方向上边界矩形的上限。
- cfloat,可选。
均匀比方法的移位参数,参见 Notes。默认为 0。
- random_state{None, int,
numpy.random.Generator
, 如果 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)
,如果(U, V)
也在A
中,则返回V/U + c
,可以直接验证。如果对于任何常数 k > 0,用 k * pdf 替换 pdf,则该算法不会改变。 因此,通常通过删除不必要的归一化因子来使用与概率密度函数成比例的函数。
直观上,如果
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
之类的测试作为检查。参考文献
[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 and J.F. Monahan, “Computer Generation of Random Variables Using the Ratio of Uniform Deviates”, 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
指数分布提供了另一个可以明确确定边界矩形的示例。
>>> 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