基于Hermite插值的CDF逆变换（HINV）#
必需：CDF
可选：PDF, dPDF
速度：
设置：（非常）慢
采样：（非常）快
HINV是数值逆变换的一种变体，其中逆CDF通过Hermite插值进行近似，即区间[0,1]被划分为多个子区间，并且在每个子区间内，逆CDF通过在区间边界处使用CDF和PDF的值构建的多项式进行近似。这使得可以通过划分特定区间来提高精度，而无需重新计算不受影响的区间。实现了三种类型的样条：线性、三次和五次插值。对于线性插值，仅需要CDF。三次插值还需要PDF，五次插值则需要PDF及其导数。
这些样条必须在设置步骤中计算。然而，它仅适用于具有有界域的分布；对于具有无界域的分布，尾部会被截断，使得尾部区域的概率相对于给定的u分辨率很小。
该方法不是精确的，因为它仅生成近似分布的随机变量。尽管如此，在“u方向”上的最大数值误差（即``|U - CDF(X)|``，其中``X``是与分位数``U``对应的近似百分位数，即``X = approx_ppf(U)``）可以设置为所需的分辨率（在机器精度范围内）。请注意，非常小的u分辨率值是可能的，但可能会增加设置步骤的成本。
NumericalInverseHermite 使用Hermite样条近似连续统计分布的CDF的逆。可以通过传递`order`参数来指定Hermite样条的阶数。
如[1]_中所述，它首先在分布的支持范围内的一组分位数``x``处评估分布的PDF和CDF。
它使用这些结果来拟合一个Hermite样条曲线 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) # 从1%到99%的百分位数
>>> np.allclose(rng.ppf(p), norm.ppf(p))
True
根据分布的随机采样方法的实现，生成的随机变量在给定相同的随机状态时可能几乎相同。
>>> dist = genexpon(9, 16, 3)
>>> rng = NumericalInverseHermite(dist)
>>> # `seed` 确保每个 `rvs` 方法使用相同的随机流
>>> 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='真实分布')
>>> plt.hist(rvs, bins=20, density=True, alpha=0.8, label='随机变量')
>>> plt.xlabel('x')
>>> plt.ylabel('PDF(x)')
>>> plt.title('数值逆Hermite样本')
>>> plt.legend()
>>> plt.show()
给定PDF相对于变量（即``x``）的导数，我们可以通过传递`order`参数使用五次Hermite插值来近似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-error。它运行一个小型的蒙特卡罗模拟来估计u-error。默认情况下，使用100,000个样本。 used. 可以通过传递 sample_size 参数来更改此设置：
>>> rng1 = NumericalInverseHermite(dist, u_resolution=1e-10)
>>> rng1.u_error(sample_size=1000000) # 使用一百万个样本
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 # 可能会有所不同
>>> f = lambda: NumericalInverseHermite(dist, u_resolution=1e-13)
>>> timeit(f, number=1)
0.08671202100003939 # 可能会有所不同