scipy.stats.NumericalInversePolynomial

class scipy.stats.NumericalInversePolynomial(dist, mode=None, center=None, *, domain=None, order=5, u_resolution=1e-10, max_intervals=10000, keep_cdf=False, random_state=None)

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

PINV是数值反演的一种变体，其中逆CDF是使用牛顿插值公式来近似的。时间间隔 [0,1] 被分成几个子区间。在这些节点中的每一个中，逆CDF都是在节点处构造的 (CDF(x),x) 在某些方面 x 在此子间隔内。如果给定PDF，则从给定的PDF使用带5个点的自适应高斯-洛巴托积分来数值计算CDF。子间隔被拆分，直到达到所请求的精度目标。

这种方法是不精确的，因为它只产生近似分布的随机变量。不过，可以将允许的最大近似误差设置为分辨率(当然，受机器精度的限制)。我们使用u-error |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 以及可选的一个 cdf 方法。

  • pdf ：发行版的PDF格式。PDF的签名应为： def pdf(self, x: float) -> float 即，PDF应该接受Python浮点并返回Python浮点。它不需要集成到1，也就是说，PDF不需要标准化。

  • cdf ：分发的CDF。此方法是可选的。如果提供，它将启用“u-Error”的计算。看见 u_error 。必须具有与PDF相同的签名。

mode浮动，可选

(精确)分发模式。默认值为 None 。

center浮动，可选

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

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

发行版的支持。默认值为 None 。什么时候 None ：

  • 如果一个 support 方法由分发对象提供 dist ，用于设置分布的域。

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

order整型，可选

插值多项式的阶数。有效阶数介于3到17之间。阶数越高，近似的间隔越短。默认值为5。

u_resolution浮动，可选

设置最大容许u误差。u_Resolution的值最多必须为1.e-15和1.e-5。请注意，大多数均匀随机数源的分辨率为2-32=2.3e-10。因此，值1.e-10导致可以称为精确的反演算法。对于大多数仿真来说，稍微大一点的最大误差值也就足够了。默认值为1e-10。

max_intervals整型，可选

设置最大间隔数。必须在范围内 [100,1000000] 。默认值为10000。

keep_cdf布尔值，可选

如果给定了PDF，则使用自适应高斯-洛巴托积分从给定的PDF数值计算CDF。因此，CDF点表被存储，以使PDF的求值次数保持最小。通常，此表在安装完成后会被丢弃。如果 keep_cdf 为真，则保留此表，并且可以通过调用 cdf 方法。默认值为 False 。

random_state ：{无，整型， numpy.random.Generator ，{无，整型，

用于生成均匀随机数流的基础NumPy随机数生成器的NumPy随机数生成器或种子。如果 random_state 为无(或 np.random )、 numpy.random.RandomState 使用的是Singleton。如果 random_state 是一个整型、一个新的 RandomState 实例，其种子设定为 random_state 。如果 random_state 已经是一个 Generator 或 RandomState 实例，则使用该实例。

注意事项

此方法不适用于具有恒定部分的密度(例如 uniform 分布)和分段故障(如果这样的密度被传递给构造器)。建议使用合成法从这些分布中抽样。

参考文献

1

书名/作者Deflinger，Gerhard，Wolfgang Hörmann，和Josef Leydold。“在只知道密度的情况下，通过数值反演产生随机变量。”ACM建模与计算机仿真学报(TOMACS)20.4(2010)：1-25。

2

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

示例

>>> from scipy.stats import NumericalInversePolynomial
>>> from scipy.stats import norm

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

>>> 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-NumericalInversePolynomial-1_00_00.png

如果在设置期间有精确的CDF可用，则可以估计近似PPF的u误差。要执行此操作，请传递一个 dist 初始化期间具有确切分布CDF的对象：

>>> 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误差。默认情况下，使用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分辨率(允许的最大u误差)来减小u误差：

>>> 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 方法：

>>> 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-NumericalInversePolynomial-1_01_00.png
属性
intervals

获取计算中使用的间隔数。

方法：

cdf \(X)

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

ppf \(U)

给定分布的近似PPF。

rvs \([size] )

分发的样本。

set_random_state \([random_state] )

设置基础均匀随机数生成器。

u_error \([sample_size] )

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

上一页

scipy.stats.NumericalInverseHermite.rvs

下一页

scipy.stats.NumericalInversePolynomial.cdf