scipy.stats.sampling.

NumericalInversePolynomial#

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

基于 CDF 逆反 (PINV) 的多项式插值。

PINV 是数值反演的一种变体，其中使用牛顿插值公式近似逆 CDF。间隔 [0,1] 被分成多个子间隔。在每一个子间隔中，逆 CDF 在节点 (CDF(x),x) 处构建，其中 x 是该子间隔中的某个点。如果给定 PDF，则使用自适应高斯-洛巴托积分和 5 个点从给定的 PDF 中以数字方式计算 CDF。直到达到请求的精度目标，子间隔才会被分割。

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

插值多项式的阶数和 u-resolution 都可以被选择。请注意，u-resolution 的值可以非常小，但是这会增加设置步骤的成本。

插值多项式必须在设置步骤中被计算。但是，它只适用于具有有界域的分布；对于具有无界域的分布，尾部会被切掉，以便尾部区域的概率与给定的 u-resolution 相比很小。

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

对于给定的分布有一些限制

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

  • 如果 PDF 以数字方式积分，那么给定的 PDF 必须是连续的，并且应该是平滑的。

  • PDF 必须是有界的。

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

参数:
distobject

具有 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：分布的 CDF。该方法是可选的。如果提供，则可计算“u 错误”。请参阅 u_error。必须与 PDF 具有相同的签名。

mode可选的 float

分布的（精确）模式。默认值为 None。

center可选的 float

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

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

分布的支持。默认值为 None。当 None

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

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

order可选的 int

插值多项式的阶数。有效的阶数介于 3 到 17 之间。更高的阶数导致近似值的间隔更小。默认值为 5。

u_resolution浮点数，可选

设置最大的容忍 u 误差。u_resolution 的值必须至少为 1.e-15，最多为 1.e-5。请注意，大多数均匀随机数源的分辨率为 2-32= 2.3e-10。因此，值 1.e-10 会导致反转算法被称作精确算法。对于大多数模拟，最大误差的略大一些的值就足够了。默认值为 1e-10。

random_state无、整数、numpy.random.Generator,

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

引用

[1]

Derflinger, Gerhard, Wolfgang Hörmann, and Josef Leydold. “Random variate generation by numerical inversion when only the density is known.” ACM Transactions on Modeling and Computer Simulation (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-error 和平均绝对 u-error。

可以通过减小 u 分辨率（允许的最大 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-error，即近似 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
属性:
间隔

获得计算中使用到的间隔数。

方法

cdf(x)

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

ppf(u)

给定分布的近似 PPF。

qrvs([size, d, qmc_engine])

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

rvs([size, random_state])

从分布中采样。

set_random_state([random_state])

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

u_error([sample_size])

使用蒙特卡罗模拟估计近似的 u-error。
本页内容