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

  • 要求：PDF

  • 可选：CDF、MODE、CENTER

  • 速度：

    • 设置：(非常)缓慢

    • 采样：(非常)快

基于多项式插值的CDF反演(PINV)是一种只需要密度函数对分布进行采样的反演方法。它是基于多项式插值的PPF和高斯-洛巴托积分的PDF。它提供对插补误差和积分误差的控制。它的主要目的是提供非常快的采样，这种采样对于任何给定的分布几乎是相同的，但代价是中等到慢的设置时间。对于固定参数情况，它是已知最快的反演方法。

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

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

由于该方法的优点，适用于随机模拟。其中一些最具吸引力的是：

  • 它保留了均匀随机数采样器的结构特性。

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

  • 来自截断分布的简单而有效的抽样。

不幸的是，PPF很少以封闭的形式提供，或者当有的时候太慢了。对于许多发行版，CDF也不容易获得。这种方法解决了这两个缺点。用户只需提供PDF，并且可选地提供模式附近的点(称为“中心”)以及最大可接受误差的大小。它使用自适应求积和简单的Gauss-Lobatto求积相结合的方法来获得CDF，并使用牛顿插值来获得PPF。这种方法是不精确的，因为它只产生近似分布的随机变量。然而，可以将最大容许近似误差设置为接近机器精度。利用u误差的概念对误差进行测量和控制。它被定义为：

\[\Epsilon_{u}(U)=|u-F\Left(F^{-1}_{a}(U)\Right)|\]

哪里 \(u \in (0, 1)\) 是我们要测量误差的分位数， \(F^{{-1}}_a\) 是给定分布的近似PPF。

在数值计算CDF和PPF时，最大u误差是近似误差的判据。算法允许的最大u-误差称为算法的u-分辨率，由 \(\epsilon_{{u}}\) ：

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

u误差的主要优点是，如果CDF可用，则可以容易地计算u误差。我们指的是 1 以进行更详细的讨论。

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

给定的分发有一些限制：

  • 分布的支持(即，PDF严格为正的区域)必须是连接的。实际上，这意味着PDF“不太小”的区域必须连接起来。单峰密度满足这一条件。如果违反此条件，则分布的域可能会被截断。

  • 当对PDF进行数值积分时，给定的PDF必须是连续的，并且应该是平滑的。

  • PDF必须是有界的。

  • 当分布有重尾时，该算法有问题(因为这样，逆CDF在0或1处变得非常陡峭)，并且所要求的u分辨率非常小。例如，当所请求的u分辨率小于1.e-12时，柯西分布很可能显示该问题。

警告

此方法不适用于具有恒定部分的密度(例如 uniform distribution) and segmentation faults if such a density is passed to the constructor. It is recommended to use the composition method 从这样的分布中取样。

该算法在设置期间执行以下四个步骤：

  • 计算分布的终点：如果给定有限支撑，则跳过此步骤。否则，尾巴的末端被切断，使得它们下面的面积小于或等于 \(0.05\epsilon_{{u}}\) 。

  • 将区域划分为子区间以计算CDF和PPF。

  • CDF的计算采用高斯-洛巴托求积，使得积分误差最大 \(0.05I_{{0}}\epsilon_{{u}}\) 哪里 \(I_{{0}}\) 大约是PDF下的总面积。

  • 用牛顿插值公式计算具有最大插补误差的PPF \(0.9\epsilon_{{u}}\) 。

要将生成器初始化为来自标准正态分布的样本，请执行以下操作：

>>> from scipy.stats 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 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()
../../_images/sampling_pinv-1.png

最大允许误差(即u分辨率)可以通过传递 u_resolution 初始化过程中的关键字：

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

这可以更精确地近似PPF，并且生成的RV更接近精确分布。不过，请注意，这是以昂贵的设置为代价的。

设置时间主要取决于评估PDF的次数。对于难以评估的PDF来说，它的成本更高。请注意，我们可以忽略归一化常量以加快PDF的计算速度。对于以下较小的值，PDF评估在设置过程中会增加 u_resolution 。

>>> from scipy.stats 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
>>> 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
>>> 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

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

>>> from scipy.stats 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
3340
>>> rng_poly = NumericalInversePolynomial(norm(), u_resolution=1e-12)
>>> rng_poly.intervals
252

当分布的CDF精确可用时，可以通过调用 u_error 方法：

>>> 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 runs a monte carlo simulation with a given number of samples to estimate the u-error. In the above example, 100,000 samples are used by the simulation to approximate the u-error. It returns the maximum u-error (max_error) and the mean absolute u-error (mean_absolute_error) in a UError namedtuple. As we can see, max_error is below the default 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，请传递 keep_cdf=True 初始化期间：

>>> rng = NumericalInversePolynomial(dist, keep_cdf=True, random_state=urng)

现在，可以通过调用 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的调用，但由于微调步骤使用简单的Gauss-Lobatto求积，PDF被调用多次，从而降低了计算速度。

参考文献

1

书名/作者Deflinger，Gerhard，Wolfgang Hörmann，和Josef Leydold。“在只知道密度的情况下，通过数值反演产生随机变量。”ACM建模与计算机仿真学报(TOMACS)20.4(2010)：1-25。

上一页

离散别名骨灰盒(DAU)

下一页

朴素制服比(NROU)法