Python SciPy sampling.NumericalInversePolynomial用法及代码示例

本文简要介绍 python 语言中 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)#

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

PINV 是数值反演的一种变体，其中逆 CDF 使用牛顿插值公式进行近似。区间[0,1] 被分成几个子区间。在这些中的每一个中，在该子区间中的某些点 x 在节点 (CDF(x),x) 处构造逆 CDF。如果给出了 PDF，则使用 5 点的自适应Gauss-Lobatto 积分从给定的 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 为 “not too small” 的区域。单峰密度满足这个条件。如果违反此条件，则分布的域可能会被截断。

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

• PDF 必须有界。

• 当分布有重尾时(因为然后逆 CDF 在 0 或 1 处变得非常陡峭)并且请求的 u-resolution 非常小，该算法会出现问题。例如，当请求的u-resolution 小于 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：分布的 CDF。此方法是可选的。如果提供，它将启用 “u-error” 的计算。见 u_error 。必须与 PDF 具有相同的签名。

mode： 浮点数，可选

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

center： 浮点数，可选

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

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

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

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

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

order： 整数，可选

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

u_resolution： 浮点数，可选

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

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]

Derflinger、Gerhard、Wolfgang Hörmann 和 Josef Leydold。 “当仅知道密度时，通过数值反演生成随机变量。” ACM 建模和计算机仿真交易 (TOMACS) 20.4 (2010)：1-25。

[2]

UNU.RAN 参考手册，第 5.3.12 节，“PINV - CDF 的基于多项式插值的 INVersion”，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()

如果在设置期间有确切的 CDF 可用，则可以估计近似 PPF 的 u-error。为此，请在初始化期间传递具有精确分布 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。它运行Monte-Carlo 模拟来估计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-error 可以通过减少u-resolution(允许的最大值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()

属性 ：：

intervals

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