Ctrl+K

簡單均勻比率法 (SROU)#

  • 必要條件：PDF，PDF 下的面積（若非 1）

  • 可選條件：眾數，眾數下的 CDF

  • 速度

    • 設定：快速

    • 採樣：慢速

SROU 基於均勻比率法，該方法使用通用不等式來建構（通用）邊界矩形。它適用於 T 凹分佈，其中 T(x) = -1/sqrt(x)。

>>> import numpy as np
>>> from scipy.stats.sampling import SimpleRatioUniforms

假設我們有常態分佈

>>> class StdNorm:
...     def pdf(self, x):
...         return np.exp(-0.5 * x**2)

請注意，PDF 的積分不為 1。我們可以傳遞產生器初始化期間 PDF 下的確切面積，或 PDF 下確切面積的上限。此外，建議傳遞分佈的眾數以加速設定。

>>> urng = np.random.default_rng()
>>> dist = StdNorm()
>>> rng = SimpleRatioUniforms(dist, mode=0,
...                           pdf_area=np.sqrt(2*np.pi),
...                           random_state=urng)

現在，我們可以使用 rvs 方法從分佈中產生樣本

>>> rvs = rng.rvs(10)

如果眾數下的 CDF 可用，則可以設定它以提高 rvs 的效能

>>> from scipy.stats import norm
>>> rng = SimpleRatioUniforms(dist, mode=0,
...                           pdf_area=np.sqrt(2*np.pi),
...                           cdf_at_mode=norm.cdf(0),
...                           random_state=urng)
>>> rvs = rng.rvs(1000)

我們可以透過視覺化其直方圖來檢查樣本是否來自給定的分佈

>>> from scipy.stats.sampling import SimpleRatioUniforms
>>> from scipy.stats import norm
>>> import matplotlib.pyplot as plt
>>> class StdNorm:
...     def pdf(self, x):
...         return np.exp(-0.5 * x**2)
...
>>> urng = np.random.default_rng()
>>> dist = StdNorm()
>>> rng = SimpleRatioUniforms(dist, mode=0,
...                           pdf_area=np.sqrt(2*np.pi),
...                           cdf_at_mode=norm.cdf(0),
...                           random_state=urng)
>>> rvs = rng.rvs(1000)
>>> x = np.linspace(rvs.min()-0.1, rvs.max()+0.1, 1000)
>>> fx = 1/np.sqrt(2*np.pi) * dist.pdf(x)
>>> fig, ax = plt.subplots()
>>> ax.plot(x, fx, 'r-', lw=2, label='true distribution')
>>> ax.hist(rvs, bins=10, density=True, alpha=0.8, label='random variates')
>>> ax.set_xlabel('x')
>>> ax.set_ylabel('PDF(x)')
>>> ax.set_title('Simple Ratio-of-Uniforms Samples')
>>> ax.legend()
>>> plt.show()
" "

該方法的主要優點是設定快速。如果需要重複產生具有不同形狀參數的分佈的小到中等樣本，這可能很有利。在這種情況下，sampling.NumericalInverseHermite 或 sampling.NumericalInversePolynomial 的設定步驟將導致效能不佳。例如，假設我們有興趣為 Gamma 分佈產生 100 個樣本，該分佈具有 1000 個不同的形狀參數，由 np.arange(1.5, 5, 1000) 給出。

>>> import math
>>> class GammaDist:
...     def __init__(self, p):
...         self.p = p
...     def pdf(self, x):
...         return x**(self.p-1) * np.exp(-x)
...
>>> urng = np.random.default_rng()
>>> p = np.arange(1.5, 5, 1000)
>>> res = np.empty((1000, 100))
>>> for i in range(1000):
...     dist = GammaDist(p[i])
...     rng = SimpleRatioUniforms(dist, mode=p[i]-1,
...                               pdf_area=math.gamma(p[i]),
...                               random_state=urng)
...     with np.testing.suppress_warnings() as sup:
...         sup.filter(RuntimeWarning, "invalid value encountered in double_scalars")
...         sup.filter(RuntimeWarning, "overflow encountered in exp")
...         res[i] = rng.rvs(100)

詳情請參閱 [1]、[2] 和 [3]。

參考文獻#

本頁面