scipy.integrate.

qmc_quad#

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

Compute an integral in N-dimensions using Quasi-Monte Carlo quadrature.

Parameters:

funccallable: The integrand. Must accept a single argument x, an array which specifies the point(s) at which to evaluate the scalar-valued integrand, and return the value(s) of the integrand. For efficiency, the function should be vectorized to accept an array of shape (d, n_points), where d is the number of variables (i.e. the dimensionality of the function domain) and n_points is the number of quadrature points, and return an array of shape (n_points,), the integrand at each quadrature point.
a, barray-like: One-dimensional arrays specifying the lower and upper integration limits, respectively, of each of the d variables.
n_estimates, n_pointsint, optional: n_estimates (default: 8) statistically independent QMC samples, each of n_points (default: 1024) points, will be generated by qrng. The total number of points at which the integrand func will be evaluated is n_points * n_estimates. See Notes for details.
qrngQMCEngine, optional: An instance of the QMCEngine from which to sample QMC points. The QMCEngine must be initialized to a number of dimensions d corresponding with the number of variables x1, ..., xd passed to func. The provided QMCEngine is used to produce the first integral estimate. If n_estimates is greater than one, additional QMCEngines are spawned from the first (with scrambling enabled, if it is an option.) If a QMCEngine is not provided, the default scipy.stats.qmc.Halton will be initialized with the number of dimensions determine from the length of a.
logboolean, default: False: When set to True, func returns the log of the integrand, and the result object contains the log of the integral.

Returns:

resultobject

A result object with attributes:

integralfloat: The estimate of the integral.
standard_error :: The error estimate. See Notes for interpretation.

Notes

Values of the integrand at each of the n_points points of a QMC sample are used to produce an estimate of the integral. This estimate is drawn from a population of possible estimates of the integral, the value of which we obtain depends on the particular points at which the integral was evaluated. We perform this process n_estimates times, each time evaluating the integrand at different scrambled QMC points, effectively drawing i.i.d. random samples from the population of integral estimates. The sample mean \(m\) of these integral estimates is an unbiased estimator of the true value of the integral, and the standard error of the mean \(s\) of these estimates may be used to generate confidence intervals using the t distribution with n_estimates - 1 degrees of freedom. Perhaps counter-intuitively, increasing n_points while keeping the total number of function evaluation points n_points * n_estimates fixed tends to reduce the actual error, whereas increasing n_estimates tends to decrease the error estimate.

Examples

QMC quadrature is particularly useful for computing integrals in higher dimensions. An example integrand is the probability density function of a multivariate normal distribution.

>>> 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)

To compute the integral over the unit hypercube:

>>> 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)

A two-sided, 99% confidence interval for the integral may be estimated as:

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

Indeed, the value reported by scipy.stats.multivariate_normal is within this range.

>>> stats.multivariate_normal.cdf(b, mean, cov, lower_limit=a)
0.00018430867675187443