`scipy.integrate._quadrature:qmc_quad`

source: /scipy/integrate/_quadrature.py :1182

Signature

   def     qmc_quad (    func    ,    a    ,    b    ,  * ,     n_estimates    =  8   ,    n_points    =  1024   ,    qrng    =  None   ,    log    =  False     )

Summary

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

Parameters

func : callable: 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, b : array-like: One-dimensional arrays specifying the lower and upper integration limits, respectively, of each of the d variables.
n_estimates, n_points : int, 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.
qrng : `~scipy.stats.qmc.QMCEngine`, 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.
log : boolean, default: False: When set to True, func returns the log of the integrand, and the result object contains the log of the integral.

Returns

result : object

A result object with attributes:

integral: 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.

Array API Standard Support

qmc_quad has experimental support for Python Array API Standard compatible backends in addition to NumPy. Please consider testing these features by setting an environment variable SCIPY_ARRAY_API=1 and providing CuPy, PyTorch, JAX, or Dask arrays as array arguments. The following combinations of backend and device (or other capability) are supported.

====================  ====================  ====================
Library               CPU                   GPU
====================  ====================  ====================
NumPy                 ✅                     n/a                 
CuPy                  n/a                   ✅                   
PyTorch               ✅                     ✅                   
JAX                   ⚠️ no JIT             ⚠️ no JIT           
Dask                  ⛔                     n/a                 
====================  ====================  ====================

See dev-arrayapi for more information.

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

✗

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)

✗

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

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