离散引导表#
- class scipy.stats.sampling.DiscreteGuideTable(dist, *, domain=None, guide_factor=1, random_state=None)#
离散引导表方法。
离散引导表方法从任意但有限的概率向量中采样。它使用大小为 \(N\) 的概率向量或具有有限支撑的概率质量函数来从分布中生成随机数。离散引导表具有非常缓慢的设置(与向量长度线性相关),但提供非常快速的采样。
- 参数:
- distarray_like 或 object,可选
分布的概率向量 (PV)。如果 PV 不可用,则需要一个具有
pmf
方法的类的实例。PMF 的签名应为:def pmf(self, k: int) -> float
。即它应该接受一个 Python 整数并返回一个 Python 浮点数。- domainint,可选
PMF 的支撑。如果概率向量 (
pv
) 不可用,则必须给出有限域。即 PMF 必须具有有限支撑。默认为None
。当None
如果分布对象 dist 提供了
support
方法,则它用于设置分布的域。否则,假设支撑为
(0, len(pv))
。当此参数与概率向量一起传递时,domain[0]
用于将分布从(0, len(pv))
重定位到(domain[0], domain[0]+len(pv))
,而domain[1]
被忽略。有关更详细的解释,请参阅“注释”和“教程”。
- guide_factor: int,可选
引导表相对于 PV 长度的尺寸。更大的引导表会导致更快的生成时间,但需要更昂贵的设置。不建议使用大于 3 的尺寸。如果相对尺寸设置为 0,则使用顺序搜索。默认为 1。
- random_state{None, int,
numpy.random.Generator
, NumPy 随机数生成器或用于生成均匀随机数流的底层 NumPy 随机数生成器的种子。如果 random_state 为 None(或 np.random),则使用
numpy.random.RandomState
单例。如果 random_state 是一个整数,则使用一个新的RandomState
实例,并使用 random_state 对其进行播种。如果 random_state 已经是Generator
或RandomState
实例,则使用该实例。
注释
当有限概率向量可用或分布的 PMF 可用时,此方法有效。如果仅 PMF 可用,则还必须给出 PMF 的有限支撑(域)。建议首先通过在支撑中的每个点上评估 PMF 来获取概率向量,然后使用它代替。
DGT 从任意但有限的概率向量中采样。随机数是通过反演方法生成的,即
生成一个随机数 U ~ U(0,1)。
找到最小的整数 I 使得 F(I) = P(X<=I) >= U。
步骤 (2) 是关键步骤。使用顺序搜索需要 O(E(X)) 次比较,其中 E(X) 是分布的期望值。但是,索引搜索使用引导表跳到 I 附近的某个 I' <= I,以在恒定时间内找到 X。事实上,当引导表与概率向量的大小相同时,预期的比较次数会减少到 2(这是默认值)。对于更大的引导表,这个数字会更小(但始终大于 1),对于更小的表,这个数字会更大。对于表大小为 1 的极限情况,该算法只是执行顺序搜索。
另一方面,引导表的设置时间为 O(N),其中 N 表示概率向量的长度(对于大小为 1 的情况,不需要预处理)。此外,对于非常大的引导表,内存效应甚至可能降低算法的速度。因此,我们不建议使用比给定概率向量大三倍以上的引导表。如果只需要生成几个随机数,则使用更小的表尺寸更好。引导表相对于给定概率向量长度的尺寸可以通过
guide_factor
参数设置。如果给出了概率向量,它必须是一个一维数组,包含非负浮点数,不包含任何
inf
或nan
值。此外,必须至少有一个非零条目,否则会引发异常。默认情况下,概率向量的索引从 0 开始。但是,可以通过传递
domain
参数来更改此设置。当domain
与 PV 一起传递时,它会将分布从(0, len(pv))
重定位到(domain[0], domain[0] + len(pv))
。在这种情况下,domain[1]
被忽略。参考文献
[1]UNU.RAN 参考手册,第 5.8.4 节,“DGT - (离散)引导表方法(索引搜索)” https://statmath.wu.ac.at/unuran/doc/unuran.html#DGT
[2]H.C. Chen 和 Y. Asau (1974)。关于从经验分布生成随机变量,AIIE Trans. 6,第 163-166 页。
示例
>>> from scipy.stats.sampling import DiscreteGuideTable >>> import numpy as np
要使用概率向量创建随机数生成器,请使用
>>> pv = [0.1, 0.3, 0.6] >>> urng = np.random.default_rng() >>> rng = DiscreteGuideTable(pv, random_state=urng)
RNG 已设置。现在,我们可以使用
rvs
方法从分布中生成样本>>> rvs = rng.rvs(size=1000)
为了验证随机变量是否遵循给定分布,我们可以使用卡方检验(作为拟合优度度量)
>>> from scipy.stats import chisquare >>> _, freqs = np.unique(rvs, return_counts=True) >>> freqs = freqs / np.sum(freqs) >>> freqs array([0.092, 0.355, 0.553]) >>> chisquare(freqs, pv).pvalue 0.9987382966178464
由于 p 值非常高,因此我们不能拒绝观察到的频率与预期频率相同的零假设。因此,我们可以安全地假设这些变量是从给定分布中生成的。请注意,这只是表明算法的正确性,而不是样本的质量。
如果 PV 不可用,也可以传递一个具有 PMF 方法和有限域的类的实例。
>>> urng = np.random.default_rng() >>> from scipy.stats import binom >>> n, p = 10, 0.2 >>> dist = binom(n, p) >>> rng = DiscreteGuideTable(dist, random_state=urng)
现在,我们可以使用
rvs
方法从分布中采样,并测量样本的拟合优度>>> rvs = rng.rvs(1000) >>> _, freqs = np.unique(rvs, return_counts=True) >>> freqs = freqs / np.sum(freqs) >>> obs_freqs = np.zeros(11) # some frequencies may be zero. >>> obs_freqs[:freqs.size] = freqs >>> pv = [dist.pmf(i) for i in range(0, 11)] >>> pv = np.asarray(pv) / np.sum(pv) >>> chisquare(obs_freqs, pv).pvalue 0.9999999999999989
为了检查样本是否已从正确的分布中抽取,我们可以可视化样本的直方图
>>> import matplotlib.pyplot as plt >>> rvs = rng.rvs(1000) >>> fig = plt.figure() >>> ax = fig.add_subplot(111) >>> x = np.arange(0, n+1) >>> fx = dist.pmf(x) >>> fx = fx / fx.sum() >>> ax.plot(x, fx, 'bo', label='true distribution') >>> ax.vlines(x, 0, fx, lw=2) >>> ax.hist(rvs, bins=np.r_[x, n+1]-0.5, density=True, alpha=0.5, ... color='r', label='samples') >>> ax.set_xlabel('x') >>> ax.set_ylabel('PMF(x)') >>> ax.set_title('Discrete Guide Table Samples') >>> plt.legend() >>> plt.show()
要设置引导表的大小,请使用 guide_factor 关键字参数。这会设置引导表相对于概率向量的尺寸
>>> rng = DiscreteGuideTable(pv, guide_factor=1, random_state=urng)
要计算 \(n=4\) 和 \(p=0.1\) 的二项式分布的 PPF:我们可以按照以下步骤设置引导表
>>> n, p = 4, 0.1 >>> dist = binom(n, p) >>> rng = DiscreteGuideTable(dist, random_state=42) >>> rng.ppf(0.5) 0.0
方法
ppf
(u)给定分布的 PPF。
rvs
([size, random_state])从分布中采样。
set_random_state
([random_state])设置底层均匀随机数生成器。