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)
, whered
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.- qrng
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 variablesx1, ..., 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 defaultscipy.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 pointsn_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