scipy.stats.sampling.

NumericalInverseHermite#

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

基于埃尔米特插值的 CDF 的逆 (HINV)

HINV 是数值反演的一种变体，其中逆 CDF 使用埃尔米特插值来近似，也就是说，[0,1] 区间被分成几个区间，并且在每个区间内逆 CDF 由使用区间边界处的 CDF 和 PDF 的值构建的多项式逼近。这使得通过在不受影响的区间内进行重新计算的情况下，通过分割特定区间来提高准确性成为可能。已实现三种类型的样条：线性、三次和五次插值。对于线性插值，仅需要 CDF。三次插值还需要 PDF，五次插值需要 PDF 及其导出。

这些样条曲线必须在设置步骤中进行计算。但是，它仅适用于范围有界的分布；对于范围无界的分布，会将尾部切去，以便与给定 u 分辨率相比，尾部区域的概率很小。

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

参数:

distobject

具有 cdf 方法（以及可选的 pdf 和 dpdf 方法）的类的实例。

cdf: 分布的 CDF。CDF 的签名预计如下：def cdf(self, x: float) -> float。即 CDF 应接受 Python 浮点数并返回 Python 浮点数。
pdf: 分布的 PDF。当 order=1 时，此方法是可选的。必须具有与 PDF 相同的签名。
dpdf: PDF 关于变量的导数（即 x）。当 order=1 或 order=3 时，此方法是可选的。必须具有与 CDF 相同的签名。

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

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

如果分布对象 dist 提供 support 方法，则使用它设置分布的域。
否则，假定支持为 \((-\infty, \infty)\)。

orderint，默认值：3

设置 Hermite 插值的阶次。有效阶次为 1、3 和 5。有效阶次为 1、3 和 5。请注意，大于 1 的阶次需要分布的密度，而大于 3 的阶次甚至需要密度的导数。对于大多数分布，使用阶次 1 会产生大量间隔，因此不推荐使用。如果 u 方向的最大误差很小（比如小于 1.e-10），建议使用阶次 5，因为它会产生明显更少的设置点，只要没有极点或重尾。

u_resolutionfloat，默认值：1e-12

设置允许的最大 u 误差。请注意，大多数均匀随机数源的分辨率为 2-32= 2.3e-10。因此，值 1.e-10 会产生一个可以称之为精密的反演算法。对于大多数仿真来说，将最大误差略微调高一些也足够了。默认值为 1e-12。

construction_pointsarray_like，可选

设置 Hermite 插法的起始设置点（节点）。由于在设置中只估计了可能的最大误差，因此可能需要为计算 Hermite 插值设置一些特殊的设计点，以确保最大 u 误差不超过所需值。这样的点是不连续或具有局部极值的点。

random_state{None，int，numpy.random.Generator,

numpy.random.RandomState}, 可选

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

注释

NumericalInverseHermite 使用埃尔米特样条逼近连续统计分布的CDF的逆函数。可以通过传递order参数来指定埃尔米特样条的阶数。

如[1]所述，它首先在分布的支持内的量分位数网格x处评估分布的PDF和CDF。它使用结果拟合埃尔米特样条H，使得H(p) == x，其中p是与量分位数x相对应的百分位数数组。因此，样条在百分位数p处以机器精度逼近分布的CDF的逆函数，但通常情况下，样条在百分位数点之间的中点处不会如此精确

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

因此，根据需要细化量分位数的网格以减少最大“u-误差”

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

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

参考文献

[1]

Hörmann, Wolfgang，和Josef Leydold。“快速数值反演法产生连续随机变量”。ACM建模和计算机模拟杂志（TOMACS）13.4（2003）：347-362。

[2]

UNU.RAN参考手册，第5.3.5节，“HINV - 基于埃尔米特插值的CDF反演”，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()

../../_images/scipy-stats-sampling-NumericalInverseHermite-1_00_00.png

给定PDF相对于变量（即x）的导数，我们可以通过传递order参数来使用五次埃尔米特插值来逼近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 误差。它运行小型蒙特卡罗模拟来估算 u 误差。默认情况下，使用 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 误差和平均绝对 u 误差的命名元组。

可以通过减小 u 分辨率（允许的最大 u 误差）来减小 u 误差

>>> 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 误差，即近似 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-NumericalInverseHermite-1_01_00.png

属性：

intervals: 获取生成器对象中用于 Hermite 插值中的节点数（设计点）。
midpoint_error

方法

`ppf`(u)	给定分布的近似 PPF。
`qrvs`([size, d, qmc_engine])	给定 RV 的准随机变量。
`rvs`([size, random_state])	从分布中进行采样。
`set_random_state`([random_state])	设置底层均匀随机数生成器。
`u_error`([sample_size])	使用蒙特卡罗模拟估算近似的 u 误差。