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

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

用法:

class  scipy.stats.sampling.TransformedDensityRejection(dist, *, mode=None, center=None, domain=None, c=-0.5, construction_points=30, use_dars=True, max_squeeze_hat_ratio=0.99, random_state=None)#

變換密度抑製 (TDR) 方法。

TDR是一種接受/拒絕方法，它使用變換密度的凹度來構造帽子函數並自動擠壓。與專門針對該分布的算法相比，大多數通用算法都非常慢。快速的算法設置緩慢並且需要大表。這種通用方法的目的是提供一種不太慢且隻需要較短設置的算法。此方法可應用於具有T-concave 密度函數的單變量和單峰連續分布。有關詳細信息，請參閱 [1] 和 [2]。

參數 ：：

dist： 對象

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

• pdf：分發的 PDF。 PDF 的簽名應為：def pdf(self, x: float) -> float。即 PDF 應該接受 Python 浮點數並返回一個 Python 浮點數。它不需要集成到 1，即 PDF 不需要被規範化。

• dpdf：PDF w.r.t x 的導數(即變量)。必須與 PDF 具有相同的簽名。

mode： 浮點數，可選

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

center： 浮點數，可選

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

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

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

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

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

c： {-0.5, 0.}，可選

為轉換函數 T 設置參數 c 。默認值為 -0.5。為了構造帽子函數，PDF 的變換必須是凹的。這樣的 PDF 稱為T-concave。目前支持以下轉換：

\[\begin{split}c = 0.: T(x) &= \log(x)\\ c = -0.5: T(x) &= \frac{1}{\sqrt{x}} \text{ (Default)}\end{split}\]

construction_points： int 或 數組，可選

如果是整數，則定義構造點的數量。如果它是類似數組的，則數組的元素將用作構造點。默認值為 30。

use_dars： 布爾型，可選

如果為 True，則在設置中使用“去隨機化自適應拒絕采樣”(DARS)。有關 DARS 算法的詳細信息，請參見 [1]。默認為真。

max_squeeze_hat_ratio： 浮點數，可選

設置比率的上限(擠壓下方的區域)/(帽子下方的區域)。它必須是介於 0 和 1 之間的數字。默認值為 0.99。

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] (1,2)

UNU.RAN 參考手冊，第 5.3.16 節，“TDR - 變換密度抑製”，http://statmath.wu.ac.at/software/unuran/doc/unuran.html#TDR

[2]

霍曼，沃爾夫岡。 “一種從 T-concave 分布中抽樣的拒絕技術。” ACM 數學軟件交易 (TOMS) 21.2 (1995): 182-193

[3]

W.R. Gilks 和 P. Wild (1992)。吉布斯采樣的自適應拒絕采樣，應用統計 41，第 337-348 頁。

例子：

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

假設我們有一個密度：

\[\begin{split}f(x) = \begin{cases} 1 - x^2, & -1 \leq x \leq 1 \\ 0, & \text{otherwise} \end{cases}\end{split}\]

這個密度函數的導數是：

\[\begin{split}\frac{df(x)}{dx} = \begin{cases} -2x, & -1 \leq x \leq 1 \\ 0, & \text{otherwise} \end{cases}\end{split}\]

請注意，PDF 沒有集成到 1。由於這是一種基於拒絕的方法，我們不需要標準化的 PDF。要初始化生成器，我們可以使用：

>>> urng = np.random.default_rng()
>>> class MyDist:
...     def pdf(self, x):
...         return 1-x*x
...     def dpdf(self, x):
...         return -2*x
...
>>> dist = MyDist()
>>> rng = TransformedDensityRejection(dist, domain=(-1, 1),
...                                   random_state=urng)

域對於截斷分布非常有用，但為了避免每次都將其傳遞給構造函數，可以通過在分布對象 (dist) 中提供支持方法來設置默認域：

>>> class MyDist:
...     def pdf(self, x):
...         return 1-x*x
...     def dpdf(self, x):
...         return -2*x
...     def support(self):
...         return (-1, 1)
...
>>> dist = MyDist()
>>> rng = TransformedDensityRejection(dist, random_state=urng)

現在，我們可以使用 rvs 方法從分布中生成樣本：

>>> rvs = rng.rvs(1000)

我們可以通過可視化其直方圖來檢查樣本是否來自給定分布：

>>> import matplotlib.pyplot as plt
>>> x = np.linspace(-1, 1, 1000)
>>> fx = 3/4 * dist.pdf(x)  # 3/4 is the normalizing constant
>>> 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('Transformed Density Rejection Samples')
>>> plt.legend()
>>> plt.show()

屬性 ：：

hat_area

獲取生成器帽子下方的區域。

squeeze_area

獲取生成器擠壓下方的區域。

squeeze_hat_ratio

獲取生成器的電流比率(擠壓下方的麵積)/(帽子下方的麵積)。