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

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

用法:

class  scipy.stats.sampling.NumericalInversePolynomial(dist, *, mode=None, center=None, domain=None, order=5, u_resolution=1e-10, random_state=None)#

基於多項式插值的 CDF (PINV) 的 INVersion。

PINV 是數值反演的一種變體，其中逆 CDF 使用牛頓插值公式進行近似。區間[0,1] 被分成幾個子區間。在這些中的每一個中，在該子區間中的某些點 x 在節點 (CDF(x),x) 處構造逆 CDF。如果給出了 PDF，則使用 5 點的自適應Gauss-Lobatto 積分從給定的 PDF 數值計算 CDF。子區間被分割，直到達到要求的準確度目標。

該方法並不精確，因為它隻產生近似分布的隨機變量。盡管如此，最大容許近似誤差可以設置為分辨率(但當然，受機器精度的限製)。我們使用 u-error |U - CDF(X)| 來測量誤差，其中 X 是對應於分位數 U 即 X = approx_ppf(U) 的近似百分位數。我們將算法的最大容許u-error稱為u-resolution。

可以選擇插值多項式的階數和u-resolution。請注意，u-resolution 的非常小的值是可能的，但會增加設置步驟的成本。

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

插值多項式的構造僅在 PDF 為單峰或 PDF 在兩種模式之間不消失時才有效。

給定的分布有一些限製：

• 必須連接分布的支持(即 PDF 嚴格為正的區域)。實際上，這意味著必須連接 PDF 為 “not too small” 的區域。單峰密度滿足這個條件。如果違反此條件，則分布的域可能會被截斷。

• 當 PDF 被數值積分時，給定的 PDF 必須是連續的並且應該是平滑的。

• PDF 必須有界。

• 當分布有重尾時(因為然後逆 CDF 在 0 或 1 處變得非常陡峭)並且請求的 u-resolution 非常小，該算法會出現問題。例如，當請求的u-resolution 小於 1.e-12 時，柯西分布可能會顯示此問題。

參數 ：：

dist： 對象

具有 pdf 或 logpdf 方法的類的實例，也可以選擇 cdf 方法。

• pdf：發行版的 PDF。 PDF 的簽名預計為： def pdf(self, x: float) -> float ，即 PDF 應接受 Python 浮點數並返回 Python 浮點數。它不需要積分為1，即PDF不需要標準化。此方法是可選的，但需要指定 pdf 或 logpdf。如果兩者都給出，則使用logpdf。

• logpdf：發行版 PDF 的日誌。簽名與 pdf 相同。類似地，PDF的歸一化常數的對數可以被忽略。此方法是可選的，但需要指定 pdf 或 logpdf。如果兩者都給出，則使用logpdf。

• cdf：分布的 CDF。此方法是可選的。如果提供，它將啟用 “u-error” 的計算。見 u_error 。必須與 PDF 具有相同的簽名。

mode： 浮點數，可選

(精確)分布模式。默認為 None 。

center： 浮點數，可選

模式的近似位置或分布的平均值。此位置提供有關 PDF 主要部分的一些信息，用於避免數字問題。默認為 None 。

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

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

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

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

order： 整數，可選

插值多項式的階數。有效階數介於 3 到 17 之間。階數越高，近似值的間隔就越小。默認值為 5。

u_resolution： 浮點數，可選

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

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實例然後使用該實例。

參考：

[1]

Derflinger、Gerhard、Wolfgang Hörmann 和 Josef Leydold。 “當僅知道密度時，通過數值反演生成隨機變量。” ACM 建模和計算機仿真交易 (TOMACS) 20.4 (2010)：1-25。

[2]

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

例子：

>>> from scipy.stats.sampling import NumericalInversePolynomial
>>> from scipy.stats import norm
>>> import numpy as np

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

>>> 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()
-1.5244996276336318

要檢查隨機變量是否緊密遵循給定分布，我們可以查看它的直方圖：

>>> import matplotlib.pyplot as plt
>>> 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 Polynomial Samples')
>>> plt.legend()
>>> plt.show()

如果在設置期間有確切的 CDF 可用，則可以估計近似 PPF 的 u-error。為此，請在初始化期間傳遞具有精確分布 CDF 的 dist 對象：

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

現在，可以通過調用 u_error 方法來估計u-error。它運行Monte-Carlo 模擬來估計u-error。默認情況下，使用 100000 個樣本。要更改這一點，您可以將樣本數作為參數傳遞：

>>> rng.u_error(sample_size=1000000)  # uses one million samples
UError(max_error=8.785994154436594e-11, mean_absolute_error=2.930890027826552e-11)

這將返回一個包含最大 u-error 和平均絕對值 u-error 的命名元組。

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

>>> urng = np.random.default_rng()
>>> rng = NumericalInversePolynomial(dist, u_resolution=1.e-12, random_state=urng)
>>> rng.u_error(sample_size=1000000)
UError(max_error=9.07496300328603e-13, mean_absolute_error=3.5255644517257716e-13)

請注意，這是以增加設置時間為代價的。

可以通過調用 ppf 方法來評估近似的 PPF：

>>> rng.ppf(0.975)
1.9599639857012559
>>> norm.ppf(0.975)
1.959963984540054

由於正態分布的 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

獲取計算中使用的間隔數。