基于多项式插值的CDF逆变换（PINV）#
必需：PDF
可选：CDF、众数、中心
速度：
设置：（非常）慢
采样：（非常）快
基于多项式插值的CDF逆变换（PINV）是一种仅需要密度函数即可从分布中采样的逆变换方法。它基于PPF的多项式插值和PDF的高斯-洛巴托积分。它提供了对插值误差和积分误差的控制。其主要目的是以中等至较慢的设置时间为代价，提供非常快速的采样，几乎适用于任何给定的分布。在固定参数情况下，它是已知的最快的逆变换方法。
逆变换方法是采样非均匀随机变量的最简单和最灵活的方法。对于具有CDF \(F\) 的目标分布和从 \(\text{Uniform}(0, 1)\) 采样的均匀随机变量 \(U\)，随机变量 X 通过使用分布的PPF（逆CDF）对均匀随机变量 \(U\) 进行变换生成：
这种方法适用于随机模拟，因为它具有以下优点：
它保留了均匀随机数采样器的结构特性。
将均匀随机变量 \(U\) 一对一地转换为非均匀随机变量 \(X\)。
易于且高效地从截断分布中采样。
不幸的是，PPF很少以封闭形式提供，或者在可用时速度太慢。对于许多分布，CDF也不容易获得。这种方法解决了这两个缺点。用户只需提供PDF，并可选择性地提供一个接近众数的点（称为“中心”）以及最大可接受误差的大小。它结合了自适应和 简单的Gauss-Lobatto积分法来获得CDF，以及牛顿插值法来获得PPF。该方法并不精确，因为它仅生成近似分布的随机变量。然而，最大容许的近似误差可以设置得接近机器精度。使用u-误差的概念来测量和控制误差。它定义为：
其中 \(u \in (0, 1)\) 是我们想要测量误差的分位数，而 \(F^{-1}_a\) 是给定分布的近似PPF。
最大u-误差是计算CDF和PPF时近似误差的准则。算法所能容忍的最大u-误差称为算法的u-分辨率，记作 \(\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。
使用Gauss-Lobatto积分计算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("使用PINV方法抽取的样本。")
>>> plt.legend()
>>> plt.show()
最大容忍误差（即u分辨率）可以通过在初始化时传递``u_resolution``关键字来更改：
>>> rng = NumericalInversePolynomial(dist, u_resolution=1e-12,
... random_state=urng)
这会导致对PPF的更精确近似，并且生成的随机变量更紧密地遵循精确分布。尽管如此，需要注意的是，这会带来昂贵的设置成本。
设置时间主要取决于PDF被评估的次数。对于难以评估的PDF，其成本更高。请注意，我们可以忽略归一化常数以加快PDF的评估速度。在设置过程中，对于较小的``u_resolution``值，PDF评估次数会增加。
>>> 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
>>> # => 需要较少的PDF评估
>>> # => 更快的设置
>>> rng = NumericalInversePolynomial(dist, u_resolution=1e-8,
... random_state=urng)
>>> dist.callbacks
4095 # 可能会有所不同
>>> dist.callbacks = 0 # 重置回调次数
>>> # u_resolution = 10^-10（默认）
>>> # => 需要更多PDF评估
>>> # => 设置较慢
>>> rng = NumericalInversePolynomial(dist, u_resolution=1e-10,
... random_state=urng)
>>> dist.callbacks
11454 # 可能会有所不同
>>> dist.callbacks = 0 # 重置回调次数
>>> # u_resolution = 10^-12
>>> # => 需要大量PDF评估
>>> # => 设置非常慢
>>> rng = NumericalInversePolynomial(dist, u_resolution=1e-12,
... random_state=urng)
13902 # 可能会有所不同
正如我们所见，所需的PDF评估次数非常高，快速的PDF对算法至关重要。尽管如此，这有助于减少为达到误差目标所需的子区间数量，从而节省内存并使采样速度更快。`NumericalInverseHermite`是一种类似的反演方法，它基于Hermite插值反演CDF，并通过u-resolution控制最大容忍误差。但它使用的区间比`NumericalInversePolynomial`多得多：
>>> from scipy.stats.sampling import NumericalInverseHermite
>>> # NumericalInverseHermite接受一个`tol`参数来设置生成器的u-resolution。
>>> rng_hermite = NumericalInverseHermite(norm(), tol=1e-12)
>>> rng_hermite.intervals
3000
>>> rng_poly = NumericalInversePolynomial(norm(), u_resolution=1e-12)
>>> rng_poly.intervals
252
当分布的精确CDF可用时，可以通过调用:func:`~NumericalInversePolynomial.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\) 是 :math:`sqrt{2pi}`（标准正态分布的归一化常数）：
>>> 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仍会被多次调用，从而减慢了计算速度。