基于多项式插值的CDF反演 (PINV)#
必需: PDF
可选: CDF、众数、中心
速度
设置: (非常) 慢
采样: (非常) 快
基于多项式插值的CDF反演 (PINV) 是一种仅需要密度函数来从分布中采样的反演方法。它基于 PPF 的多项式插值和 PDF 的高斯-勒让德积分。它提供对插值误差和积分误差的控制。它的主要目的是提供非常快的采样,对于任何给定分布而言,速度几乎相同,但代价是设置时间适中到缓慢。对于固定参数的情况,它是已知最快的反演方法。
反演方法是采样非均匀随机变量最简单、最灵活的方法。对于具有 CDF \(F\) 的目标分布和从 \(\text{Uniform}(0, 1)\) 中采样的均匀随机变量 \(U\),随机变量 X 是通过使用该分布的 PPF(逆 CDF)转换均匀随机变量 \(U\) 来生成的。
由于其优势,这种方法适合于随机模拟。一些最具吸引力的优势包括
它保留了均匀随机数生成器的结构属性。
将均匀随机变量 \(U\) 一对一地转换为非均匀随机变量 \(X\)。
从截断分布中轻松高效地采样。
不幸的是,PPF 很少以闭合形式给出,或者在可用时速度太慢。对于许多分布,CDF 也难以获得。此方法解决了这两个缺点。用户只需提供 PDF 以及可选的靠近众数的点(称为“中心”) 以及最大可接受误差的大小。它结合了自适应和简单的高斯-勒让德求积法来获得 CDF,并使用牛顿插值法来获得 PPF。该方法不精确,因为它只生成近似分布的随机变量。然而,最大容许近似误差可以设置得接近机器精度。u 误差的概念用于衡量和控制误差。它定义为
其中 \(u \in (0, 1)\) 是我们要测量误差的 quantile,\(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。
使用高斯-勒让德求积法计算 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 被调用了几次,从而减慢了计算速度。