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

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

用法:

class  scipy.stats.sampling.RatioUniforms(pdf, *, umax, vmin, vmax, c=0, random_state=None)#

使用 ratio-of-uniforms 方法從概率密度函數生成隨機樣本。

參數 ：：

pdf： 可調用的

具有簽名 pdf(x) 的函數，與分布的概率密度函數成比例。

umax： 浮點數

u-direction 中邊界矩形的上限。

vmin： 浮點數

v-direction 中邊界矩形的下界。

vmax： 浮點數

v-direction 中邊界矩形的上限。

c： 浮點數，可選。

ratio-of-uniforms 方法的移位參數，請參見注釋。默認值為 0。

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

numpy.random.RandomState}, optional

如果種子是無(或np.random)， 這numpy.random.RandomState使用單例。如果種子是一個 int，一個新的RandomState使用實例，播種種子.如果種子已經是一個Generator或者RandomState實例然後使用該實例。

注意：

給定一個單變量概率密度函數pdf和一個常數c, 定義集合A = {(u, v) : 0 < u <= sqrt(pdf(v/u + c))}.如果(U, V)是一個均勻分布的隨機向量A， 然後V/U + c遵循一個分布pdf.

上述結果(見[1],[2]) 可用於僅使用 PDF 對隨機變量進行采樣，即不需要 CDF 的反轉。典型選擇c是零或模式pdf.套裝A是矩形的子集R = [0, umax] x [vmin, vmax]其中

• umax = sup sqrt(pdf(x))

• vmin = inf (x - c) sqrt(pdf(x))

• vmax = sup (x - c) sqrt(pdf(x))

特別是，這些值是有限的，如果pdf是有界的並且x**2 * pdf(x)是有界的(即次二次尾部)。一個可以生成(U, V)一致地R並返回V/U + c如果(U, V)也在A可以直接驗證。

如果對於任何常數 k > 0，將 pdf 替換為 k * pdf，則該算法不會改變。因此，通過刪除不必要的歸一化因子，使用與概率密度函數成比例的函數通常很方便。

直觀地說，如果A填滿大部分封閉矩形，使得概率很高(U, V)在於A每當它位於R因為否則所需的迭代次數會變得太大。更準確地說，請注意要繪製的預期迭代次數(U, V)均勻分布在R這樣(U, V)也在A由比率給出area(R) / area(A) = 2 * umax * (vmax - vmin) / area(pdf)，其中麵積(pdf)是積分pdf(如果使用概率密度函數，則等於一，但如果使用與密度成比例的函數，則可以取其他值)。等式成立，因為麵積A等於0.5 * area(pdf)(定理 7.1 中[1])。如果采樣在 50000 次迭代後未能生成單個隨機變量(即沒有一次抽簽在A)，引發異常。

如果未正確指定邊界矩形(即，如果它不包含A)，算法從不同於下式給出的分布中采樣pdf.因此，建議執行測試，例如scipy.stats.kstest作為支票。

參考：

[1] (1,2)

L. Devroye，“非均勻隨機變量生成”，Springer-Verlag，1986 年。

[2]

W. Hoermann 和 J. Leydold，“生成廣義逆高斯隨機變量”，統計與計算，24(4)，p。 547-557，2014 年。

[3]

A.J. Kinderman 和 J.F. Monahan，“使用均勻偏差比率的隨機變量的計算機生成”，ACM Transactions on Mathematical Software，3(3)，p。 257-260，1977。

例子：

>>> import numpy as np
>>> from scipy import stats

>>> from scipy.stats.sampling import RatioUniforms
>>> rng = np.random.default_rng()

模擬正態分布的隨機變量。在這種情況下，很容易顯式地計算邊界矩形。為簡單起見，我們刪除了密度的歸一化因子。

>>> f = lambda x: np.exp(-x**2 / 2)
>>> v = np.sqrt(f(np.sqrt(2))) * np.sqrt(2)
>>> umax = np.sqrt(f(0))
>>> gen = RatioUniforms(f, umax=umax, vmin=-v, vmax=v, random_state=rng)
>>> r = gen.rvs(size=2500)

K-S 檢驗確認隨機變量確實是正態分布的(正態性在 5% 顯著性水平上不被拒絕)：

>>> stats.kstest(r, 'norm')[1]
0.250634764150542

index 分布提供了另一個可以明確確定邊界矩形的示例。

>>> gen = RatioUniforms(lambda x: np.exp(-x), umax=1, vmin=0,
...                     vmax=2*np.exp(-1), random_state=rng)
>>> r = gen.rvs(1000)
>>> stats.kstest(r, 'expon')[1]
0.21121052054580314