scipy.integrate.

qmc_quad#

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

使用拟蒙特卡罗积分法在N维空间中计算积分。

参数:
函数可调用

被积函数。必须接受一个参数 x，这是一个指定在哪些点上计算标量值被积函数的数组，并返回被积函数的值。为了提高效率，函数应向量化以接受形状为 (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``个点上进行评估。详情请参见注释。

qrng : QMCEngine, 可选QMCEngine, 可选

一个用于采样QMC点的QMCEngine实例。QMCEngine必须初始化为与传递给`func`的变量``x1, …, xd``数量相对应的维度``d``。提供的QMCEngine用于生成第一个积分估计。如果`n_estimates`大于一，将从第一个QMCEngine（如果启用了混沌，则启用混沌）生成额外的QMCEngine。如果没有提供QMCEngine，默认将使用`scipy.stats.qmc.Halton`，其维度数量由`a`的长度决定。

日志布尔值，默认：False

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

返回:
结果对象

一个带有属性的结果对象：

积分浮动

积分的估计值。

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
On this page