scipy.integrate.

qmc_quad#

scipy.integrate.qmc_quad(func, a, b, *, n_estimates=8, n_points=1024, qrng=None, log=False)[source]#

使用拟蒙特卡罗积分法计算 N 维积分。

参数:
funccallable

被积函数。必须接受一个参数 x，一个指定要评估标量值被积函数的点的数组，并返回被积函数的值。为了提高效率，该函数应该被矢量化以接受形状为 (d, n_points) 的数组，其中 d 是变量的数量（即函数域的维数），n_points 是积分点的数量，并返回形状为 (n_points,) 的数组，每个积分点的被积函数值。

a, barray-like

一维数组，分别指定每个 d 变量的积分下限和上限。

n_estimates, n_pointsint, 可选

n_estimates（默认：8）个统计独立的 QMC 样本，每个样本有 n_points（默认：1024）个点，将由 qrng 生成。被积函数 func 将被评估的点的总数为 n_points * n_estimates。有关详细信息，请参见注释。

qrngQMCEngine, 可选

一个 QMCEngine 实例，从中采样 QMC 点。QMCEngine 必须初始化为与传递给 func 的变量数量 x1, ..., xd 相对应的维度数量 d。提供的 QMCEngine 用于生成第一个积分估计。如果 n_estimates 大于 1，则从第一个 QMCEngine 衍生出额外的 QMCEngine（如果它是选项，则启用随机化）。如果未提供 QMCEngine，则将使用默认的 scipy.stats.qmc.Halton 进行初始化，其维度数量由 a 的长度决定。

logboolean, default: False

当设置为 True 时，func 返回被积函数的对数，结果对象包含积分的对数。

返回值:
resultobject

一个具有以下属性的结果对象

integralfloat

积分估计值。

standard_error

误差估计值。有关解释，请参见注释。

注释

QMC 样本的每个 n_points 个点的被积函数值用于生成积分估计值。该估计值来自积分估计值的可能估计值总体，我们获得的值取决于评估积分的特定点。我们重复此过程 n_estimates 次，每次都评估不同随机化 QMC 点的被积函数，实际上是从积分估计值总体中抽取独立同分布的随机样本。这些积分估计值的样本均值 \(m\) 是积分真实值的无偏估计量，这些估计值的均值标准误差 \(s\) 可用于使用具有 n_estimates - 1 个自由度的 t 分布生成置信区间。也许与直觉相反，在保持函数评估点总数 n_points * n_estimates 不变的情况下增加 n_points 往往会减少实际误差，而增加 n_estimates 往往会减少误差估计值。

示例

QMC 积分法在计算更高维度的积分方面特别有用。一个示例被积函数是多元正态分布的概率密度函数。

>>> import numpy as np
>>> from scipy import stats
>>> dim = 8
>>> mean = np.zeros(dim)
>>> cov = np.eye(dim)
>>> def func(x):
...     # `multivariate_normal` expects the _last_ axis to correspond with
...     # the dimensionality of the space, so `x` must be transposed
...     return stats.multivariate_normal.pdf(x.T, mean, cov)

为了计算单位超立方体上的积分

>>> from scipy.integrate import qmc_quad
>>> a = np.zeros(dim)
>>> b = np.ones(dim)
>>> rng = np.random.default_rng()
>>> qrng = stats.qmc.Halton(d=dim, seed=rng)
>>> n_estimates = 8
>>> res = qmc_quad(func, a, b, n_estimates=n_estimates, qrng=qrng)
>>> res.integral, res.standard_error
(0.00018429555666024108, 1.0389431116001344e-07)

积分的双边 99% 置信区间可估计为

>>> t = stats.t(df=n_estimates-1, loc=res.integral,
...             scale=res.standard_error)
>>> t.interval(0.99)
(0.0001839319802536469, 0.00018465913306683527)

实际上，scipy.stats.multivariate_normal 报告的值在这个范围内。

>>> stats.multivariate_normal.cdf(b, mean, cov, lower_limit=a)
0.00018430867675187443
本页内容