scipy.stats.sampling.

NumericalInversePolynomial#

class scipy.stats.sampling.NumericalInversePolynomial(dist, *, mode=None, center=None, domain=None, order=5, u_resolution=1e-10, random_state=None)#

基于多项式插值的累积分布函数逆变换 (PINV)。

PINV 是数值反演的一种变体，其中逆CDF使用牛顿插值公式进行近似。区间 [0,1] 被分割成几个子区间。在这些子区间中的每一个，逆CDF在节点 (CDF(x),x) 处构造，其中 x 是该子区间内的一些点。如果给出了PDF，则使用5点自适应Gauss-Lobatto积分从给定的PDF中数值计算CDF。子区间被分割，直到达到请求的精度目标。

该方法并不精确，因为它仅生成近似分布的随机变量。尽管如此，最大容忍的近似误差可以设置为分辨率（当然，受限于机器精度）。我们使用 u-误差 |U - CDF(X)| 来测量误差，其中 X 是对应于分位数 U 的近似百分位数，即 X = approx_ppf(U)。我们将最大容忍的 u-误差称为算法的 u-分辨率。

插值多项式的顺序和 u-分辨率都可以选择。请注意，u-分辨率非常小的值是可能的，但会增加设置步骤的成本。

插值多项式必须在设置步骤中计算。然而，它仅适用于具有有界域的分布；对于具有无界域的分布，尾部被截断，使得尾部区域的概率与给定的u分辨率相比很小。

插值多项式的构造仅在PDF是单峰的或当PDF在两个模式之间不消失时有效。

给定的分布有一些限制：

  • 分布的支持（即PDF严格为正的区域）必须是连通的。实际上，这意味着PDF“不太小”的区域必须是连通的。单峰密度满足这一条件。如果违反这一条件，则分布的定义域可能会被截断。

  • 当PDF被数值积分时，给定的PDF必须是连续的，并且应该是平滑的。

  • PDF 必须装订。

  • 当分布具有重尾（此时逆CDF在0或1处变得非常陡峭）且请求的u分辨率非常小时，该算法会出现问题。例如，当请求的u分辨率小于1.e-12时，柯西分布很可能会显示出这个问题。

参数:
dist对象

一个具有 pdf 或 logpdf 方法的类实例，可选地包含 cdf 方法。

  • pdf: 分布的PDF。PDF的签名应为：def pdf(self, x: float) -> float，即PDF应接受一个Python浮点数并返回一个Python浮点数。它不需要积分为1，即PDF不需要归一化。此方法是可选的，但需要指定``pdf``或``logpdf``。如果两者都给出，则使用``logpdf``。

  • logpdf: 分布的PDF的对数。签名与 pdf 相同。同样，可以忽略PDF的归一化常数的对数。此方法是可选的，但需要指定 pdf 或 logpdf 之一。如果两者都给出，则使用 logpdf。

  • cdf: 分布的累积分布函数。此方法是可选的。如果提供，它将启用“u-误差”的计算。参见 u_error。必须与PDF具有相同的签名。

模式float, 可选

(精确) 分布的模式。默认是 None。

中心float, 可选

分布的模式或均值的近似位置。此位置提供了关于PDF主要部分的一些信息，并用于避免数值问题。默认值为 None。

领域长度为2的列表或元组，可选

对分布的支持。默认是 None。当 None 时：

  • 如果分布对象 dist 提供了 support 方法，它将用于设置分布的域。

  • 否则，支持集被假定为 \((-\infty, \infty)\)。

顺序int, 可选

插值多项式的阶数。有效的阶数在3到17之间。更高的阶数会导致近似的区间更少。默认值是5。

u_resolutionfloat, 可选

设置最大容忍的 u-误差。u_resolution 的值必须至少为 1.e-15 且最多为 1.e-5。注意，大多数均匀随机数源的分辨率为 2-32= 2.3e-10。因此，1.e-10 的值会导致一个可以称为精确的反演算法。对于大多数模拟，稍微大一点的值对于最大误差来说也是足够的。默认值为 1e-10。

random_state{None, int,}

用于生成均匀随机数流的底层 NumPy 随机数生成器的 NumPy 随机数生成器或种子。如果 random_state 是 None（或 np.random），则使用 numpy.random.RandomState 单例。如果 random_state 是整数，则使用一个新的 RandomState 实例，并以 random_state 为种子。如果 random_state 已经是 Generator 或 RandomState 实例，则使用该实例。

