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

获取生成器的电流比率(挤压下方的面积)/(帽子下方的面积)。