DiscreteAliasUrn#
- class scipy.stats.sampling.DiscreteAliasUrn(dist, *, domain=None, urn_factor=1, random_state=None)#
离散别名-瓮方法。
此方法用于从具有有限域的单变量离散分布中采样。 它使用大小为 \(N\) 的概率向量或具有有限支持的概率质量函数,从分布中生成随机数。
- 参数:
- dist类数组或对象,可选
分布的概率向量(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]
被忽略。 有关更详细的说明,请参见“注释”和教程。
- urn_factorfloat,可选
相对于概率向量大小的瓮表大小。 它不能小于 1。较大的表会导致更快的生成时间,但需要更昂贵的设置。 默认为 1。
- random_state{None, int,
numpy.random.Generator
, 用于生成均匀随机数流的底层 NumPy 随机数生成器或种子。 如果 random_state 为 None (或 np.random),则使用
numpy.random.RandomState
单例。 如果 random_state 是一个 int,则使用一个新的RandomState
实例,并使用 random_state 作为种子。 如果 random_state 已经是一个Generator
或RandomState
实例,则使用该实例。
方法
rvs
([size, random_state])从分布中采样。
set_random_state
([random_state])设置底层均匀随机数生成器。
备注
当有限概率向量可用或分布的 PMF 可用时,此方法有效。 如果只有 PMF 可用,则还必须给出 PMF 的*有限*支持(域)。 建议首先通过在支持中的每个点评估 PMF 来获得概率向量,然后改用它。
如果给出了概率向量,则它必须是一个非负浮点数的一维数组,没有任何
inf
或nan
值。 此外,必须至少有一个非零条目,否则会引发异常。默认情况下,概率向量从 0 开始索引。但是,可以通过传递
domain
参数来更改此设置。 当domain
与 PV 结合使用时,其效果是将分布从(0, len(pv))
重新定位到(domain[0]
,domain[0] + len(pv))
。 在这种情况下,domain[1]
被忽略。可以增加参数
urn_factor
以加快生成速度,但代价是增加设置时间。 此方法使用一个表进行随机变量生成。urn_factor
控制此表相对于概率向量的大小(或者在 PV 不可用时,控制支持的宽度)。 由于此表是在设置时计算的,因此线性增加此参数会线性增加设置所需的时间。 建议将此参数保持在 2 以下。参考文献
[1]UNU.RAN 参考手册,第 5.8.2 节,“DAU - (Discrete) Alias-Urn method”,http://statmath.wu.ac.at/software/unuran/doc/unuran.html#DAU
[2]A.J. Walker (1977)。 用于生成具有一般分布的离散随机变量的有效方法,ACM Trans。 数学 软件 3,第 253-256 页。
示例
>>> from scipy.stats.sampling import DiscreteAliasUrn >>> import numpy as np
要使用概率向量创建随机数生成器,请使用
>>> pv = [0.1, 0.3, 0.6] >>> urng = np.random.default_rng() >>> rng = DiscreteAliasUrn(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.292, 0.616]) >>> chisquare(freqs, pv).pvalue 0.9993602047563164
由于 p 值非常高,我们无法拒绝观察到的频率与预期频率相同的原假设。 因此,我们可以安全地假设变量是从给定分布生成的。 请注意,这仅给出了算法的正确性,而不是样本的质量。
如果 PV 不可用,也可以传递带有 PMF 方法和有限域的类的实例。
>>> urng = np.random.default_rng() >>> class Binomial: ... def __init__(self, n, p): ... self.n = n ... self.p = p ... def pmf(self, x): ... # note that the pmf doesn't need to be normalized. ... return self.p**x * (1-self.p)**(self.n-x) ... def support(self): ... return (0, self.n) ... >>> n, p = 10, 0.2 >>> dist = Binomial(n, p) >>> rng = DiscreteAliasUrn(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.9999999999999999
要检查样本是否是从正确的分布中提取的,我们可以可视化样本的直方图
>>> 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 Alias Urn Samples') >>> plt.legend() >>> plt.show()
要设置
urn_factor
,请使用>>> rng = DiscreteAliasUrn(pv, urn_factor=2, random_state=urng)
这使用两倍于概率向量大小的表来从分布中生成随机变量。