离散别名骨灰盒(DAU)

  • 必需:概率向量(PV)或PMF以及有限域

  • 速度:

    • 设置:慢速(与矢量长度成线性)

    • 采样:非常快

DAU样本来自长度为N的任意但有限的概率向量(PV)分布。该算法基于A.J.Walker的一种巧妙的方法,需要一个(至少)N的表。它需要一个随机数,并且每个生成的随机变量只需要一次比较。构造表的设置时间为O(N)。

>>> import numpy as np
>>> from scipy.stats import DiscreteAliasUrn
>>>
>>> pv = [0.18, 0.02, 0.8]
>>> urng = np.random.default_rng()
>>> rng = DiscreteAliasUrn(pv, random_state=urng)
>>> rng.rvs()
0

默认情况下,概率向量从0开始编制索引。但是,这可以通过传递一个 domain 参数。什么时候 domain 与pv结合使用时,它的效果是将分布从 (0, len(pv))(domain[0]domain[0] + len(pv))domain[1] 在这种情况下被忽略。

>>> rng = DiscreteAliasUrn(pv, domain=(10, 13), random_state=urng)
>>> rng.rvs()
12

当只给出PMF而不给出概率向量时,该方法同样有效。在这种情况下,还必须通过将 domain 参数显式设置或通过提供 support 分发对象中的方法:

>>> class Distribution:
...     def __init__(self, c):
...         self.c = c
...     def pmf(self, x):
...         return x**self.c
...     def support(self):
...         return (0, 10)
...
>>> dist = Distribution(2)
>>> rng = DiscreteAliasUrn(dist, random_state=urng)
>>> rng.rvs()
10
>>> import matplotlib.pyplot as plt
>>> from scipy.stats import DiscreteAliasUrn
>>> class Distribution:
...     def __init__(self, c):
...         self.c = c
...     def pmf(self, x):
...         return x**self.c
...     def support(self):
...         return (0, 10)
...
>>> dist = Distribution(2)
>>> urng = np.random.default_rng()
>>> rng = DiscreteAliasUrn(dist, random_state=urng)
>>> rvs = rng.rvs(1000)
>>> fig = plt.figure()
>>> ax = fig.add_subplot(111)
>>> x = np.arange(1, 11)
>>> 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, 11]-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()
../../_images/sampling_dau-1.png

注解

作为 DiscreteAliasUrn 需要带签名的PMF def pmf(self, x: float) -> float ,它首先对PMF进行矢量化,使用 np.vectorize 然后对域中的所有点进行评估。但是,如果PMF已经矢量化,那么只在域中的每个点对其求值并将获得的PV随域一起传递会快得多。例如, pmf SciPy离散分布的方法是矢量化的,并且可以通过执行以下操作来获得PV:

>>> from scipy.stats import binom, DiscreteAliasUrn
>>> dist = binom(10, 0.2)  # distribution object
>>> domain = dist.support()  # the domain of your distribution
>>> x = np.arange(domain[0], domain[1] + 1)
>>> pv = dist.pmf(x)
>>> rng = DiscreteAliasUrn(pv, domain=domain)

此处需要域才能重新定位分发。

设置已用表的大小可能会对性能产生轻微影响,该大小可以通过传递 urn_factor 参数。

>>> # use a table twice the length of PV.
>>> urn_factor = 2
>>> rng = DiscreteAliasUrn(pv, urn_factor=urn_factor, random_state=urng)
>>> rng.rvs()
2

注解

建议将此参数控制在2以下。

请看 12 有关此方法的更多详细信息,请参见。

参考文献

1

联合国参考手册,第5.8.2节,“DAU-(离散)别名-URN方法”,http://statmath.wu.ac.at/software/unuran/doc/unuran.html#DAU

2

A.J.沃克(1977)。一种产生具有一般分布的离散随机变量的有效方法ACM Trans。数学课。软件3,第253-256页。