变换密度拒绝法 (TDR)#

  • 必需: T-凹 PDF, dPDF

  • 可选: 众数, 中心

  • 速度:

    • 设置: 慢

    • 采样: 快

TDR 是一种接受/拒绝方法，利用变换密度的凹性自动构建帽函数和挤压。这种PDF称为T-凹。目前实现了以下变换：

\[\begin{split}c = 0.: T(x) &= \log(x)\\ c = -0.5: T(x) &= \frac{1}{\sqrt{x}} \text{ (默认)}\end{split}\]

除了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):
...         # 注意不需要归一化常数
...         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

可以通过传递一个域参数来截断分布：

>>> 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  # 设置期间的评估次数
139
>>> dist1.callbacks = 0  # 不考虑设置期间的评估次数
>>> rvs = rng1.rvs(100000)
>>> dist1.callbacks  # 采样期间的评估次数
527
>>> dist2 = StandardNormal()
>>> # 使用相同的均匀随机数流
>>> rng2 = TransformedDensityRejection(dist2, max_squeeze_hat_ratio=0.999,
...                                    random_state=urng2)
>>> dist2.callbacks  # 设置期间的评估次数
467
>>> dist2.callbacks = 0  # 设置阶段不考虑评估
>>> rvs = rng2.rvs(100000)
>>> dist2.callbacks  # 采样期间的评估
84      # 可能会有所不同

正如我们所见，当我们增加 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主要部分的一些信息，并用于避免数值问题。

>>> # 我们分布的模式 = 0
>>> # 如果无法获得精确模式，请传递 'center' 参数
>>> rng = TransformedDensityRejection(dist, mode=0.)

默认情况下，该方法使用30个构造点来构建hat函数。这可以通过传递 construction_points 参数来更改，该参数可以是构造点数组或表示要使用的构造点数量的整数。

>>> rng = TransformedDensityRejection(dist,
...                                   construction_points=[-5, 0, 5])

此方法接受许多其他设置参数。请参阅文档以获取独占列表。有关参数和方法的更多信息，请参见 UNU.RAN用户手册的第5.3.16节。

请参阅 [1] 和 [2] 以获取有关此方法的更多详细信息。

参考文献#

On this page