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 凹密度函数的单变量和单峰连续分布。有关更多详细信息,请参见 [1] 和 [2]。
- 参数:
- dist对象
具有
pdf
和dpdf
方法的类的实例。pdf
:分布的 PDF。PDF 的签名应为:def pdf(self, x: float) -> float
。即 PDF 应接受 Python float 并返回 Python float。它不需要积分到 1,即 PDF 不需要标准化。dpdf
:PDF 相对于 x(即变量)的导数。必须与 PDF 具有相同的签名。
- modefloat,可选
分布的(精确)模式。默认值为
None
。- centerfloat,可选
分布的近似模式或均值位置。此位置提供有关 PDF 主要部分的一些信息,并用于避免数值问题。默认值为
None
。- domain长度为 2 的列表或元组,可选
分布的支持。默认值为
None
。当None
时如果分布对象 dist 提供
support
方法,则使用该方法设置分布的域。否则,支持被假定为 \((-\infty, \infty)\)。
- c{-0.5, 0.},可选
为变换函数
T
设置参数c
。默认值为 -0.5。为了构造帽子函数,PDF 的变换必须是凹的。这种 PDF 称为 T-凹函数。目前支持以下变换\[\begin{split}c = 0.: T(x) &= \log(x)\\ c = -0.5: T(x) &= \frac{1}{\sqrt{x}} \text{ (Default)}\end{split}\]- construction_pointsint 或 array_like,可选
如果为整数,则定义构造点的数量。如果为类数组,则数组的元素用作构造点。默认值为 30。
- use_darsbool,可选
如果为 True,则在设置中使用“去随机化自适应拒绝抽样” (DARS)。有关 DARS 算法的详细信息,请参见 [1]。默认值为 True。
- max_squeeze_hat_ratiofloat,可选
设置比率 (挤压下方的面积) / (帽子下方的面积) 的上限。它必须是介于 0 和 1 之间的数字。默认值为 0.99。
- random_state{None, int,
numpy.random.Generator
, 用于生成均匀随机数流的底层 NumPy 随机数生成器或种子。如果 random_state 为 None(或 np.random),则使用
numpy.random.RandomState
单例。如果 random_state 是一个 int,则使用一个新的RandomState
实例,并用 random_state 进行播种。如果 random_state 已经是一个Generator
或RandomState
实例,则使用该实例。
- 属性:
hat_area
获取生成器的帽子下方的面积。
squeeze_area
获取生成器的挤压下方的面积。
squeeze_hat_ratio
获取生成器的当前比率(挤压下方的面积)/(帽子下方的面积)。
方法
ppf_hat
(u)在 u 处评估帽子分布的 CDF 的反函数。
rvs
([size, random_state])从分布中采样。
set_random_state
([random_state])设置底层均匀随机数生成器。
参考文献
[1] (1,2)UNU.RAN 参考手册,第 5.3.16 节,“TDR - 变换密度拒绝”,http://statmath.wu.ac.at/software/unuran/doc/unuran.html#TDR
[2]Hörmann, Wolfgang。“一种从 T 凹分布中采样的拒绝技术。” ACM Transactions on Mathematical Software (TOMS) 21.2 (1995): 182-193
[3]W.R. Gilks 和 P. Wild (1992)。Gibbs 采样的自适应拒绝抽样,应用统计 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) 中提供 support 方法来设置默认域
>>> 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()