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

本文简要介绍 python 语言中 scipy.stats.sampling.NumericalInverseHermite 的用法。

用法:

class  scipy.stats.sampling.NumericalInverseHermite(dist, *, domain=None, order=3, u_resolution=1e-12, construction_points=None, random_state=None)#

基于 Hermite 插值的 CDF 的 INVersion (HINV)。

HINV 是数值反演的一种变体，其中逆 CDF 使用 Hermite 插值来近似，即区间 [0,1] 被分成几个区间，并且在每个区间中，逆 CDF 由由区间边界处的 CDF 和 PDF。这使得可以通过在不受影响的间隔中重新计算而拆分特定间隔来提高准确性。实现了三种类型的样条：线性、三次和五次插值。对于线性插值，仅需要 CDF。三次插值还需要 PDF 和五次插值 PDF 及其导数。

这些样条必须在设置步骤中计算。但是，它仅适用于有界域的分布；对于具有无界域的分布，尾部被截断，使得尾部区域的概率与给定的 u-resolution 相比较小。

该方法并不精确，因为它只产生近似分布的随机变量。尽管如此，“u-direction” 中的最大数值误差(即 |U - CDF(X)| 其中 X 是对应于分位数 U 即 X = approx_ppf(U) 的近似百分位数)可以设置为所需的分辨率(在机器精度范围内)。请注意，u-resolution 的非常小的值是可能的，但可能会增加设置步骤的成本。

参数 ：：

dist： 对象

具有cdf 和可选的pdf 和dpdf 方法的类的实例。

• cdf：分布的 CDF。 CDF 的签名预计为：def cdf(self, x: float) -> float。即 CDF 应该接受一个 Python 浮点数并返回一个 Python 浮点数。

• pdf：分发的 PDF。当 order=1 时，此方法是可选的。必须与 PDF 具有相同的签名。

• dpdf ：PDF w.r.t 变量的导数(即 x )。此方法对于 order=1 或 order=3 是可选的。必须与 CDF 具有相同的签名。

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

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

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

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

order： 整数，默认值：3

设置 Hermite 插值的顺序。有效阶数为 1、3 和 5。有效阶数为 1、3 和 5。请注意，大于 1 的阶数需要分布的密度，大于 3 的阶数甚至需要密度的导数。在大量间隔中对大多数分布使用 1 阶结果，因此不建议使用。如果 u-direction 中的最大误差非常小(比如小于 1.e-10)，则建议使用 5 阶，因为只要没有极点或重尾，它会导致设计点少得多。

u_resolution： 浮点数，默认：1e-12

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

construction_points： 数组，可选

设置 Hermite 插值的起始构造点(节点)。由于可能的最大误差仅在设置中进行估计，因此可能需要设置一些特殊的设计点来计算 Hermite 插值，以确保最大 u-error 不能大于期望值。这些点是密度不可微分或具有局部极值的点。

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实例然后使用该实例。

注意：

NumericalInverseHermite用 Hermite 样条逼近连续统计分布的 CDF 的逆。 Hermite 样条的阶数可以通过传递次序范围。

如 [1] 中所述，它首先在分布的支持范围内以分位数x 的网格评估分布的 PDF 和 CDF。它使用结果来拟合 Hermite 样条 H 使得 H(p) == x ，其中 p 是与分位数 x 相对应的百分位数数组。因此，样条曲线在百分位 p 处将分布的 CDF 的倒数近似为机器精度，但通常，样条曲线在百分位点之间的中点处不会那么准确：

p_mid = (p[:-1] + p[1:])/2

因此分位数的网格根据需要进行细化以减少最大值“u-error”：

u_error = np.max(np.abs(dist.cdf(H(p_mid)) - p_mid))

低于指定的公差 u_resolution。当达到所需的容差或下一次细化后的网格间隔数可能超过允许的最大间隔数(即 100000)时，细化将停止。

参考：

[1]

霍曼、沃尔夫冈和约瑟夫·莱多德。 “通过快速数值反演连续生成随机变量。” ACM 建模和计算机仿真交易 (TOMACS) 13.4 (2003): 347-362。

