RatioUniforms#
- class scipy.stats.sampling.RatioUniforms(pdf, *, umax, vmin, vmax, c=0, random_state=None)[source]#
使用比率均匀方法从概率密度函数生成随机样本。
- 参数:
- pdf可调用
签名函数为 pdf(x) 的函数,它与分布的概率密度函数成正比。
- umax浮点数
u 方向定界矩形的上限。
- vmin浮点数
v 方向定界矩形的下限。
- vmax浮点数
v 方向定界矩形的上限。
- c浮点数,可选。
比率均匀方法的位移参数,详情见注释。默认为 0。
- random_state{None, int,
numpy.random.Generator
, numpy.random.RandomState
}, 默认为空如果seed 为 None (或 np.random),将使用
numpy.random.RandomState
单例。如果seed 为 int,将使用新的RandomState
实例,其中 seed 为种子。如果seed 已是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)
既在R
中,亦在A
中,那么(U, V)
位于A
中的可能性会很高,但如果不在A
中,就需要进行过多的迭代。更确切一点讲,请注意,要想画出(U, V)
,使之均匀分布于R
,又使(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 和 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])随机变量的采样