变换密度抑制(TDR)¶
要求：T形凹面PDF、dPDF
可选：模式、中心
速度：
设置：慢速
采样：快速
TDR是一种利用变换密度的凹性来构造HAT函数并自动压缩的接受/拒绝方法。这样的PDF被称为T形凹面。目前实现了以下转换：
除了PDF之外，它还需要PDF w.r.t的派生 x (即变量)。这些函数必须作为Python对象的方法出现，然后可以传递给生成器以实例化它们的对象。
此方法有三种变体可用：
GW：在施工点之间挤压。
PS：与帽子函数成正比的挤压。(默认)
IA：与变体PS相同，但在挤压下的区域使用“立即接受”的合成方法。
这可以通过传递一个 variant 参数。
使用此方法的示例如下所示：
>>> from scipy.stats 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不需要矢量化。它应该接受并返回一个标量。
直接计算HAT分布的CDF的逆也是可能的 ppf_hat 方法。
>>> 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意味着HAT和挤压函数将分布紧紧包围，生成样本所需的PDF求值更少。密度的预期求值次数受以下限制
(1/squeeze_hat_ratio) - 1每个样品。默认情况下，它保持在0.99以上，但可以通过传递max_squeeze_hat_ratio参数。“HAT_Area”：生成器的HAT下方区域。
“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
可以通过传递一个域参数来截断分布：
>>> 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])
如果未指定域，则 support 的方法。 dist 对象用于确定域：
>>> 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,
... max_intervals=1000, random_state=urng)
>>> rng.squeeze_hat_ratio
0.999364900465214
请注意，我们需要增加 max_intervals 当我们需要更高的 squeeze_hat_ratio 。这是因为需要更多的构造点来更紧密地拟合分布。
让我们看看这对分发的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,
... max_intervals=1000, 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
正如我们所看到的，当我们增加 squeeze_hat_ratio 。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])
还可以通过将 variant 参数：
>>> rng = TransformedDensityRejection(dist, variant='ia')
此方法接受许多其他设置参数。有关独家列表，请参阅文档。有关参数和方法的更多信息，请参阅 Section 5.3.16 of the UNU.RAN user manual 。
请看 1, 2, 和 3 有关此方法的更多详细信息，请参见。
参考文献¶
- 1
UNU.RAN参考手册，第5.3.16节，“TDR-变换密度拒绝”，http://statmath.wu.ac.at/software/unuran/doc/unuran.html#TDR
- 2
霍尔曼，沃尔夫冈。“一种从T凹形分布中抽样的拒收技术。”ACM数学软件学报(TOMS)21.2(1995)：182-193
- 3
W.R.Gilks和P.Wild(1992)。“吉布斯抽样的自适应拒绝抽样”，应用统计学41，第337-348页。