qmc_quad#
- scipy.integrate.qmc_quad(func, a, b, *, n_estimates=8, n_points=1024, qrng=None, log=False)[源代码]#
使用准蒙特卡洛(Quasi-Monte Carlo)积分计算N维积分。
- 参数:
- func可调用对象
被积函数。必须接受一个参数
x
,一个指定用于评估标量被积函数的一个或多个点的数组,并返回被积函数的值。为了效率,函数应向量化以接受形状为(d, n_points)
的数组,其中d
是变量的数量(即函数域的维度),n_points 是积分点的数量,并返回一个形状为(n_points,)
的数组,即每个积分点处的被积函数值。- a, b类数组
一维数组,分别指定
d
个变量的积分下限和上限。- n_estimates, n_points整数,可选
n_estimates(默认:8)个统计独立的 QMC 样本,每个样本包含 n_points(默认:1024)个点,将由 qrng 生成。被积函数 func 将被评估的总点数为
n_points * n_estimates
。详见备注。- qrng
QMCEngine
, 可选 用于采样 QMC 点的 QMCEngine 实例。QMCEngine 必须初始化为与传递给 func 的变量数量
x1, ..., xd
对应的维度d
。所提供的 QMCEngine 用于生成第一个积分估计。如果 n_estimates 大于 1,则会从第一个 QMCEngine 衍生出额外的 QMCEngine(如果支持,会启用乱序)。如果未提供 QMCEngine,则默认的scipy.stats.qmc.Halton
将根据 a 的长度确定维度进行初始化。- log布尔值,默认值:False
当设置为 True 时,func 返回被积函数的对数,结果对象包含积分的对数。
- 返回:
- result对象
一个包含以下属性的结果对象
- integral浮点数
积分的估计值。
- 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