基于厄米特插值的累积分布函数反演 (HINV)#
必需: CDF
可选: PDF, dPDF
速度
设置: (非常) 慢
采样: (非常) 快
HINV 是数值反演的一种变体，其中使用厄米特插值逼近逆 CDF，即，将区间 [0,1] 分成几个子区间，在每个子区间中，逆 CDF 被用 CDF 和 PDF 在子区间边界处的数值构造的插值多项式逼近。这样就可以通过分割特定子区间来提高精度，而不会重新计算未受影响的子区间。实现了三种类型的样条曲线: 线性、三次和五次插值。对于线性插值，只需要 CDF。三次插值还需要 PDF，而五次插值则需要 PDF 及其导数。
这些样条曲线必须在设置步骤中进行计算。但是，它只适用于具有有界域的分布；对于具有无界域的分布，将截断尾部，使尾部区域的概率与给定的 u 分辨率相比很小。
该方法并不精确，因为它只生成逼近分布的随机变量。然而，"u 方向" (即 |U - CDF(X)|，其中 X 是对应于分位数 U 的近似百分位数，即 X = approx_ppf(U)) 中的最大数值误差可以设置为所需的分辨率 (在机器精度内)。请注意，u 分辨率的非常小的值是可能的，但可能会增加设置步骤的成本。
NumericalInverseHermite 使用厄米特样条曲线逼近连续统计分布的 CDF 的逆。可以通过传递 order 参数来指定厄米特样条曲线的阶数。
如 [1] 中所述，它首先在分布的支持范围内对一组分位数 x 处的分布的 PDF 和 CDF 进行评估。它利用这些结果来拟合一个厄米特样条曲线 H，使得 H(p) == x，其中 p 是与分位数 x 相对应的百分位数数组。因此，该样条曲线在百分位数 p 处以机器精度逼近分布的 CDF 的逆，但通常情况下，样条曲线在百分位数点之间的中点处不会那么准确。
p_mid = (p[:-1] + p[1:])/2
因此，需要根据需要细化分位数网格，以将最大 "u 误差"。
u_error = np.max(np.abs(dist.cdf(H(p_mid)) - p_mid))
降低到指定的容差 u_resolution 以下。当达到所需容差或在下一次细化后网格区间的数量可能超过允许的最大区间数量 (100000) 时，细化停止。
>>> import numpy as np
>>> from scipy.stats.sampling import NumericalInverseHermite
>>> from scipy.stats import norm, genexpon
>>> from scipy.special import ndtr
要创建一个从标准正态分布中采样的生成器，请执行以下操作
>>> class StandardNormal:
... def pdf(self, x):
... return 1/np.sqrt(2*np.pi) * np.exp(-x**2 / 2)
... def cdf(self, x):
... return ndtr(x)
...
>>> dist = StandardNormal()
>>> urng = np.random.default_rng()
>>> rng = NumericalInverseHermite(dist, random_state=urng)
NumericalInverseHermite 具有一个方法，该方法可以逼近分布的 PPF。
>>> rng = NumericalInverseHermite(dist)
>>> p = np.linspace(0.01, 0.99, 99) # percentiles from 1% to 99%
>>> np.allclose(rng.ppf(p), norm.ppf(p))
True
根据分布的随机采样方法的实现，在给定相同的随机状态下，生成的随机变量可能几乎相同。
>>> dist = genexpon(9, 16, 3)
>>> rng = NumericalInverseHermite(dist)
>>> # `seed` ensures identical random streams are used by each `rvs` method
>>> seed = 500072020
>>> rvs1 = dist.rvs(size=100, random_state=np.random.default_rng(seed))
>>> rvs2 = rng.rvs(size=100, random_state=np.random.default_rng(seed))
>>> np.allclose(rvs1, rvs2)
True
为了检查随机变量是否紧密地遵循给定的分布，我们可以查看它的直方图
>>> import matplotlib.pyplot as plt
>>> dist = StandardNormal()
>>> urng = np.random.default_rng()
>>> rng = NumericalInverseHermite(dist, random_state=urng)
>>> rvs = rng.rvs(10000)
>>> x = np.linspace(rvs.min()-0.1, rvs.max()+0.1, 1000)
>>> fx = norm.pdf(x)
>>> plt.plot(x, fx, 'r-', lw=2, label='true distribution')
>>> plt.hist(rvs, bins=20, density=True, alpha=0.8, label='random variates')
>>> plt.xlabel('x')
>>> plt.ylabel('PDF(x)')
>>> plt.title('Numerical Inverse Hermite Samples')
>>> plt.legend()
>>> plt.show()
给定 PDF 相对于变量 (即 x) 的导数，我们可以通过传递 order 参数来使用五次厄米特插值来逼近 PPF
>>> class StandardNormal:
... def pdf(self, x):
... return 1/np.sqrt(2*np.pi) * np.exp(-x**2 / 2)
... def dpdf(self, x):
... return -1/np.sqrt(2*np.pi) * x * np.exp(-x**2 / 2)
... def cdf(self, x):
... return ndtr(x)
...
>>> dist = StandardNormal()
>>> urng = np.random.default_rng()
>>> rng = NumericalInverseHermite(dist, order=5, random_state=urng)
较高的阶数会导致较少的区间数量
>>> rng3 = NumericalInverseHermite(dist, order=3)
>>> rng5 = NumericalInverseHermite(dist, order=5)
>>> rng3.intervals, rng5.intervals
(3000, 522)
可以通过调用 u_error 方法来估计 u 误差。它运行一个小型的蒙特卡罗模拟来估计 u 误差。默认情况下，使用 100,000 个样本。这可以通过传递 sample_size 参数来更改
>>> rng1 = NumericalInverseHermite(dist, u_resolution=1e-10)
>>> rng1.u_error(sample_size=1000000) # uses one million samples
UError(max_error=9.53167544892608e-11, mean_absolute_error=2.2450136432146864e-11)
这将返回一个命名元组，其中包含最大 u 误差和平均绝对 u 误差。
可以通过降低 u 分辨率 (允许的最大 u 误差) 来降低 u 误差
>>> rng2 = NumericalInverseHermite(dist, u_resolution=1e-13)
>>> rng2.u_error(sample_size=1000000)
UError(max_error=9.32027892364129e-14, mean_absolute_error=1.5194172675685075e-14)
请注意，这会导致计算时间的增加，因为设置时间和区间数量的增加。
>>> rng1.intervals
1022
>>> rng2.intervals
5687
>>> from timeit import timeit
>>> f = lambda: NumericalInverseHermite(dist, u_resolution=1e-10)
>>> timeit(f, number=1)
0.017409582000254886 # may vary
>>> f = lambda: NumericalInverseHermite(dist, u_resolution=1e-13)
>>> timeit(f, number=1)
0.08671202100003939 # may vary