Skip to main content
`Ctrl`+`K`
scipy.integrate.

# qmc_quad

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

  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 scalar or 1-d 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 `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 bool, 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:

      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.

  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 Support for the array API standard 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
      (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