[2]

UNU.RAN 参考手册，第 5.3.5 节，“HINV - 基于 Hermite 插值的 CDF 的 INVersion”，https://statmath.wu.ac.at/software/unuran/doc/unuran.html#HINV

例子：

>>> from scipy.stats.sampling import NumericalInverseHermite
>>> from scipy.stats import norm, genexpon
>>> from scipy.special import ndtr
>>> import numpy as np

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

>>> class StandardNormal:
...     def pdf(self, x):
...        return 1/np.sqrt(2*np.pi) * np.exp(-x**2 / 2)
...     def cdf(self, x):
...        return ndtr(x)
...
>>> dist = StandardNormal()
>>> urng = np.random.default_rng()
>>> rng = NumericalInverseHermite(dist, random_state=urng)

NumericalInverseHermite 有一种近似分布 PPF 的方法。

>>> rng = NumericalInverseHermite(dist)
>>> p = np.linspace(0.01, 0.99, 99) # percentiles from 1% to 99%
>>> np.allclose(rng.ppf(p), norm.ppf(p))
True

根据分布的随机抽样方法的实现，在相同的随机状态下，生成的随机变量可能几乎相同。

>>> dist = genexpon(9, 16, 3)
>>> rng = NumericalInverseHermite(dist)
>>> # `seed` ensures identical random streams are used by each `rvs` method
>>> seed = 500072020
>>> rvs1 = dist.rvs(size=100, random_state=np.random.default_rng(seed))
>>> rvs2 = rng.rvs(size=100, random_state=np.random.default_rng(seed))
>>> np.allclose(rvs1, rvs2)
True

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

>>> import matplotlib.pyplot as plt
>>> dist = StandardNormal()
>>> rng = NumericalInverseHermite(dist)
>>> 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 Hermite Samples')
>>> plt.legend()
>>> plt.show()

给定 PDF w.r.t 的导数变量(即x)，我们可以使用五次 Hermite 插值来近似 PPF，通过传递次序范围：

>>> class StandardNormal:
...     def pdf(self, x):
...        return 1/np.sqrt(2*np.pi) * np.exp(-x**2 / 2)
...     def dpdf(self, x):
...        return -1/np.sqrt(2*np.pi) * x * np.exp(-x**2 / 2)
...     def cdf(self, x):
...        return ndtr(x)
...
>>> dist = StandardNormal()
>>> urng = np.random.default_rng()
>>> rng = NumericalInverseHermite(dist, order=5, random_state=urng)

阶数越高，间隔数越少：

>>> rng3 = NumericalInverseHermite(dist, order=3)
>>> rng5 = NumericalInverseHermite(dist, order=5)
>>> rng3.intervals, rng5.intervals
(3000, 522)

u-error 可以通过调用u_error方法。它运行一个小的Monte-Carlo 模拟来估计u-error。默认情况下，使用 100,000 个样本。这可以通过传递sample_size争论：

>>> rng1 = NumericalInverseHermite(dist, u_resolution=1e-10)
>>> rng1.u_error(sample_size=1000000)  # uses one million samples
UError(max_error=9.53167544892608e-11, mean_absolute_error=2.2450136432146864e-11)

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

u-error 可以通过减少u-resolution(允许的最大值u-error)来减少：

>>> rng2 = NumericalInverseHermite(dist, u_resolution=1e-13)
>>> rng2.u_error(sample_size=1000000)
UError(max_error=9.32027892364129e-14, mean_absolute_error=1.5194172675685075e-14)

请注意，这是以增加设置时间和间隔数量为代价的。

>>> rng1.intervals
1022
>>> rng2.intervals
5687
>>> from timeit import timeit
>>> f = lambda: NumericalInverseHermite(dist, u_resolution=1e-10)
>>> timeit(f, number=1)
0.017409582000254886  # may vary
>>> f = lambda: NumericalInverseHermite(dist, u_resolution=1e-13)
>>> timeit(f, number=1)
0.08671202100003939  # may vary

由于正态分布的 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

获取生成器对象中用于 Hermite 插值的节点数(设计点)。

midpoint_error：