Ctrl+K
scipy.integrate.

qmc_quad#

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

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

参数:
func可调用对象

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

a, b类数组

一维数组，分别指定每个 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 的长度确定的维数进行初始化。

log布尔值, 默认值: False

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

返回:
result对象

具有以下属性的结果对象

integralfloat

积分的估计值。

standard_error

误差估计。详见“注释”以了解解释。

注释

在 QMC 样本的每个 n_points 个点处的被积函数的值用于生成积分的估计。此估计是从积分的可能估计值总体中提取的，我们获得的值取决于计算积分的特定点。我们执行此过程 n_estimates 次，每次在不同的加扰 QMC 点计算被积函数，有效地从积分估计值总体中抽取 i.i.d. 随机样本。这些积分估计值的样本均值 \(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
在本页