Python SciPy sampling.NumericalInverseHermite用法及代碼示例

本文簡要介紹 python 語言中 scipy.stats.sampling.NumericalInverseHermite 的用法。

用法:

class  scipy.stats.sampling.NumericalInverseHermite(dist, *, domain=None, order=3, u_resolution=1e-12, construction_points=None, random_state=None)#

基於 Hermite 插值的 CDF 的 INVersion (HINV)。

HINV 是數值反演的一種變體，其中逆 CDF 使用 Hermite 插值來近似，即區間 [0,1] 被分成幾個區間，並且在每個區間中，逆 CDF 由由區間邊界處的 CDF 和 PDF。這使得可以通過在不受影響的間隔中重新計算而拆分特定間隔來提高準確性。實現了三種類型的樣條：線性、三次和五次插值。對於線性插值，僅需要 CDF。三次插值還需要 PDF 和五次插值 PDF 及其導數。

這些樣條必須在設置步驟中計算。但是，它僅適用於有界域的分布；對於具有無界域的分布，尾部被截斷，使得尾部區域的概率與給定的 u-resolution 相比較小。

該方法並不精確，因為它隻產生近似分布的隨機變量。盡管如此，“u-direction” 中的最大數值誤差(即 |U - CDF(X)| 其中 X 是對應於分位數 U 即 X = approx_ppf(U) 的近似百分位數)可以設置為所需的分辨率(在機器精度範圍內)。請注意，u-resolution 的非常小的值是可能的，但可能會增加設置步驟的成本。

參數 ：：

dist： 對象

具有cdf 和可選的pdf 和dpdf 方法的類的實例。

• cdf：分布的 CDF。 CDF 的簽名預計為：def cdf(self, x: float) -> float。即 CDF 應該接受一個 Python 浮點數並返回一個 Python 浮點數。

• pdf：分發的 PDF。當 order=1 時，此方法是可選的。必須與 PDF 具有相同的簽名。

• dpdf ：PDF w.r.t 變量的導數(即 x )。此方法對於 order=1 或 order=3 是可選的。必須與 CDF 具有相同的簽名。

domain： 長度為 2 的列表或元組，可選

分布的支持。默認為 None 。當 None ：

• 如果一個support方法由分發對象提供距離，用於設置分布的域。

• 否則，假定支持為 \((-\infty, \infty)\) 。

order： 整數，默認值：3

設置 Hermite 插值的順序。有效階數為 1、3 和 5。有效階數為 1、3 和 5。請注意，大於 1 的階數需要分布的密度，大於 3 的階數甚至需要密度的導數。在大量間隔中對大多數分布使用 1 階結果，因此不建議使用。如果 u-direction 中的最大誤差非常小(比如小於 1.e-10)，則建議使用 5 階，因為隻要沒有極點或重尾，它會導致設計點少得多。

u_resolution： 浮點數，默認：1e-12

設置最大容許u-error。請注意，大多數均勻隨機數源的分辨率是 2-32= 2.3e-10。因此，1.e-10 的值導致了一種可以稱為精確的反演算法。對於大多數模擬，最大誤差稍大的值也足夠了。默認值為 1e-12。

construction_points： 數組，可選

設置 Hermite 插值的起始構造點(節點)。由於可能的最大誤差僅在設置中進行估計，因此可能需要設置一些特殊的設計點來計算 Hermite 插值，以確保最大 u-error 不能大於期望值。這些點是密度不可微分或具有局部極值的點。

random_state： {無，整數， numpy.random.Generator ，

numpy.random.RandomState}, optional

NumPy 隨機數生成器或底層NumPy 隨機數生成器的種子，用於生成統一隨機數流。如果random_state是無(或np.random)， 這numpy.random.RandomState使用單例。如果random_state是一個 int，一個新的RandomState使用實例，播種random_state.如果random_state已經是一個Generator或者RandomState實例然後使用該實例。

注意：

NumericalInverseHermite用 Hermite 樣條逼近連續統計分布的 CDF 的逆。 Hermite 樣條的階數可以通過傳遞次序範圍。

如 [1] 中所述，它首先在分布的支持範圍內以分位數x 的網格評估分布的 PDF 和 CDF。它使用結果來擬合 Hermite 樣條 H 使得 H(p) == x ，其中 p 是與分位數 x 相對應的百分位數數組。因此，樣條曲線在百分位 p 處將分布的 CDF 的倒數近似為機器精度，但通常，樣條曲線在百分位點之間的中點處不會那麽準確：

p_mid = (p[:-1] + p[1:])/2

因此分位數的網格根據需要進行細化以減少最大值“u-error”：

u_error = np.max(np.abs(dist.cdf(H(p_mid)) - p_mid))

低於指定的公差 u_resolution。當達到所需的容差或下一次細化後的網格間隔數可能超過允許的最大間隔數(即 100000)時，細化將停止。

參考：

[1]

霍曼、沃爾夫岡和約瑟夫·萊多德。 “通過快速數值反演連續生成隨機變量。” ACM 建模和計算機仿真交易 (TOMACS) 13.4 (2003): 347-362。

[2]

UNU.RAN 參考手冊，第 5.3.5 節，“HINV - 基於 Hermite 插值的 CDF 的 INVersion”，https://statmath.wu.ac.at/software/unuran/doc/unuran.html#HINV

例子：

>>> from scipy.stats.sampling import NumericalInverseHermite
>>> from scipy.stats import norm, genexpon
>>> from scipy.special import ndtr
>>> import numpy as np

要創建從標準正態分布中采樣的生成器，請執行以下操作：

>>> 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()
>>> rng = NumericalInverseHermite(dist)
>>> 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 w.r.t 的導數變量(即x)，我們可以使用五次 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方法。它運行一個小的Monte-Carlo 模擬來估計u-error。默認情況下，使用 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-error 和平均絕對值 u-error 的命名元組。

u-error 可以通過減少u-resolution(允許的最大值u-error)來減少：

>>> 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

由於正態分布的 PPF 可作為特殊函數使用，我們還可以檢查x-error，即近似 PPF 和精確 PPF 之間的差異：

>>> import matplotlib.pyplot as plt
>>> u = np.linspace(0.01, 0.99, 1000)
>>> approxppf = rng.ppf(u)
>>> exactppf = norm.ppf(u)
>>> error = np.abs(exactppf - approxppf)
>>> plt.plot(u, error)
>>> plt.xlabel('u')
>>> plt.ylabel('error')
>>> plt.title('Error between exact and approximated PPF (x-error)')
>>> plt.show()

屬性 ：：

intervals

獲取生成器對象中用於 Hermite 插值的節點數(設計點)。

midpoint_error：