属性:
intervals

获取计算中使用的区间数量。

方法

cdf(x)

给定分布的近似累积分布函数。

ppf(u)

给定分布的近似PPF。

qrvs([size, d, qmc_engine])

给定随机变量的拟随机变量。

rvs([size, random_state])

来自发行版的示例。

set_random_state([random_state])

设置底层均匀随机数生成器。

u_error([sample_size])

使用蒙特卡罗模拟估计近似的u误差。

参考文献

[1]

Derflinger, Gerhard, Wolfgang Hörmann, 和 Josef Leydold. “仅知密度时的数值反演随机变量生成.” ACM 建模与计算机模拟交易 (TOMACS) 20.4 (2010): 1-25.

[2]

UNU.RAN 参考手册, 第5.3.12节, “PINV - 基于多项式插值的CDF反演”, https://statmath.wu.ac.at/software/unuran/doc/unuran.html#PINV

示例

>>> from scipy.stats.sampling import NumericalInversePolynomial
>>> from scipy.stats import norm
>>> import numpy as np

要创建一个从标准正态分布中采样的生成器，请执行以下操作：

>>> class StandardNormal:
...    def pdf(self, x):
...        return np.exp(-0.5 * x*x)
...
>>> dist = StandardNormal()
>>> urng = np.random.default_rng()
>>> rng = NumericalInversePolynomial(dist, random_state=urng)

一旦生成器被创建，可以通过调用 rvs 方法从分布中抽取样本：

>>> rng.rvs()
-1.5244996276336318

为了检查随机变量是否紧密遵循给定的分布，我们可以查看其直方图：

>>> import matplotlib.pyplot as plt
>>> rvs = rng.rvs(10000)
>>> x = np.linspace(rvs.min()-0.1, rvs.max()+0.1, 1000)
>>> fx = norm.pdf(x)
>>> 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('Numerical Inverse Polynomial Samples')
>>> plt.legend()
>>> plt.show()
../../_images/scipy-stats-sampling-NumericalInversePolynomial-1_00_00.png

如果在设置过程中可以获得精确的CDF，则可以估计近似PPF的u-误差。为此，请在初始化时传递一个包含分布精确CDF的`dist`对象：

>>> from scipy.special import ndtr
>>> class StandardNormal:
...    def pdf(self, x):
...        return np.exp(-0.5 * x*x)
...    def cdf(self, x):
...        return ndtr(x)
...
>>> dist = StandardNormal()
>>> urng = np.random.default_rng()
>>> rng = NumericalInversePolynomial(dist, random_state=urng)

现在，可以通过调用 u_error 方法来估计 u-error。它运行蒙特卡罗模拟来估计 u-error。默认情况下，使用 100000 个样本。要更改此设置，可以将样本数量作为参数传递：

>>> rng.u_error(sample_size=1000000)  # uses one million samples
UError(max_error=8.785994154436594e-11, mean_absolute_error=2.930890027826552e-11)

这将返回一个包含最大 u-误差和平均绝对 u-误差的命名元组。

通过降低 u-resolution（最大允许的 u-error）可以减少 u-error。

>>> urng = np.random.default_rng()
>>> rng = NumericalInversePolynomial(dist, u_resolution=1.e-12, random_state=urng)
>>> rng.u_error(sample_size=1000000)
UError(max_error=9.07496300328603e-13, mean_absolute_error=3.5255644517257716e-13)

请注意，这会增加设置时间。

可以通过调用 ppf 方法来评估近似的 PPF：

>>> rng.ppf(0.975)
1.9599639857012559
>>> norm.ppf(0.975)
1.959963984540054

由于正态分布的PPF（百分点函数）可以作为特殊函数使用，我们也可以检查x误差，即近似PPF与精确PPF之间的差异:

>>> import matplotlib.pyplot as plt
>>> u = np.linspace(0.01, 0.99, 1000)
>>> approxppf = rng.ppf(u)
>>> exactppf = norm.ppf(u)
>>> error = np.abs(exactppf - approxppf)
>>> plt.plot(u, error)
>>> plt.xlabel('u')
>>> plt.ylabel('error')
>>> plt.title('Error between exact and approximated PPF (x-error)')
>>> plt.show()
../../_images/scipy-stats-sampling-NumericalInversePolynomial-1_01_00.png
On this page