Python SciPy integrate.qmc_quad用法及代码示例

本文简要介绍 python 语言中 scipy.integrate.qmc_quad 的用法。

用法:

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

使用 Quasi-Monte 卡洛求积计算 N-dimensions 中的积分。

参数 ：：

func： 可调用的

被积函数。必须接受单个参数x，一个数组，指定计算 scalar-valued 被积函数的点，并返回被积函数的值。为了提高效率，函数应该被向量化以接受形状数组(d, n_points)，其中d是变量的数量(即函数域的维数)并且n_points是正交点数，返回形状数组(n_points,)，每个正交点的被积函数。

a, b： 类数组

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

n_estimates, n_points： 整数，可选

n_estimates(默认值：8)统计上独立的 QMC 样本，每个样本n_points(默认：1024)点，将由以下方式生成qrng。被积函数的总点数函数将被评估为n_points * n_estimates。详细信息请参见注释。

qrng： QMCEngine ，可选

从中采样 QMC 点的 QMCEngine 实例。 QMCEngine 必须初始化为多个维度d与变量的数量相对应x1, ..., xd传递给函数。提供的 QMCEngine 用于生成第一个积分估计。如果n_estimates大于 1，则从第一个生成额外的 QMCEngine(如果有选项，则启用加扰。)如果未提供 QMCEngine，则使用默认值scipy.stats.qmc.Halton将使用根据长度确定的维数进行初始化a.

log： 布尔值，默认：False

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

返回 ：：

result： 对象

具有属性的结果对象：

不可缺少的 浮点数

积分的估计。

standard_error :

误差估计。解释见注释。

注意：

每个被积函数的值n_pointsQMC 样本的点用于生成积分的估计。该估计值是从积分的可能估计值总体中得出的，我们获得的值取决于评估积分的特定点。我们执行这个过程n_estimates次，每次评估不同扰乱 QMC 点处的被积函数，有效地绘制独立同分布。来自整体估计总体的随机样本。样本均值\(m\)这些积分估计值是积分真实值的无偏估计量，以及平均值的标准误差\(s\)这些估计值可用于使用 t 分布生成置信区间n_estimates - 1自由度。也许与直觉相反，增加n_points同时保留函数评价点的总数n_points * n_estimates固定往往会减少实际误差，而增加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