变换密度舍选法 (TDR)#
必需:T-凹 PDF,dPDF
可选:mode,center
速度
设置:慢
采样:快
TDR 是一种接受/拒绝方法,它使用变换密度的凹度来自动构建帽子函数和挤压函数。这种 PDF 称为 T-凹 PDF。目前实现了以下转换
除了 PDF 之外,它还需要 PDF 相对于 x
(即变量)的导数。这些函数必须作为 Python 对象的方法存在,然后可以将其传递给生成器以实例化其对象。实现的变体使用与帽子函数成比例的挤压函数 ([1])。
下面显示了使用此方法的一个示例
>>> import numpy as np
>>> from scipy.stats.sampling import TransformedDensityRejection
>>> from scipy.stats import norm
>>>
>>> class StandardNormal:
... def pdf(self, x):
... # note that the normalization constant is not required
... return np.exp(-0.5 * x*x)
... def dpdf(self, x):
... return -x * np.exp(-0.5 * x*x)
...
>>> dist = StandardNormal()
>>>
>>> urng = np.random.default_rng()
>>> rng = TransformedDensityRejection(dist, random_state=urng)
>>> rng.rvs()
-1.526829048388144
在上面的示例中,我们使用了 TDR 方法从标准正态分布中进行采样。请注意,在计算 PDF 时,我们可以删除归一化常数。这通常有助于加快采样阶段。另请注意,PDF 不需要向量化。它应该接受并返回一个标量。
也可以使用 ppf_hat
方法直接评估帽子分布的 CDF 的逆函数。
>>> rng.ppf_hat(0.5)
-0.00018050266342362759
>>> norm.ppf(0.5)
0.0
>>> u = np.linspace(0, 1, num=10)
>>> rng.ppf_hat(u)
array([ -inf, -1.22227372, -0.7656556 , -0.43135731, -0.14002921,
0.13966423, 0.43096141, 0.76517113, 1.22185606, inf])
>>> norm.ppf(u)
array([ -inf, -1.22064035, -0.76470967, -0.4307273 , -0.1397103 ,
0.1397103 , 0.4307273 , 0.76470967, 1.22064035, inf])
除了 PPF 方法之外,还可以访问其他属性,以查看生成器如何很好地拟合给定的分布。这些是
“squeeze_hat_ratio”:生成器的(挤压函数下面的面积)/(帽子函数下面的面积)。它是一个介于 0 和 1 之间的数字。越接近 1,表示帽子函数和挤压函数越紧密地包围分布,生成样本所需的 PDF 评估就越少。密度评估的预期数量受到每个样本
(1/squeeze_hat_ratio) - 1
的限制。默认情况下,它保持在 0.99 以上,但可以通过传递max_squeeze_hat_ratio
参数来更改。“hat_area”:生成器的帽子函数下面的面积。
“squeeze_area”:生成器的挤压函数下面的面积。
>>> rng.squeeze_hat_ratio 0.9947024204884917 >>> rng.hat_area 2.510253139791547 >>> rng.squeeze_area 2.4969548741894876 >>> rng.squeeze_hat_ratio == rng.squeeze_area / rng.hat_area True
可以通过传递 domain 参数来截断分布
>>> urng = np.random.default_rng()
>>> rng = TransformedDensityRejection(dist, domain=[0, 1], random_state=urng)
>>> rng.rvs(10)
array([0.05452512, 0.97251362, 0.49955877, 0.82789729, 0.33048885,
0.55558548, 0.23168323, 0.13423275, 0.73176575, 0.35739799])
如果未指定域,则使用 dist
对象的 support
方法来确定域
>>> class StandardNormal:
... def pdf(self, x):
... return np.exp(-0.5 * x*x)
... def dpdf(self, x):
... return -x * np.exp(-0.5 * x*x)
... def support(self):
... return -np.inf, np.inf
...
>>> dist = StandardNormal()
>>>
>>> urng = np.random.default_rng()
>>> rng = TransformedDensityRejection(dist, random_state=urng)
>>> rng.rvs(10)
array([-1.52682905, 2.06206883, 0.15205036, 1.11587367, -0.30775562,
0.29879802, -0.61858268, -1.01049115, 0.78853694, -0.23060766])
如果 dist
对象未提供 support
方法,则假定域为 (-np.inf, np.inf)
。
要增加 squeeze_hat_ratio
,请传递 max_squeeze_hat_ratio
>>> dist = StandardNormal()
>>> rng = TransformedDensityRejection(dist, max_squeeze_hat_ratio=0.999,
... random_state=urng)
>>> rng.squeeze_hat_ratio
0.999364900465214
让我们看看这如何影响分布的 PDF 方法的回调
>>> from copy import copy
>>> class StandardNormal:
... def __init__(self):
... self.callbacks = 0
... def pdf(self, x):
... self.callbacks += 1
... return np.exp(-0.5 * x*x)
... def dpdf(self, x):
... return -x * np.exp(-0.5 * x*x)
...
>>> dist1 = StandardNormal()
>>> urng1 = np.random.default_rng()
>>> urng2 = copy(urng1)
>>> rng1 = TransformedDensityRejection(dist1, random_state=urng1)
>>> dist1.callbacks # evaluations during setup
139
>>> dist1.callbacks = 0 # don't consider evaluations during setup
>>> rvs = rng1.rvs(100000)
>>> dist1.callbacks # evaluations during sampling
527
>>> dist2 = StandardNormal()
>>> # use the same stream of uniform random numbers
>>> rng2 = TransformedDensityRejection(dist2, max_squeeze_hat_ratio=0.999,
... random_state=urng2)
>>> dist2.callbacks # evaluations during setup
467
>>> dist2.callbacks = 0 # don't consider evaluations during setup
>>> rvs = rng2.rvs(100000)
>>> dist2.callbacks # evaluations during sampling
84 # may vary
我们可以看到,当我们增加 squeeze_hat_ratio
时,采样期间所需的 PDF 评估要少得多。PPF-hat 函数也更准确
>>> abs(norm.ppf(0.975) - rng1.ppf_hat(0.975))
0.0027054565421578136
>>> abs(norm.ppf(0.975) - rng2.ppf_hat(0.975))
0.00047824084476300044
但是,请注意,这是以设置期间增加的 PDF 评估为代价的。
对于模式不接近 0 的密度,建议通过传递 mode
或 center
参数来设置分布的模式或中心。后者是模式或分布均值的近似位置。此位置提供有关 PDF 主要部分的一些信息,并用于避免数值问题。
>>> # mode = 0 for our distribution
>>> # if exact mode is not available, pass 'center' parameter instead
>>> rng = TransformedDensityRejection(dist, mode=0.)
默认情况下,该方法使用 30 个构造点来构造帽子函数。可以通过传递 construction_points
参数来更改此设置,该参数可以是构造点数组,也可以是表示要使用的构造点数量的整数。
>>> rng = TransformedDensityRejection(dist,
... construction_points=[-5, 0, 5])
此方法接受许多其他设置参数。有关独占列表,请参阅文档。有关参数和该方法的更多信息,请参见 UNU.RAN 用户手册的第 5.3.16 节。