基於 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 的反函數。Hermite 樣條曲線的階數可以通過傳遞 order 參數來指定。
如 [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) # 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 參數來使用五次 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 誤差可以通過呼叫 u_error 方法來估計。它運行一個小的蒙地卡羅模擬來估計 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)
這會傳回一個 namedtuple，其中包含最大 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