基于多项式插值的CDF反演 (PINV)#

  • 必需: PDF

  • 可选: CDF、众数、中心

  • 速度

    • 设置: (非常) 慢

    • 采样: (非常) 快

基于多项式插值的CDF反演 (PINV) 是一种仅需要密度函数来从分布中采样的反演方法。它基于 PPF 的多项式插值和 PDF 的高斯-勒让德积分。它提供对插值误差和积分误差的控制。它的主要目的是提供非常快的采样,对于任何给定分布而言,速度几乎相同,但代价是设置时间适中到缓慢。对于固定参数的情况,它是已知最快的反演方法。

反演方法是采样非均匀随机变量最简单、最灵活的方法。对于具有 CDF \(F\) 的目标分布和从 \(\text{Uniform}(0, 1)\) 中采样的均匀随机变量 \(U\),随机变量 X 是通过使用该分布的 PPF(逆 CDF)转换均匀随机变量 \(U\) 来生成的。

\[X = F^{-1}(U)\]

由于其优势,这种方法适合于随机模拟。一些最具吸引力的优势包括

  • 它保留了均匀随机数生成器的结构属性。

  • 将均匀随机变量 \(U\) 一对一地转换为非均匀随机变量 \(X\)

  • 从截断分布中轻松高效地采样。

不幸的是,PPF 很少以闭合形式给出,或者在可用时速度太慢。对于许多分布,CDF 也难以获得。此方法解决了这两个缺点。用户只需提供 PDF 以及可选的靠近众数的点(称为“中心”) 以及最大可接受误差的大小。它结合了自适应和简单的高斯-勒让德求积法来获得 CDF,并使用牛顿插值法来获得 PPF。该方法不精确,因为它只生成近似分布的随机变量。然而,最大容许近似误差可以设置得接近机器精度。u 误差的概念用于衡量和控制误差。它定义为

\[\epsilon_{u}(u) = | u - F\left(F^{-1}_{a}(u)\right) |\]

其中 \(u \in (0, 1)\) 是我们要测量误差的 quantile,\(F^{-1}_a\) 是给定分布的近似 PPF。

最大 u 误差是在数值计算 CDF 和 PPF 时用于近似误差的标准。算法的最大容许 u 误差称为算法的 u 分辨率,用 \(\epsilon_{u}\) 表示

\[\sup_{u \in (0,1)} | u - F\left(F^{-1}_{a}(u)\right) | \le \epsilon_{u}\]

u 误差的主要优势是,如果 CDF 可用,则可以轻松计算。关于更详细的讨论,请参考 [1]

此外,该方法仅适用于有界分布。对于无限尾部,尾部末端会被截断,使得它们下的面积小于或等于 \(0.05\epsilon_{u}\)

给定分布存在一些限制

  • 分布的支持范围(即 PDF 严格为正的区域)必须是连通的。在实践中,这意味着 PDF “不太小” 的区域必须是连通的。单峰密度满足此条件。如果此条件不满足,则该分布的域可能会被截断。

  • 当数值积分 PDF 时,给定的 PDF 必须是连续的,并且应是平滑的。

  • PDF 必须是有界的。

  • 当分布具有重尾时(因为此时逆 CDF 在 0 或 1 处变得非常陡峭),并且请求的 u 分辨率非常小时,该算法会出现问题。例如,当请求的 u 分辨率小于 1.e-12 时,柯西分布可能会出现此问题。

在设置过程中,算法会执行以下四个步骤

  • 计算分布的端点: 如果给定了有限的支持范围,则跳过此步骤。否则,尾部末端会被截断,使得它们下的面积小于或等于 \(0.05\epsilon_{u}\)

  • 该域被划分为子区间,用于计算 CDF 和 PPF。

  • 使用高斯-勒让德求积法计算 CDF,使得积分误差最多为 \(0.05I_{0}\epsilon_{u}\),其中 \(I_{0}\) 近似等于 PDF 下的总面积。

  • 使用牛顿插值公式计算 PPF,最大插值误差为 \(0.9\epsilon_{u}\)

要将生成器初始化为从标准正态分布中采样,请执行以下操作

>>> import numpy as np
>>> from scipy.stats.sampling import NumericalInversePolynomial
>>> class StandardNormal:
...     def pdf(self, x):
...         return np.exp(-0.5 * x*x)
...
>>> dist = StandardNormal()
>>> urng = np.random.default_rng()
>>> rng = NumericalInversePolynomial(dist, random_state=urng)

生成器已设置,我们可以开始从我们的分布中采样

>>> rng.rvs((5, 3))
array([[-1.52449963,  1.31933688,  2.05884468],
       [ 0.48883931,  0.15207903, -0.02150773],
       [ 1.11486463,  1.95449597, -0.30724928],
       [ 0.9854643 ,  0.29867424,  0.7560304 ],
       [-0.61776203,  0.16033378, -1.00933003]])

我们可以查看随机变量的直方图,以检查它们与我们的分布的匹配程度

>>> import matplotlib.pyplot as plt
>>> from scipy.stats import norm
>>> from scipy.stats.sampling import NumericalInversePolynomial
>>> class StandardNormal:
...     def pdf(self, x):
...         return np.exp(-0.5 * x*x)
...
>>> dist = StandardNormal()
>>> urng = np.random.default_rng()
>>> rng = NumericalInversePolynomial(dist, random_state=urng)
>>> rvs = rng.rvs(10000)
>>> x = np.linspace(rvs.min()-0.1, rvs.max()+0.1, num=10000)
>>> fx = norm.pdf(x)
>>> plt.plot(x, fx, "r-", label="pdf")
>>> plt.hist(rvs, bins=50, density=True, alpha=0.8, label="rvs")
>>> plt.xlabel("x")
>>> plt.ylabel("PDF(x)")
>>> plt.title("Samples drawn using PINV method.")
>>> plt.legend()
>>> plt.show()
" "

