Coverage for /usr/lib/python3/dist-packages/scipy/integrate/_quadrature.py: 9%
444 statements
« prev ^ index » next coverage.py v7.9.1, created at 2025-06-14 15:55 +0200
1from __future__ import annotations
2from typing import TYPE_CHECKING, Callable, Any, cast
3import numpy as np
4import math
5import warnings
6from collections import namedtuple
8from scipy.special import roots_legendre
9from scipy.special import gammaln, logsumexp
10from scipy._lib._util import _rng_spawn
13__all__ = ['fixed_quad', 'quadrature', 'romberg', 'romb',
14 'trapezoid', 'trapz', 'simps', 'simpson',
15 'cumulative_trapezoid', 'cumtrapz', 'newton_cotes',
16 'qmc_quad', 'AccuracyWarning']
19def trapezoid(y, x=None, dx=1.0, axis=-1):
20 r"""
21 Integrate along the given axis using the composite trapezoidal rule.
23 If `x` is provided, the integration happens in sequence along its
24 elements - they are not sorted.
26 Integrate `y` (`x`) along each 1d slice on the given axis, compute
27 :math:`\int y(x) dx`.
28 When `x` is specified, this integrates along the parametric curve,
29 computing :math:`\int_t y(t) dt =
30 \int_t y(t) \left.\frac{dx}{dt}\right|_{x=x(t)} dt`.
32 Parameters
33 ----------
34 y : array_like
35 Input array to integrate.
36 x : array_like, optional
37 The sample points corresponding to the `y` values. If `x` is None,
38 the sample points are assumed to be evenly spaced `dx` apart. The
39 default is None.
40 dx : scalar, optional
41 The spacing between sample points when `x` is None. The default is 1.
42 axis : int, optional
43 The axis along which to integrate.
45 Returns
46 -------
47 trapezoid : float or ndarray
48 Definite integral of `y` = n-dimensional array as approximated along
49 a single axis by the trapezoidal rule. If `y` is a 1-dimensional array,
50 then the result is a float. If `n` is greater than 1, then the result
51 is an `n`-1 dimensional array.
53 See Also
54 --------
55 cumulative_trapezoid, simpson, romb
57 Notes
58 -----
59 Image [2]_ illustrates trapezoidal rule -- y-axis locations of points
60 will be taken from `y` array, by default x-axis distances between
61 points will be 1.0, alternatively they can be provided with `x` array
62 or with `dx` scalar. Return value will be equal to combined area under
63 the red lines.
65 References
66 ----------
67 .. [1] Wikipedia page: https://en.wikipedia.org/wiki/Trapezoidal_rule
69 .. [2] Illustration image:
70 https://en.wikipedia.org/wiki/File:Composite_trapezoidal_rule_illustration.png
72 Examples
73 --------
74 Use the trapezoidal rule on evenly spaced points:
76 >>> import numpy as np
77 >>> from scipy import integrate
78 >>> integrate.trapezoid([1, 2, 3])
79 4.0
81 The spacing between sample points can be selected by either the
82 ``x`` or ``dx`` arguments:
84 >>> integrate.trapezoid([1, 2, 3], x=[4, 6, 8])
85 8.0
86 >>> integrate.trapezoid([1, 2, 3], dx=2)
87 8.0
89 Using a decreasing ``x`` corresponds to integrating in reverse:
91 >>> integrate.trapezoid([1, 2, 3], x=[8, 6, 4])
92 -8.0
94 More generally ``x`` is used to integrate along a parametric curve. We can
95 estimate the integral :math:`\int_0^1 x^2 = 1/3` using:
97 >>> x = np.linspace(0, 1, num=50)
98 >>> y = x**2
99 >>> integrate.trapezoid(y, x)
100 0.33340274885464394
102 Or estimate the area of a circle, noting we repeat the sample which closes
103 the curve:
105 >>> theta = np.linspace(0, 2 * np.pi, num=1000, endpoint=True)
106 >>> integrate.trapezoid(np.cos(theta), x=np.sin(theta))
107 3.141571941375841
109 ``trapezoid`` can be applied along a specified axis to do multiple
110 computations in one call:
112 >>> a = np.arange(6).reshape(2, 3)
113 >>> a
114 array([[0, 1, 2],
115 [3, 4, 5]])
116 >>> integrate.trapezoid(a, axis=0)
117 array([1.5, 2.5, 3.5])
118 >>> integrate.trapezoid(a, axis=1)
119 array([2., 8.])
120 """
121 # Future-proofing, in case NumPy moves from trapz to trapezoid for the same
122 # reasons as SciPy
123 if hasattr(np, 'trapezoid'):
124 return np.trapezoid(y, x=x, dx=dx, axis=axis)
125 else:
126 return np.trapz(y, x=x, dx=dx, axis=axis)
129# Note: alias kept for backwards compatibility. Rename was done
130# because trapz is a slur in colloquial English (see gh-12924).
131def trapz(y, x=None, dx=1.0, axis=-1):
132 """An alias of `trapezoid`.
134 `trapz` is kept for backwards compatibility. For new code, prefer
135 `trapezoid` instead.
136 """
137 return trapezoid(y, x=x, dx=dx, axis=axis)
140class AccuracyWarning(Warning):
141 pass
144if TYPE_CHECKING:
145 # workaround for mypy function attributes see:
146 # https://github.com/python/mypy/issues/2087#issuecomment-462726600
147 from typing import Protocol
149 class CacheAttributes(Protocol):
150 cache: dict[int, tuple[Any, Any]]
151else:
152 CacheAttributes = Callable
155def cache_decorator(func: Callable) -> CacheAttributes:
156 return cast(CacheAttributes, func)
159@cache_decorator
160def _cached_roots_legendre(n):
161 """
162 Cache roots_legendre results to speed up calls of the fixed_quad
163 function.
164 """
165 if n in _cached_roots_legendre.cache:
166 return _cached_roots_legendre.cache[n]
168 _cached_roots_legendre.cache[n] = roots_legendre(n)
169 return _cached_roots_legendre.cache[n]
172_cached_roots_legendre.cache = dict()
175def fixed_quad(func, a, b, args=(), n=5):
176 """
177 Compute a definite integral using fixed-order Gaussian quadrature.
179 Integrate `func` from `a` to `b` using Gaussian quadrature of
180 order `n`.
182 Parameters
183 ----------
184 func : callable
185 A Python function or method to integrate (must accept vector inputs).
186 If integrating a vector-valued function, the returned array must have
187 shape ``(..., len(x))``.
188 a : float
189 Lower limit of integration.
190 b : float
191 Upper limit of integration.
192 args : tuple, optional
193 Extra arguments to pass to function, if any.
194 n : int, optional
195 Order of quadrature integration. Default is 5.
197 Returns
198 -------
199 val : float
200 Gaussian quadrature approximation to the integral
201 none : None
202 Statically returned value of None
204 See Also
205 --------
206 quad : adaptive quadrature using QUADPACK
207 dblquad : double integrals
208 tplquad : triple integrals
209 romberg : adaptive Romberg quadrature
210 quadrature : adaptive Gaussian quadrature
211 romb : integrators for sampled data
212 simpson : integrators for sampled data
213 cumulative_trapezoid : cumulative integration for sampled data
214 ode : ODE integrator
215 odeint : ODE integrator
217 Examples
218 --------
219 >>> from scipy import integrate
220 >>> import numpy as np
221 >>> f = lambda x: x**8
222 >>> integrate.fixed_quad(f, 0.0, 1.0, n=4)
223 (0.1110884353741496, None)
224 >>> integrate.fixed_quad(f, 0.0, 1.0, n=5)
225 (0.11111111111111102, None)
226 >>> print(1/9.0) # analytical result
227 0.1111111111111111
229 >>> integrate.fixed_quad(np.cos, 0.0, np.pi/2, n=4)
230 (0.9999999771971152, None)
231 >>> integrate.fixed_quad(np.cos, 0.0, np.pi/2, n=5)
232 (1.000000000039565, None)
233 >>> np.sin(np.pi/2)-np.sin(0) # analytical result
234 1.0
236 """
237 x, w = _cached_roots_legendre(n)
238 x = np.real(x)
239 if np.isinf(a) or np.isinf(b):
240 raise ValueError("Gaussian quadrature is only available for "
241 "finite limits.")
242 y = (b-a)*(x+1)/2.0 + a
243 return (b-a)/2.0 * np.sum(w*func(y, *args), axis=-1), None
246def vectorize1(func, args=(), vec_func=False):
247 """Vectorize the call to a function.
249 This is an internal utility function used by `romberg` and
250 `quadrature` to create a vectorized version of a function.
252 If `vec_func` is True, the function `func` is assumed to take vector
253 arguments.
255 Parameters
256 ----------
257 func : callable
258 User defined function.
259 args : tuple, optional
260 Extra arguments for the function.
261 vec_func : bool, optional
262 True if the function func takes vector arguments.
264 Returns
265 -------
266 vfunc : callable
267 A function that will take a vector argument and return the
268 result.
270 """
271 if vec_func:
272 def vfunc(x):
273 return func(x, *args)
274 else:
275 def vfunc(x):
276 if np.isscalar(x):
277 return func(x, *args)
278 x = np.asarray(x)
279 # call with first point to get output type
280 y0 = func(x[0], *args)
281 n = len(x)
282 dtype = getattr(y0, 'dtype', type(y0))
283 output = np.empty((n,), dtype=dtype)
284 output[0] = y0
285 for i in range(1, n):
286 output[i] = func(x[i], *args)
287 return output
288 return vfunc
291def quadrature(func, a, b, args=(), tol=1.49e-8, rtol=1.49e-8, maxiter=50,
292 vec_func=True, miniter=1):
293 """
294 Compute a definite integral using fixed-tolerance Gaussian quadrature.
296 Integrate `func` from `a` to `b` using Gaussian quadrature
297 with absolute tolerance `tol`.
299 Parameters
300 ----------
301 func : function
302 A Python function or method to integrate.
303 a : float
304 Lower limit of integration.
305 b : float
306 Upper limit of integration.
307 args : tuple, optional
308 Extra arguments to pass to function.
309 tol, rtol : float, optional
310 Iteration stops when error between last two iterates is less than
311 `tol` OR the relative change is less than `rtol`.
312 maxiter : int, optional
313 Maximum order of Gaussian quadrature.
314 vec_func : bool, optional
315 True or False if func handles arrays as arguments (is
316 a "vector" function). Default is True.
317 miniter : int, optional
318 Minimum order of Gaussian quadrature.
320 Returns
321 -------
322 val : float
323 Gaussian quadrature approximation (within tolerance) to integral.
324 err : float
325 Difference between last two estimates of the integral.
327 See Also
328 --------
329 romberg : adaptive Romberg quadrature
330 fixed_quad : fixed-order Gaussian quadrature
331 quad : adaptive quadrature using QUADPACK
332 dblquad : double integrals
333 tplquad : triple integrals
334 romb : integrator for sampled data
335 simpson : integrator for sampled data
336 cumulative_trapezoid : cumulative integration for sampled data
337 ode : ODE integrator
338 odeint : ODE integrator
340 Examples
341 --------
342 >>> from scipy import integrate
343 >>> import numpy as np
344 >>> f = lambda x: x**8
345 >>> integrate.quadrature(f, 0.0, 1.0)
346 (0.11111111111111106, 4.163336342344337e-17)
347 >>> print(1/9.0) # analytical result
348 0.1111111111111111
350 >>> integrate.quadrature(np.cos, 0.0, np.pi/2)
351 (0.9999999999999536, 3.9611425250996035e-11)
352 >>> np.sin(np.pi/2)-np.sin(0) # analytical result
353 1.0
355 """
356 if not isinstance(args, tuple):
357 args = (args,)
358 vfunc = vectorize1(func, args, vec_func=vec_func)
359 val = np.inf
360 err = np.inf
361 maxiter = max(miniter+1, maxiter)
362 for n in range(miniter, maxiter+1):
363 newval = fixed_quad(vfunc, a, b, (), n)[0]
364 err = abs(newval-val)
365 val = newval
367 if err < tol or err < rtol*abs(val):
368 break
369 else:
370 warnings.warn(
371 "maxiter (%d) exceeded. Latest difference = %e" % (maxiter, err),
372 AccuracyWarning)
373 return val, err
376def tupleset(t, i, value):
377 l = list(t)
378 l[i] = value
379 return tuple(l)
382# Note: alias kept for backwards compatibility. Rename was done
383# because cumtrapz is a slur in colloquial English (see gh-12924).
384def cumtrapz(y, x=None, dx=1.0, axis=-1, initial=None):
385 """An alias of `cumulative_trapezoid`.
387 `cumtrapz` is kept for backwards compatibility. For new code, prefer
388 `cumulative_trapezoid` instead.
389 """
390 return cumulative_trapezoid(y, x=x, dx=dx, axis=axis, initial=initial)
393def cumulative_trapezoid(y, x=None, dx=1.0, axis=-1, initial=None):
394 """
395 Cumulatively integrate y(x) using the composite trapezoidal rule.
397 Parameters
398 ----------
399 y : array_like
400 Values to integrate.
401 x : array_like, optional
402 The coordinate to integrate along. If None (default), use spacing `dx`
403 between consecutive elements in `y`.
404 dx : float, optional
405 Spacing between elements of `y`. Only used if `x` is None.
406 axis : int, optional
407 Specifies the axis to cumulate. Default is -1 (last axis).
408 initial : scalar, optional
409 If given, insert this value at the beginning of the returned result.
410 Typically this value should be 0. Default is None, which means no
411 value at ``x[0]`` is returned and `res` has one element less than `y`
412 along the axis of integration.
414 Returns
415 -------
416 res : ndarray
417 The result of cumulative integration of `y` along `axis`.
418 If `initial` is None, the shape is such that the axis of integration
419 has one less value than `y`. If `initial` is given, the shape is equal
420 to that of `y`.
422 See Also
423 --------
424 numpy.cumsum, numpy.cumprod
425 quad : adaptive quadrature using QUADPACK
426 romberg : adaptive Romberg quadrature
427 quadrature : adaptive Gaussian quadrature
428 fixed_quad : fixed-order Gaussian quadrature
429 dblquad : double integrals
430 tplquad : triple integrals
431 romb : integrators for sampled data
432 ode : ODE integrators
433 odeint : ODE integrators
435 Examples
436 --------
437 >>> from scipy import integrate
438 >>> import numpy as np
439 >>> import matplotlib.pyplot as plt
441 >>> x = np.linspace(-2, 2, num=20)
442 >>> y = x
443 >>> y_int = integrate.cumulative_trapezoid(y, x, initial=0)
444 >>> plt.plot(x, y_int, 'ro', x, y[0] + 0.5 * x**2, 'b-')
445 >>> plt.show()
447 """
448 y = np.asarray(y)
449 if x is None:
450 d = dx
451 else:
452 x = np.asarray(x)
453 if x.ndim == 1:
454 d = np.diff(x)
455 # reshape to correct shape
456 shape = [1] * y.ndim
457 shape[axis] = -1
458 d = d.reshape(shape)
459 elif len(x.shape) != len(y.shape):
460 raise ValueError("If given, shape of x must be 1-D or the "
461 "same as y.")
462 else:
463 d = np.diff(x, axis=axis)
465 if d.shape[axis] != y.shape[axis] - 1:
466 raise ValueError("If given, length of x along axis must be the "
467 "same as y.")
469 nd = len(y.shape)
470 slice1 = tupleset((slice(None),)*nd, axis, slice(1, None))
471 slice2 = tupleset((slice(None),)*nd, axis, slice(None, -1))
472 res = np.cumsum(d * (y[slice1] + y[slice2]) / 2.0, axis=axis)
474 if initial is not None:
475 if not np.isscalar(initial):
476 raise ValueError("`initial` parameter should be a scalar.")
478 shape = list(res.shape)
479 shape[axis] = 1
480 res = np.concatenate([np.full(shape, initial, dtype=res.dtype), res],
481 axis=axis)
483 return res
486def _basic_simpson(y, start, stop, x, dx, axis):
487 nd = len(y.shape)
488 if start is None:
489 start = 0
490 step = 2
491 slice_all = (slice(None),)*nd
492 slice0 = tupleset(slice_all, axis, slice(start, stop, step))
493 slice1 = tupleset(slice_all, axis, slice(start+1, stop+1, step))
494 slice2 = tupleset(slice_all, axis, slice(start+2, stop+2, step))
496 if x is None: # Even-spaced Simpson's rule.
497 result = np.sum(y[slice0] + 4.0*y[slice1] + y[slice2], axis=axis)
498 result *= dx / 3.0
499 else:
500 # Account for possibly different spacings.
501 # Simpson's rule changes a bit.
502 h = np.diff(x, axis=axis)
503 sl0 = tupleset(slice_all, axis, slice(start, stop, step))
504 sl1 = tupleset(slice_all, axis, slice(start+1, stop+1, step))
505 h0 = h[sl0].astype(float, copy=False)
506 h1 = h[sl1].astype(float, copy=False)
507 hsum = h0 + h1
508 hprod = h0 * h1
509 h0divh1 = np.true_divide(h0, h1, out=np.zeros_like(h0), where=h1 != 0)
510 tmp = hsum/6.0 * (y[slice0] *
511 (2.0 - np.true_divide(1.0, h0divh1,
512 out=np.zeros_like(h0divh1),
513 where=h0divh1 != 0)) +
514 y[slice1] * (hsum *
515 np.true_divide(hsum, hprod,
516 out=np.zeros_like(hsum),
517 where=hprod != 0)) +
518 y[slice2] * (2.0 - h0divh1))
519 result = np.sum(tmp, axis=axis)
520 return result
523# Note: alias kept for backwards compatibility. simps was renamed to simpson
524# because the former is a slur in colloquial English (see gh-12924).
525def simps(y, x=None, dx=1.0, axis=-1, even=None):
526 """An alias of `simpson`.
528 `simps` is kept for backwards compatibility. For new code, prefer
529 `simpson` instead.
530 """
531 return simpson(y, x=x, dx=dx, axis=axis, even=even)
534def simpson(y, x=None, dx=1.0, axis=-1, even=None):
535 """
536 Integrate y(x) using samples along the given axis and the composite
537 Simpson's rule. If x is None, spacing of dx is assumed.
539 If there are an even number of samples, N, then there are an odd
540 number of intervals (N-1), but Simpson's rule requires an even number
541 of intervals. The parameter 'even' controls how this is handled.
543 Parameters
544 ----------
545 y : array_like
546 Array to be integrated.
547 x : array_like, optional
548 If given, the points at which `y` is sampled.
549 dx : float, optional
550 Spacing of integration points along axis of `x`. Only used when
551 `x` is None. Default is 1.
552 axis : int, optional
553 Axis along which to integrate. Default is the last axis.
554 even : {None, 'simpson', 'avg', 'first', 'last'}, optional
555 'avg' : Average two results:
556 1) use the first N-2 intervals with
557 a trapezoidal rule on the last interval and
558 2) use the last
559 N-2 intervals with a trapezoidal rule on the first interval.
561 'first' : Use Simpson's rule for the first N-2 intervals with
562 a trapezoidal rule on the last interval.
564 'last' : Use Simpson's rule for the last N-2 intervals with a
565 trapezoidal rule on the first interval.
567 None : equivalent to 'simpson' (default)
569 'simpson' : Use Simpson's rule for the first N-2 intervals with the
570 addition of a 3-point parabolic segment for the last
571 interval using equations outlined by Cartwright [1]_.
572 If the axis to be integrated over only has two points then
573 the integration falls back to a trapezoidal integration.
575 .. versionadded:: 1.11.0
577 .. versionchanged:: 1.11.0
578 The newly added 'simpson' option is now the default as it is more
579 accurate in most situations.
581 .. deprecated:: 1.11.0
582 Parameter `even` is deprecated and will be removed in SciPy
583 1.13.0. After this time the behaviour for an even number of
584 points will follow that of `even='simpson'`.
586 Returns
587 -------
588 float
589 The estimated integral computed with the composite Simpson's rule.
591 See Also
592 --------
593 quad : adaptive quadrature using QUADPACK
594 romberg : adaptive Romberg quadrature
595 quadrature : adaptive Gaussian quadrature
596 fixed_quad : fixed-order Gaussian quadrature
597 dblquad : double integrals
598 tplquad : triple integrals
599 romb : integrators for sampled data
600 cumulative_trapezoid : cumulative integration for sampled data
601 ode : ODE integrators
602 odeint : ODE integrators
604 Notes
605 -----
606 For an odd number of samples that are equally spaced the result is
607 exact if the function is a polynomial of order 3 or less. If
608 the samples are not equally spaced, then the result is exact only
609 if the function is a polynomial of order 2 or less.
611 References
612 ----------
613 .. [1] Cartwright, Kenneth V. Simpson's Rule Cumulative Integration with
614 MS Excel and Irregularly-spaced Data. Journal of Mathematical
615 Sciences and Mathematics Education. 12 (2): 1-9
617 Examples
618 --------
619 >>> from scipy import integrate
620 >>> import numpy as np
621 >>> x = np.arange(0, 10)
622 >>> y = np.arange(0, 10)
624 >>> integrate.simpson(y, x)
625 40.5
627 >>> y = np.power(x, 3)
628 >>> integrate.simpson(y, x)
629 1640.5
630 >>> integrate.quad(lambda x: x**3, 0, 9)[0]
631 1640.25
633 >>> integrate.simpson(y, x, even='first')
634 1644.5
636 """
637 y = np.asarray(y)
638 nd = len(y.shape)
639 N = y.shape[axis]
640 last_dx = dx
641 first_dx = dx
642 returnshape = 0
643 if x is not None:
644 x = np.asarray(x)
645 if len(x.shape) == 1:
646 shapex = [1] * nd
647 shapex[axis] = x.shape[0]
648 saveshape = x.shape
649 returnshape = 1
650 x = x.reshape(tuple(shapex))
651 elif len(x.shape) != len(y.shape):
652 raise ValueError("If given, shape of x must be 1-D or the "
653 "same as y.")
654 if x.shape[axis] != N:
655 raise ValueError("If given, length of x along axis must be the "
656 "same as y.")
658 # even keyword parameter is deprecated
659 if even is not None:
660 warnings.warn(
661 "The 'even' keyword is deprecated as of SciPy 1.11.0 and will be "
662 "removed in SciPy 1.13.0",
663 DeprecationWarning, stacklevel=2
664 )
666 if N % 2 == 0:
667 val = 0.0
668 result = 0.0
669 slice_all = (slice(None),) * nd
671 # default is 'simpson'
672 even = even if even is not None else "simpson"
674 if even not in ['avg', 'last', 'first', 'simpson']:
675 raise ValueError(
676 "Parameter 'even' must be 'simpson', "
677 "'avg', 'last', or 'first'."
678 )
680 if N == 2:
681 # need at least 3 points in integration axis to form parabolic
682 # segment. If there are two points then any of 'avg', 'first',
683 # 'last' should give the same result.
684 slice1 = tupleset(slice_all, axis, -1)
685 slice2 = tupleset(slice_all, axis, -2)
686 if x is not None:
687 last_dx = x[slice1] - x[slice2]
688 val += 0.5 * last_dx * (y[slice1] + y[slice2])
690 # calculation is finished. Set `even` to None to skip other
691 # scenarios
692 even = None
694 if even == 'simpson':
695 # use Simpson's rule on first intervals
696 result = _basic_simpson(y, 0, N-3, x, dx, axis)
698 slice1 = tupleset(slice_all, axis, -1)
699 slice2 = tupleset(slice_all, axis, -2)
700 slice3 = tupleset(slice_all, axis, -3)
702 h = np.asfarray([dx, dx])
703 if x is not None:
704 # grab the last two spacings from the appropriate axis
705 hm2 = tupleset(slice_all, axis, slice(-2, -1, 1))
706 hm1 = tupleset(slice_all, axis, slice(-1, None, 1))
708 diffs = np.float64(np.diff(x, axis=axis))
709 h = [np.squeeze(diffs[hm2], axis=axis),
710 np.squeeze(diffs[hm1], axis=axis)]
712 # This is the correction for the last interval according to
713 # Cartwright.
714 # However, I used the equations given at
715 # https://en.wikipedia.org/wiki/Simpson%27s_rule#Composite_Simpson's_rule_for_irregularly_spaced_data
716 # A footnote on Wikipedia says:
717 # Cartwright 2017, Equation 8. The equation in Cartwright is
718 # calculating the first interval whereas the equations in the
719 # Wikipedia article are adjusting for the last integral. If the
720 # proper algebraic substitutions are made, the equation results in
721 # the values shown.
722 num = 2 * h[1] ** 2 + 3 * h[0] * h[1]
723 den = 6 * (h[1] + h[0])
724 alpha = np.true_divide(
725 num,
726 den,
727 out=np.zeros_like(den),
728 where=den != 0
729 )
731 num = h[1] ** 2 + 3.0 * h[0] * h[1]
732 den = 6 * h[0]
733 beta = np.true_divide(
734 num,
735 den,
736 out=np.zeros_like(den),
737 where=den != 0
738 )
740 num = 1 * h[1] ** 3
741 den = 6 * h[0] * (h[0] + h[1])
742 eta = np.true_divide(
743 num,
744 den,
745 out=np.zeros_like(den),
746 where=den != 0
747 )
749 result += alpha*y[slice1] + beta*y[slice2] - eta*y[slice3]
751 # The following code (down to result=result+val) can be removed
752 # once the 'even' keyword is removed.
754 # Compute using Simpson's rule on first intervals
755 if even in ['avg', 'first']:
756 slice1 = tupleset(slice_all, axis, -1)
757 slice2 = tupleset(slice_all, axis, -2)
758 if x is not None:
759 last_dx = x[slice1] - x[slice2]
760 val += 0.5*last_dx*(y[slice1]+y[slice2])
761 result = _basic_simpson(y, 0, N-3, x, dx, axis)
762 # Compute using Simpson's rule on last set of intervals
763 if even in ['avg', 'last']:
764 slice1 = tupleset(slice_all, axis, 0)
765 slice2 = tupleset(slice_all, axis, 1)
766 if x is not None:
767 first_dx = x[tuple(slice2)] - x[tuple(slice1)]
768 val += 0.5*first_dx*(y[slice2]+y[slice1])
769 result += _basic_simpson(y, 1, N-2, x, dx, axis)
770 if even == 'avg':
771 val /= 2.0
772 result /= 2.0
773 result = result + val
774 else:
775 result = _basic_simpson(y, 0, N-2, x, dx, axis)
776 if returnshape:
777 x = x.reshape(saveshape)
778 return result
781def romb(y, dx=1.0, axis=-1, show=False):
782 """
783 Romberg integration using samples of a function.
785 Parameters
786 ----------
787 y : array_like
788 A vector of ``2**k + 1`` equally-spaced samples of a function.
789 dx : float, optional
790 The sample spacing. Default is 1.
791 axis : int, optional
792 The axis along which to integrate. Default is -1 (last axis).
793 show : bool, optional
794 When `y` is a single 1-D array, then if this argument is True
795 print the table showing Richardson extrapolation from the
796 samples. Default is False.
798 Returns
799 -------
800 romb : ndarray
801 The integrated result for `axis`.
803 See Also
804 --------
805 quad : adaptive quadrature using QUADPACK
806 romberg : adaptive Romberg quadrature
807 quadrature : adaptive Gaussian quadrature
808 fixed_quad : fixed-order Gaussian quadrature
809 dblquad : double integrals
810 tplquad : triple integrals
811 simpson : integrators for sampled data
812 cumulative_trapezoid : cumulative integration for sampled data
813 ode : ODE integrators
814 odeint : ODE integrators
816 Examples
817 --------
818 >>> from scipy import integrate
819 >>> import numpy as np
820 >>> x = np.arange(10, 14.25, 0.25)
821 >>> y = np.arange(3, 12)
823 >>> integrate.romb(y)
824 56.0
826 >>> y = np.sin(np.power(x, 2.5))
827 >>> integrate.romb(y)
828 -0.742561336672229
830 >>> integrate.romb(y, show=True)
831 Richardson Extrapolation Table for Romberg Integration
832 ======================================================
833 -0.81576
834 4.63862 6.45674
835 -1.10581 -3.02062 -3.65245
836 -2.57379 -3.06311 -3.06595 -3.05664
837 -1.34093 -0.92997 -0.78776 -0.75160 -0.74256
838 ======================================================
839 -0.742561336672229 # may vary
841 """
842 y = np.asarray(y)
843 nd = len(y.shape)
844 Nsamps = y.shape[axis]
845 Ninterv = Nsamps-1
846 n = 1
847 k = 0
848 while n < Ninterv:
849 n <<= 1
850 k += 1
851 if n != Ninterv:
852 raise ValueError("Number of samples must be one plus a "
853 "non-negative power of 2.")
855 R = {}
856 slice_all = (slice(None),) * nd
857 slice0 = tupleset(slice_all, axis, 0)
858 slicem1 = tupleset(slice_all, axis, -1)
859 h = Ninterv * np.asarray(dx, dtype=float)
860 R[(0, 0)] = (y[slice0] + y[slicem1])/2.0*h
861 slice_R = slice_all
862 start = stop = step = Ninterv
863 for i in range(1, k+1):
864 start >>= 1
865 slice_R = tupleset(slice_R, axis, slice(start, stop, step))
866 step >>= 1
867 R[(i, 0)] = 0.5*(R[(i-1, 0)] + h*y[slice_R].sum(axis=axis))
868 for j in range(1, i+1):
869 prev = R[(i, j-1)]
870 R[(i, j)] = prev + (prev-R[(i-1, j-1)]) / ((1 << (2*j))-1)
871 h /= 2.0
873 if show:
874 if not np.isscalar(R[(0, 0)]):
875 print("*** Printing table only supported for integrals" +
876 " of a single data set.")
877 else:
878 try:
879 precis = show[0]
880 except (TypeError, IndexError):
881 precis = 5
882 try:
883 width = show[1]
884 except (TypeError, IndexError):
885 width = 8
886 formstr = "%%%d.%df" % (width, precis)
888 title = "Richardson Extrapolation Table for Romberg Integration"
889 print(title, "=" * len(title), sep="\n", end="\n")
890 for i in range(k+1):
891 for j in range(i+1):
892 print(formstr % R[(i, j)], end=" ")
893 print()
894 print("=" * len(title))
896 return R[(k, k)]
898# Romberg quadratures for numeric integration.
899#
900# Written by Scott M. Ransom <ransom@cfa.harvard.edu>
901# last revision: 14 Nov 98
902#
903# Cosmetic changes by Konrad Hinsen <hinsen@cnrs-orleans.fr>
904# last revision: 1999-7-21
905#
906# Adapted to SciPy by Travis Oliphant <oliphant.travis@ieee.org>
907# last revision: Dec 2001
910def _difftrap(function, interval, numtraps):
911 """
912 Perform part of the trapezoidal rule to integrate a function.
913 Assume that we had called difftrap with all lower powers-of-2
914 starting with 1. Calling difftrap only returns the summation
915 of the new ordinates. It does _not_ multiply by the width
916 of the trapezoids. This must be performed by the caller.
917 'function' is the function to evaluate (must accept vector arguments).
918 'interval' is a sequence with lower and upper limits
919 of integration.
920 'numtraps' is the number of trapezoids to use (must be a
921 power-of-2).
922 """
923 if numtraps <= 0:
924 raise ValueError("numtraps must be > 0 in difftrap().")
925 elif numtraps == 1:
926 return 0.5*(function(interval[0])+function(interval[1]))
927 else:
928 numtosum = numtraps/2
929 h = float(interval[1]-interval[0])/numtosum
930 lox = interval[0] + 0.5 * h
931 points = lox + h * np.arange(numtosum)
932 s = np.sum(function(points), axis=0)
933 return s
936def _romberg_diff(b, c, k):
937 """
938 Compute the differences for the Romberg quadrature corrections.
939 See Forman Acton's "Real Computing Made Real," p 143.
940 """
941 tmp = 4.0**k
942 return (tmp * c - b)/(tmp - 1.0)
945def _printresmat(function, interval, resmat):
946 # Print the Romberg result matrix.
947 i = j = 0
948 print('Romberg integration of', repr(function), end=' ')
949 print('from', interval)
950 print('')
951 print('%6s %9s %9s' % ('Steps', 'StepSize', 'Results'))
952 for i in range(len(resmat)):
953 print('%6d %9f' % (2**i, (interval[1]-interval[0])/(2.**i)), end=' ')
954 for j in range(i+1):
955 print('%9f' % (resmat[i][j]), end=' ')
956 print('')
957 print('')
958 print('The final result is', resmat[i][j], end=' ')
959 print('after', 2**(len(resmat)-1)+1, 'function evaluations.')
962def romberg(function, a, b, args=(), tol=1.48e-8, rtol=1.48e-8, show=False,
963 divmax=10, vec_func=False):
964 """
965 Romberg integration of a callable function or method.
967 Returns the integral of `function` (a function of one variable)
968 over the interval (`a`, `b`).
970 If `show` is 1, the triangular array of the intermediate results
971 will be printed. If `vec_func` is True (default is False), then
972 `function` is assumed to support vector arguments.
974 Parameters
975 ----------
976 function : callable
977 Function to be integrated.
978 a : float
979 Lower limit of integration.
980 b : float
981 Upper limit of integration.
983 Returns
984 -------
985 results : float
986 Result of the integration.
988 Other Parameters
989 ----------------
990 args : tuple, optional
991 Extra arguments to pass to function. Each element of `args` will
992 be passed as a single argument to `func`. Default is to pass no
993 extra arguments.
994 tol, rtol : float, optional
995 The desired absolute and relative tolerances. Defaults are 1.48e-8.
996 show : bool, optional
997 Whether to print the results. Default is False.
998 divmax : int, optional
999 Maximum order of extrapolation. Default is 10.
1000 vec_func : bool, optional
1001 Whether `func` handles arrays as arguments (i.e., whether it is a
1002 "vector" function). Default is False.
1004 See Also
1005 --------
1006 fixed_quad : Fixed-order Gaussian quadrature.
1007 quad : Adaptive quadrature using QUADPACK.
1008 dblquad : Double integrals.
1009 tplquad : Triple integrals.
1010 romb : Integrators for sampled data.
1011 simpson : Integrators for sampled data.
1012 cumulative_trapezoid : Cumulative integration for sampled data.
1013 ode : ODE integrator.
1014 odeint : ODE integrator.
1016 References
1017 ----------
1018 .. [1] 'Romberg's method' https://en.wikipedia.org/wiki/Romberg%27s_method
1020 Examples
1021 --------
1022 Integrate a gaussian from 0 to 1 and compare to the error function.
1024 >>> from scipy import integrate
1025 >>> from scipy.special import erf
1026 >>> import numpy as np
1027 >>> gaussian = lambda x: 1/np.sqrt(np.pi) * np.exp(-x**2)
1028 >>> result = integrate.romberg(gaussian, 0, 1, show=True)
1029 Romberg integration of <function vfunc at ...> from [0, 1]
1031 ::
1033 Steps StepSize Results
1034 1 1.000000 0.385872
1035 2 0.500000 0.412631 0.421551
1036 4 0.250000 0.419184 0.421368 0.421356
1037 8 0.125000 0.420810 0.421352 0.421350 0.421350
1038 16 0.062500 0.421215 0.421350 0.421350 0.421350 0.421350
1039 32 0.031250 0.421317 0.421350 0.421350 0.421350 0.421350 0.421350
1041 The final result is 0.421350396475 after 33 function evaluations.
1043 >>> print("%g %g" % (2*result, erf(1)))
1044 0.842701 0.842701
1046 """
1047 if np.isinf(a) or np.isinf(b):
1048 raise ValueError("Romberg integration only available "
1049 "for finite limits.")
1050 vfunc = vectorize1(function, args, vec_func=vec_func)
1051 n = 1
1052 interval = [a, b]
1053 intrange = b - a
1054 ordsum = _difftrap(vfunc, interval, n)
1055 result = intrange * ordsum
1056 resmat = [[result]]
1057 err = np.inf
1058 last_row = resmat[0]
1059 for i in range(1, divmax+1):
1060 n *= 2
1061 ordsum += _difftrap(vfunc, interval, n)
1062 row = [intrange * ordsum / n]
1063 for k in range(i):
1064 row.append(_romberg_diff(last_row[k], row[k], k+1))
1065 result = row[i]
1066 lastresult = last_row[i-1]
1067 if show:
1068 resmat.append(row)
1069 err = abs(result - lastresult)
1070 if err < tol or err < rtol * abs(result):
1071 break
1072 last_row = row
1073 else:
1074 warnings.warn(
1075 "divmax (%d) exceeded. Latest difference = %e" % (divmax, err),
1076 AccuracyWarning)
1078 if show:
1079 _printresmat(vfunc, interval, resmat)
1080 return result
1083# Coefficients for Newton-Cotes quadrature
1084#
1085# These are the points being used
1086# to construct the local interpolating polynomial
1087# a are the weights for Newton-Cotes integration
1088# B is the error coefficient.
1089# error in these coefficients grows as N gets larger.
1090# or as samples are closer and closer together
1092# You can use maxima to find these rational coefficients
1093# for equally spaced data using the commands
1094# a(i,N) := integrate(product(r-j,j,0,i-1) * product(r-j,j,i+1,N),r,0,N) / ((N-i)! * i!) * (-1)^(N-i);
1095# Be(N) := N^(N+2)/(N+2)! * (N/(N+3) - sum((i/N)^(N+2)*a(i,N),i,0,N));
1096# Bo(N) := N^(N+1)/(N+1)! * (N/(N+2) - sum((i/N)^(N+1)*a(i,N),i,0,N));
1097# B(N) := (if (mod(N,2)=0) then Be(N) else Bo(N));
1098#
1099# pre-computed for equally-spaced weights
1100#
1101# num_a, den_a, int_a, num_B, den_B = _builtincoeffs[N]
1102#
1103# a = num_a*array(int_a)/den_a
1104# B = num_B*1.0 / den_B
1105#
1106# integrate(f(x),x,x_0,x_N) = dx*sum(a*f(x_i)) + B*(dx)^(2k+3) f^(2k+2)(x*)
1107# where k = N // 2
1108#
1109_builtincoeffs = {
1110 1: (1,2,[1,1],-1,12),
1111 2: (1,3,[1,4,1],-1,90),
1112 3: (3,8,[1,3,3,1],-3,80),
1113 4: (2,45,[7,32,12,32,7],-8,945),
1114 5: (5,288,[19,75,50,50,75,19],-275,12096),
1115 6: (1,140,[41,216,27,272,27,216,41],-9,1400),
1116 7: (7,17280,[751,3577,1323,2989,2989,1323,3577,751],-8183,518400),
1117 8: (4,14175,[989,5888,-928,10496,-4540,10496,-928,5888,989],
1118 -2368,467775),
1119 9: (9,89600,[2857,15741,1080,19344,5778,5778,19344,1080,
1120 15741,2857], -4671, 394240),
1121 10: (5,299376,[16067,106300,-48525,272400,-260550,427368,
1122 -260550,272400,-48525,106300,16067],
1123 -673175, 163459296),
1124 11: (11,87091200,[2171465,13486539,-3237113, 25226685,-9595542,
1125 15493566,15493566,-9595542,25226685,-3237113,
1126 13486539,2171465], -2224234463, 237758976000),
1127 12: (1, 5255250, [1364651,9903168,-7587864,35725120,-51491295,
1128 87516288,-87797136,87516288,-51491295,35725120,
1129 -7587864,9903168,1364651], -3012, 875875),
1130 13: (13, 402361344000,[8181904909, 56280729661, -31268252574,
1131 156074417954,-151659573325,206683437987,
1132 -43111992612,-43111992612,206683437987,
1133 -151659573325,156074417954,-31268252574,
1134 56280729661,8181904909], -2639651053,
1135 344881152000),
1136 14: (7, 2501928000, [90241897,710986864,-770720657,3501442784,
1137 -6625093363,12630121616,-16802270373,19534438464,
1138 -16802270373,12630121616,-6625093363,3501442784,
1139 -770720657,710986864,90241897], -3740727473,
1140 1275983280000)
1141 }
1144def newton_cotes(rn, equal=0):
1145 r"""
1146 Return weights and error coefficient for Newton-Cotes integration.
1148 Suppose we have (N+1) samples of f at the positions
1149 x_0, x_1, ..., x_N. Then an N-point Newton-Cotes formula for the
1150 integral between x_0 and x_N is:
1152 :math:`\int_{x_0}^{x_N} f(x)dx = \Delta x \sum_{i=0}^{N} a_i f(x_i)
1153 + B_N (\Delta x)^{N+2} f^{N+1} (\xi)`
1155 where :math:`\xi \in [x_0,x_N]`
1156 and :math:`\Delta x = \frac{x_N-x_0}{N}` is the average samples spacing.
1158 If the samples are equally-spaced and N is even, then the error
1159 term is :math:`B_N (\Delta x)^{N+3} f^{N+2}(\xi)`.
1161 Parameters
1162 ----------
1163 rn : int
1164 The integer order for equally-spaced data or the relative positions of
1165 the samples with the first sample at 0 and the last at N, where N+1 is
1166 the length of `rn`. N is the order of the Newton-Cotes integration.
1167 equal : int, optional
1168 Set to 1 to enforce equally spaced data.
1170 Returns
1171 -------
1172 an : ndarray
1173 1-D array of weights to apply to the function at the provided sample
1174 positions.
1175 B : float
1176 Error coefficient.
1178 Notes
1179 -----
1180 Normally, the Newton-Cotes rules are used on smaller integration
1181 regions and a composite rule is used to return the total integral.
1183 Examples
1184 --------
1185 Compute the integral of sin(x) in [0, :math:`\pi`]:
1187 >>> from scipy.integrate import newton_cotes
1188 >>> import numpy as np
1189 >>> def f(x):
1190 ... return np.sin(x)
1191 >>> a = 0
1192 >>> b = np.pi
1193 >>> exact = 2
1194 >>> for N in [2, 4, 6, 8, 10]:
1195 ... x = np.linspace(a, b, N + 1)
1196 ... an, B = newton_cotes(N, 1)
1197 ... dx = (b - a) / N
1198 ... quad = dx * np.sum(an * f(x))
1199 ... error = abs(quad - exact)
1200 ... print('{:2d} {:10.9f} {:.5e}'.format(N, quad, error))
1201 ...
1202 2 2.094395102 9.43951e-02
1203 4 1.998570732 1.42927e-03
1204 6 2.000017814 1.78136e-05
1205 8 1.999999835 1.64725e-07
1206 10 2.000000001 1.14677e-09
1208 """
1209 try:
1210 N = len(rn)-1
1211 if equal:
1212 rn = np.arange(N+1)
1213 elif np.all(np.diff(rn) == 1):
1214 equal = 1
1215 except Exception:
1216 N = rn
1217 rn = np.arange(N+1)
1218 equal = 1
1220 if equal and N in _builtincoeffs:
1221 na, da, vi, nb, db = _builtincoeffs[N]
1222 an = na * np.array(vi, dtype=float) / da
1223 return an, float(nb)/db
1225 if (rn[0] != 0) or (rn[-1] != N):
1226 raise ValueError("The sample positions must start at 0"
1227 " and end at N")
1228 yi = rn / float(N)
1229 ti = 2 * yi - 1
1230 nvec = np.arange(N+1)
1231 C = ti ** nvec[:, np.newaxis]
1232 Cinv = np.linalg.inv(C)
1233 # improve precision of result
1234 for i in range(2):
1235 Cinv = 2*Cinv - Cinv.dot(C).dot(Cinv)
1236 vec = 2.0 / (nvec[::2]+1)
1237 ai = Cinv[:, ::2].dot(vec) * (N / 2.)
1239 if (N % 2 == 0) and equal:
1240 BN = N/(N+3.)
1241 power = N+2
1242 else:
1243 BN = N/(N+2.)
1244 power = N+1
1246 BN = BN - np.dot(yi**power, ai)
1247 p1 = power+1
1248 fac = power*math.log(N) - gammaln(p1)
1249 fac = math.exp(fac)
1250 return ai, BN*fac
1253def _qmc_quad_iv(func, a, b, n_points, n_estimates, qrng, log):
1255 # lazy import to avoid issues with partially-initialized submodule
1256 if not hasattr(qmc_quad, 'qmc'):
1257 from scipy import stats
1258 qmc_quad.stats = stats
1259 else:
1260 stats = qmc_quad.stats
1262 if not callable(func):
1263 message = "`func` must be callable."
1264 raise TypeError(message)
1266 # a, b will be modified, so copy. Oh well if it's copied twice.
1267 a = np.atleast_1d(a).copy()
1268 b = np.atleast_1d(b).copy()
1269 a, b = np.broadcast_arrays(a, b)
1270 dim = a.shape[0]
1272 try:
1273 func((a + b) / 2)
1274 except Exception as e:
1275 message = ("`func` must evaluate the integrand at points within "
1276 "the integration range; e.g. `func( (a + b) / 2)` "
1277 "must return the integrand at the centroid of the "
1278 "integration volume.")
1279 raise ValueError(message) from e
1281 try:
1282 func(np.array([a, b]).T)
1283 vfunc = func
1284 except Exception as e:
1285 message = ("Exception encountered when attempting vectorized call to "
1286 f"`func`: {e}. For better performance, `func` should "
1287 "accept two-dimensional array `x` with shape `(len(a), "
1288 "n_points)` and return an array of the integrand value at "
1289 "each of the `n_points.")
1290 warnings.warn(message, stacklevel=3)
1292 def vfunc(x):
1293 return np.apply_along_axis(func, axis=-1, arr=x)
1295 n_points_int = np.int64(n_points)
1296 if n_points != n_points_int:
1297 message = "`n_points` must be an integer."
1298 raise TypeError(message)
1300 n_estimates_int = np.int64(n_estimates)
1301 if n_estimates != n_estimates_int:
1302 message = "`n_estimates` must be an integer."
1303 raise TypeError(message)
1305 if qrng is None:
1306 qrng = stats.qmc.Halton(dim)
1307 elif not isinstance(qrng, stats.qmc.QMCEngine):
1308 message = "`qrng` must be an instance of scipy.stats.qmc.QMCEngine."
1309 raise TypeError(message)
1311 if qrng.d != a.shape[0]:
1312 message = ("`qrng` must be initialized with dimensionality equal to "
1313 "the number of variables in `a`, i.e., "
1314 "`qrng.random().shape[-1]` must equal `a.shape[0]`.")
1315 raise ValueError(message)
1317 rng_seed = getattr(qrng, 'rng_seed', None)
1318 rng = stats._qmc.check_random_state(rng_seed)
1320 if log not in {True, False}:
1321 message = "`log` must be boolean (`True` or `False`)."
1322 raise TypeError(message)
1324 return (vfunc, a, b, n_points_int, n_estimates_int, qrng, rng, log, stats)
1327QMCQuadResult = namedtuple('QMCQuadResult', ['integral', 'standard_error'])
1330def qmc_quad(func, a, b, *, n_estimates=8, n_points=1024, qrng=None,
1331 log=False):
1332 """
1333 Compute an integral in N-dimensions using Quasi-Monte Carlo quadrature.
1335 Parameters
1336 ----------
1337 func : callable
1338 The integrand. Must accept a single argument ``x``, an array which
1339 specifies the point(s) at which to evaluate the scalar-valued
1340 integrand, and return the value(s) of the integrand.
1341 For efficiency, the function should be vectorized to accept an array of
1342 shape ``(d, n_points)``, where ``d`` is the number of variables (i.e.
1343 the dimensionality of the function domain) and `n_points` is the number
1344 of quadrature points, and return an array of shape ``(n_points,)``,
1345 the integrand at each quadrature point.
1346 a, b : array-like
1347 One-dimensional arrays specifying the lower and upper integration
1348 limits, respectively, of each of the ``d`` variables.
1349 n_estimates, n_points : int, optional
1350 `n_estimates` (default: 8) statistically independent QMC samples, each
1351 of `n_points` (default: 1024) points, will be generated by `qrng`.
1352 The total number of points at which the integrand `func` will be
1353 evaluated is ``n_points * n_estimates``. See Notes for details.
1354 qrng : `~scipy.stats.qmc.QMCEngine`, optional
1355 An instance of the QMCEngine from which to sample QMC points.
1356 The QMCEngine must be initialized to a number of dimensions ``d``
1357 corresponding with the number of variables ``x1, ..., xd`` passed to
1358 `func`.
1359 The provided QMCEngine is used to produce the first integral estimate.
1360 If `n_estimates` is greater than one, additional QMCEngines are
1361 spawned from the first (with scrambling enabled, if it is an option.)
1362 If a QMCEngine is not provided, the default `scipy.stats.qmc.Halton`
1363 will be initialized with the number of dimensions determine from
1364 the length of `a`.
1365 log : boolean, default: False
1366 When set to True, `func` returns the log of the integrand, and
1367 the result object contains the log of the integral.
1369 Returns
1370 -------
1371 result : object
1372 A result object with attributes:
1374 integral : float
1375 The estimate of the integral.
1376 standard_error :
1377 The error estimate. See Notes for interpretation.
1379 Notes
1380 -----
1381 Values of the integrand at each of the `n_points` points of a QMC sample
1382 are used to produce an estimate of the integral. This estimate is drawn
1383 from a population of possible estimates of the integral, the value of
1384 which we obtain depends on the particular points at which the integral
1385 was evaluated. We perform this process `n_estimates` times, each time
1386 evaluating the integrand at different scrambled QMC points, effectively
1387 drawing i.i.d. random samples from the population of integral estimates.
1388 The sample mean :math:`m` of these integral estimates is an
1389 unbiased estimator of the true value of the integral, and the standard
1390 error of the mean :math:`s` of these estimates may be used to generate
1391 confidence intervals using the t distribution with ``n_estimates - 1``
1392 degrees of freedom. Perhaps counter-intuitively, increasing `n_points`
1393 while keeping the total number of function evaluation points
1394 ``n_points * n_estimates`` fixed tends to reduce the actual error, whereas
1395 increasing `n_estimates` tends to decrease the error estimate.
1397 Examples
1398 --------
1399 QMC quadrature is particularly useful for computing integrals in higher
1400 dimensions. An example integrand is the probability density function
1401 of a multivariate normal distribution.
1403 >>> import numpy as np
1404 >>> from scipy import stats
1405 >>> dim = 8
1406 >>> mean = np.zeros(dim)
1407 >>> cov = np.eye(dim)
1408 >>> def func(x):
1409 ... # `multivariate_normal` expects the _last_ axis to correspond with
1410 ... # the dimensionality of the space, so `x` must be transposed
1411 ... return stats.multivariate_normal.pdf(x.T, mean, cov)
1413 To compute the integral over the unit hypercube:
1415 >>> from scipy.integrate import qmc_quad
1416 >>> a = np.zeros(dim)
1417 >>> b = np.ones(dim)
1418 >>> rng = np.random.default_rng()
1419 >>> qrng = stats.qmc.Halton(d=dim, seed=rng)
1420 >>> n_estimates = 8
1421 >>> res = qmc_quad(func, a, b, n_estimates=n_estimates, qrng=qrng)
1422 >>> res.integral, res.standard_error
1423 (0.00018429555666024108, 1.0389431116001344e-07)
1425 A two-sided, 99% confidence interval for the integral may be estimated
1426 as:
1428 >>> t = stats.t(df=n_estimates-1, loc=res.integral,
1429 ... scale=res.standard_error)
1430 >>> t.interval(0.99)
1431 (0.0001839319802536469, 0.00018465913306683527)
1433 Indeed, the value reported by `scipy.stats.multivariate_normal` is
1434 within this range.
1436 >>> stats.multivariate_normal.cdf(b, mean, cov, lower_limit=a)
1437 0.00018430867675187443
1439 """
1440 args = _qmc_quad_iv(func, a, b, n_points, n_estimates, qrng, log)
1441 func, a, b, n_points, n_estimates, qrng, rng, log, stats = args
1443 def sum_product(integrands, dA, log=False):
1444 if log:
1445 return logsumexp(integrands) + np.log(dA)
1446 else:
1447 return np.sum(integrands * dA)
1449 def mean(estimates, log=False):
1450 if log:
1451 return logsumexp(estimates) - np.log(n_estimates)
1452 else:
1453 return np.mean(estimates)
1455 def std(estimates, m=None, ddof=0, log=False):
1456 m = m or mean(estimates, log)
1457 if log:
1458 estimates, m = np.broadcast_arrays(estimates, m)
1459 temp = np.vstack((estimates, m + np.pi * 1j))
1460 diff = logsumexp(temp, axis=0)
1461 return np.real(0.5 * (logsumexp(2 * diff)
1462 - np.log(n_estimates - ddof)))
1463 else:
1464 return np.std(estimates, ddof=ddof)
1466 def sem(estimates, m=None, s=None, log=False):
1467 m = m or mean(estimates, log)
1468 s = s or std(estimates, m, ddof=1, log=log)
1469 if log:
1470 return s - 0.5*np.log(n_estimates)
1471 else:
1472 return s / np.sqrt(n_estimates)
1474 # The sign of the integral depends on the order of the limits. Fix this by
1475 # ensuring that lower bounds are indeed lower and setting sign of resulting
1476 # integral manually
1477 if np.any(a == b):
1478 message = ("A lower limit was equal to an upper limit, so the value "
1479 "of the integral is zero by definition.")
1480 warnings.warn(message, stacklevel=2)
1481 return QMCQuadResult(-np.inf if log else 0, 0)
1483 i_swap = b < a
1484 sign = (-1)**(i_swap.sum(axis=-1)) # odd # of swaps -> negative
1485 a[i_swap], b[i_swap] = b[i_swap], a[i_swap]
1487 A = np.prod(b - a)
1488 dA = A / n_points
1490 estimates = np.zeros(n_estimates)
1491 rngs = _rng_spawn(qrng.rng, n_estimates)
1492 for i in range(n_estimates):
1493 # Generate integral estimate
1494 sample = qrng.random(n_points)
1495 # The rationale for transposing is that this allows users to easily
1496 # unpack `x` into separate variables, if desired. This is consistent
1497 # with the `xx` array passed into the `scipy.integrate.nquad` `func`.
1498 x = stats.qmc.scale(sample, a, b).T # (n_dim, n_points)
1499 integrands = func(x)
1500 estimates[i] = sum_product(integrands, dA, log)
1502 # Get a new, independently-scrambled QRNG for next time
1503 qrng = type(qrng)(seed=rngs[i], **qrng._init_quad)
1505 integral = mean(estimates, log)
1506 standard_error = sem(estimates, m=integral, log=log)
1507 integral = integral + np.pi*1j if (log and sign < 0) else integral*sign
1508 return QMCQuadResult(integral, standard_error)