可以通过在初始化期间传递 u_resolution 关键字来更改最大容许误差(即 u 分辨率)

>>> rng = NumericalInversePolynomial(dist, u_resolution=1e-12,
...                                  random_state=urng)

这将导致对 PPF 进行更精确的近似,并且生成的 RV 将更紧密地遵循精确分布。尽管如此,请注意,这将带来昂贵的设置成本。

设置时间主要取决于 PDF 被评估的次数。对于难以评估的 PDF,它更昂贵。请注意,我们可以忽略归一化常数以加快 PDF 的评估速度。在设置过程中,PDF 评估次数会随着 u_resolution 值的减小而增加。

>>> from scipy.stats.sampling import NumericalInversePolynomial
>>> class StandardNormal:
...     def __init__(self):
...         self.callbacks = 0
...     def pdf(self, x):
...         self.callbacks += 1
...         return np.exp(-0.5 * x*x)
...
>>> dist = StandardNormal()
>>> urng = np.random.default_rng()
>>> # u_resolution = 10^-8
>>> # => less PDF evaluations required
>>> # => faster setup
>>> rng = NumericalInversePolynomial(dist, u_resolution=1e-8,
...                                  random_state=urng)
>>> dist.callbacks
4095        # may vary
>>> dist.callbacks = 0  # reset the number of callbacks
>>> # u_resolution = 10^-10 (default)
>>> # => more PDF evaluations required
>>> # => slow setup
>>> rng = NumericalInversePolynomial(dist, u_resolution=1e-10,
...                                  random_state=urng)
>>> dist.callbacks
11454       # may vary
>>> dist.callbacks = 0  # reset the number of callbacks
>>> # u_resolution = 10^-12
>>> # => lots of PDF evaluations required
>>> # => very slow setup
>>> rng = NumericalInversePolynomial(dist, u_resolution=1e-12,
...                                  random_state=urng)
13902     # may vary

正如我们所见,所需的 PDF 评估次数非常高,并且快速的 PDF 对于该算法至关重要。不过,这有助于减少实现误差目标所需的子区间数量,从而节省内存并使采样速度更快。 NumericalInverseHermite 是一种类似的反演方法,它基于 Hermite 插值来反演 CDF,并通过 u 分辨率提供对最大容许误差的控制。但与 NumericalInversePolynomial 相比,它使用了更多区间。

>>> from scipy.stats.sampling import NumericalInverseHermite
>>> # NumericalInverseHermite accepts a `tol` parameter to set the
>>> # u-resolution of the generator.
>>> rng_hermite = NumericalInverseHermite(norm(), tol=1e-12)
>>> rng_hermite.intervals
3000
>>> rng_poly = NumericalInversePolynomial(norm(), u_resolution=1e-12)
>>> rng_poly.intervals
252

当分布的精确 CDF 可用时,可以通过调用 u_error 方法来估计该算法实现的 u 误差

>>> from scipy.special import ndtr
>>> class StandardNormal:
...     def pdf(self, x):
...         return np.exp(-0.5 * x*x)
...     def cdf(self, x):
...         return ndtr(x)
...
>>> dist = StandardNormal()
>>> urng = np.random.default_rng()
>>> rng = NumericalInversePolynomial(dist, random_state=urng)
>>> rng.u_error(sample_size=100_000)
UError(max_error=8.785949745515609e-11, mean_absolute_error=2.9307548109436816e-11)

u_error 运行一个具有给定样本数的蒙特卡罗模拟,以估计 u 误差。在上面的示例中,模拟使用了 100,000 个样本来近似 u 误差。它在 UError 命名的元组中返回最大 u 误差 (max_error) 和平均绝对 u 误差 (mean_absolute_error)。正如我们所见,max_error 低于默认的 u_resolution (1e-10)。

在初始化生成器后,也可以评估给定分布的 PPF

>>> rng.ppf(0.975)
1.959963985701268
>>> norm.ppf(0.975)
1.959963984540054

例如,我们可以使用它来检查最大和平均绝对 u 误差

>>> u = np.linspace(0.001, 0.999, num=1_000_000)
>>> u_errors = np.abs(u - dist.cdf(rng.ppf(u)))
>>> u_errors.max()
8.78600525666684e-11
>>> u_errors.mean()
2.9321444940323206e-11

生成器提供的近似 PPF 方法比分布的精确 PPF 评估速度快得多。

在设置期间,会存储一个 CDF 点表,该表可用于在生成器创建后近似 CDF

>>> rng.cdf(1.959963984540054)
0.9750000000042454
>>> norm.cdf(1.959963984540054)
0.975

我们可以使用它来检查计算 CDF 时的积分误差是否超过 \(0.05I_{0}\epsilon_{u}\)。这里 \(I_0\)\(\sqrt{2\pi}\)(标准正态分布的归一化常数)

>>> x = np.linspace(-10, 10, num=100_000)
>>> x_error = np.abs(dist.cdf(x) - rng.cdf(x))
>>> x_error.max()
4.506062190046123e-12
>>> I0 = np.sqrt(2*np.pi)
>>> max_integration_error = 0.05 * I0 * 1e-10
>>> x_error.max() <= max_integration_error
True

在设置期间计算的 CDF 表用于评估 CDF,并且只需要一些进一步的微调。这减少了对 PDF 的调用,但由于微调步骤使用了简单的高斯-勒让德求积法,因此 PDF 被调用了几次,从而减慢了计算速度。

参考文献#