Coverage for /usr/lib/python3/dist-packages/scipy/stats/qmc.py: 100%
1 statements
« prev ^ index » next coverage.py v7.9.1, created at 2025-06-14 15:55 +0200
« prev ^ index » next coverage.py v7.9.1, created at 2025-06-14 15:55 +0200
1r"""
2====================================================
3Quasi-Monte Carlo submodule (:mod:`scipy.stats.qmc`)
4====================================================
6.. currentmodule:: scipy.stats.qmc
8This module provides Quasi-Monte Carlo generators and associated helper
9functions.
12Quasi-Monte Carlo
13=================
15Engines
16-------
18.. autosummary::
19 :toctree: generated/
21 QMCEngine
22 Sobol
23 Halton
24 LatinHypercube
25 PoissonDisk
26 MultinomialQMC
27 MultivariateNormalQMC
29Helpers
30-------
32.. autosummary::
33 :toctree: generated/
35 discrepancy
36 update_discrepancy
37 scale
40Introduction to Quasi-Monte Carlo
41=================================
43Quasi-Monte Carlo (QMC) methods [1]_, [2]_, [3]_ provide an
44:math:`n \times d` array of numbers in :math:`[0,1]`. They can be used in
45place of :math:`n` points from the :math:`U[0,1]^{d}` distribution. Compared to
46random points, QMC points are designed to have fewer gaps and clumps. This is
47quantified by discrepancy measures [4]_. From the Koksma-Hlawka
48inequality [5]_ we know that low discrepancy reduces a bound on
49integration error. Averaging a function :math:`f` over :math:`n` QMC points
50can achieve an integration error close to :math:`O(n^{-1})` for well
51behaved functions [2]_.
53Most QMC constructions are designed for special values of :math:`n`
54such as powers of 2 or large primes. Changing the sample
55size by even one can degrade their performance, even their
56rate of convergence [6]_. For instance :math:`n=100` points may give less
57accuracy than :math:`n=64` if the method was designed for :math:`n=2^m`.
59Some QMC constructions are extensible in :math:`n`: we can find
60another special sample size :math:`n' > n` and often an infinite
61sequence of increasing special sample sizes. Some QMC
62constructions are extensible in :math:`d`: we can increase the dimension,
63possibly to some upper bound, and typically without requiring
64special values of :math:`d`. Some QMC methods are extensible in
65both :math:`n` and :math:`d`.
67QMC points are deterministic. That makes it hard to estimate the accuracy of
68integrals estimated by averages over QMC points. Randomized QMC (RQMC) [7]_
69points are constructed so that each point is individually :math:`U[0,1]^{d}`
70while collectively the :math:`n` points retain their low discrepancy.
71One can make :math:`R` independent replications of RQMC points to
72see how stable a computation is. From :math:`R` independent values,
73a t-test (or bootstrap t-test [8]_) then gives approximate confidence
74intervals on the mean value. Some RQMC methods produce a
75root mean squared error that is actually :math:`o(1/n)` and smaller than
76the rate seen in unrandomized QMC. An intuitive explanation is
77that the error is a sum of many small ones and random errors
78cancel in a way that deterministic ones do not. RQMC also
79has advantages on integrands that are singular or, for other
80reasons, fail to be Riemann integrable.
82(R)QMC cannot beat Bahkvalov's curse of dimension (see [9]_). For
83any random or deterministic method, there are worst case functions
84that will give it poor performance in high dimensions. A worst
85case function for QMC might be 0 at all n points but very
86large elsewhere. Worst case analyses get very pessimistic
87in high dimensions. (R)QMC can bring a great improvement over
88MC when the functions on which it is used are not worst case.
89For instance (R)QMC can be especially effective on integrands
90that are well approximated by sums of functions of
91some small number of their input variables at a time [10]_, [11]_.
92That property is often a surprising finding about those functions.
94Also, to see an improvement over IID MC, (R)QMC requires a bit of smoothness of
95the integrand, roughly the mixed first order derivative in each direction,
96:math:`\partial^d f/\partial x_1 \cdots \partial x_d`, must be integral.
97For instance, a function that is 1 inside the hypersphere and 0 outside of it
98has infinite variation in the sense of Hardy and Krause for any dimension
99:math:`d = 2`.
101Scrambled nets are a kind of RQMC that have some valuable robustness
102properties [12]_. If the integrand is square integrable, they give variance
103:math:`var_{SNET} = o(1/n)`. There is a finite upper bound on
104:math:`var_{SNET} / var_{MC}` that holds simultaneously for every square
105integrable integrand. Scrambled nets satisfy a strong law of large numbers
106for :math:`f` in :math:`L^p` when :math:`p>1`. In some
107special cases there is a central limit theorem [13]_. For smooth enough
108integrands they can achieve RMSE nearly :math:`O(n^{-3})`. See [12]_
109for references about these properties.
111The main kinds of QMC methods are lattice rules [14]_ and digital
112nets and sequences [2]_, [15]_. The theories meet up in polynomial
113lattice rules [16]_ which can produce digital nets. Lattice rules
114require some form of search for good constructions. For digital
115nets there are widely used default constructions.
117The most widely used QMC methods are Sobol' sequences [17]_.
118These are digital nets. They are extensible in both :math:`n` and :math:`d`.
119They can be scrambled. The special sample sizes are powers
120of 2. Another popular method are Halton sequences [18]_.
121The constructions resemble those of digital nets. The earlier
122dimensions have much better equidistribution properties than
123later ones. There are essentially no special sample sizes.
124They are not thought to be as accurate as Sobol' sequences.
125They can be scrambled. The nets of Faure [19]_ are also widely
126used. All dimensions are equally good, but the special sample
127sizes grow rapidly with dimension :math:`d`. They can be scrambled.
128The nets of Niederreiter and Xing [20]_ have the best asymptotic
129properties but have not shown good empirical performance [21]_.
131Higher order digital nets are formed by a digit interleaving process
132in the digits of the constructed points. They can achieve higher
133levels of asymptotic accuracy given higher smoothness conditions on :math:`f`
134and they can be scrambled [22]_. There is little or no empirical work
135showing the improved rate to be attained.
137Using QMC is like using the entire period of a small random
138number generator. The constructions are similar and so
139therefore are the computational costs [23]_.
141(R)QMC is sometimes improved by passing the points through
142a baker's transformation (tent function) prior to using them.
143That function has the form :math:`1-2|x-1/2|`. As :math:`x` goes from 0 to
1441, this function goes from 0 to 1 and then back. It is very
145useful to produce a periodic function for lattice rules [14]_,
146and sometimes it improves the convergence rate [24]_.
148It is not straightforward to apply QMC methods to Markov
149chain Monte Carlo (MCMC). We can think of MCMC as using
150:math:`n=1` point in :math:`[0,1]^{d}` for very large :math:`d`, with
151ergodic results corresponding to :math:`d \to \infty`. One proposal is
152in [25]_ and under strong conditions an improved rate of convergence
153has been shown [26]_.
155Returning to Sobol' points: there are many versions depending
156on what are called direction numbers. Those are the result of
157searches and are tabulated. A very widely used set of direction
158numbers come from [27]_. It is extensible in dimension up to
159:math:`d=21201`.
161References
162----------
163.. [1] Owen, Art B. "Monte Carlo Book: the Quasi-Monte Carlo parts." 2019.
164.. [2] Niederreiter, Harald. "Random number generation and quasi-Monte Carlo
165 methods." Society for Industrial and Applied Mathematics, 1992.
166.. [3] Dick, Josef, Frances Y. Kuo, and Ian H. Sloan. "High-dimensional
167 integration: the quasi-Monte Carlo way." Acta Numerica no. 22: 133, 2013.
168.. [4] Aho, A. V., C. Aistleitner, T. Anderson, K. Appel, V. Arnol'd, N.
169 Aronszajn, D. Asotsky et al. "W. Chen et al.(eds.), "A Panorama of
170 Discrepancy Theory", Sringer International Publishing,
171 Switzerland: 679, 2014.
172.. [5] Hickernell, Fred J. "Koksma-Hlawka Inequality." Wiley StatsRef:
173 Statistics Reference Online, 2014.
174.. [6] Owen, Art B. "On dropping the first Sobol' point." :arxiv:`2008.08051`,
175 2020.
176.. [7] L'Ecuyer, Pierre, and Christiane Lemieux. "Recent advances in randomized
177 quasi-Monte Carlo methods." In Modeling uncertainty, pp. 419-474. Springer,
178 New York, NY, 2002.
179.. [8] DiCiccio, Thomas J., and Bradley Efron. "Bootstrap confidence
180 intervals." Statistical science: 189-212, 1996.
181.. [9] Dimov, Ivan T. "Monte Carlo methods for applied scientists." World
182 Scientific, 2008.
183.. [10] Caflisch, Russel E., William J. Morokoff, and Art B. Owen. "Valuation
184 of mortgage backed securities using Brownian bridges to reduce effective
185 dimension." Journal of Computational Finance: no. 1 27-46, 1997.
186.. [11] Sloan, Ian H., and Henryk Wozniakowski. "When are quasi-Monte Carlo
187 algorithms efficient for high dimensional integrals?." Journal of Complexity
188 14, no. 1 (1998): 1-33.
189.. [12] Owen, Art B., and Daniel Rudolf, "A strong law of large numbers for
190 scrambled net integration." SIAM Review, to appear.
191.. [13] Loh, Wei-Liem. "On the asymptotic distribution of scrambled net
192 quadrature." The Annals of Statistics 31, no. 4: 1282-1324, 2003.
193.. [14] Sloan, Ian H. and S. Joe. "Lattice methods for multiple integration."
194 Oxford University Press, 1994.
195.. [15] Dick, Josef, and Friedrich Pillichshammer. "Digital nets and sequences:
196 discrepancy theory and quasi-Monte Carlo integration." Cambridge University
197 Press, 2010.
198.. [16] Dick, Josef, F. Kuo, Friedrich Pillichshammer, and I. Sloan.
199 "Construction algorithms for polynomial lattice rules for multivariate
200 integration." Mathematics of computation 74, no. 252: 1895-1921, 2005.
201.. [17] Sobol', Il'ya Meerovich. "On the distribution of points in a cube and
202 the approximate evaluation of integrals." Zhurnal Vychislitel'noi Matematiki
203 i Matematicheskoi Fiziki 7, no. 4: 784-802, 1967.
204.. [18] Halton, John H. "On the efficiency of certain quasi-random sequences of
205 points in evaluating multi-dimensional integrals." Numerische Mathematik 2,
206 no. 1: 84-90, 1960.
207.. [19] Faure, Henri. "Discrepance de suites associees a un systeme de
208 numeration (en dimension s)." Acta arithmetica 41, no. 4: 337-351, 1982.
209.. [20] Niederreiter, Harold, and Chaoping Xing. "Low-discrepancy sequences and
210 global function fields with many rational places." Finite Fields and their
211 applications 2, no. 3: 241-273, 1996.
212.. [21] Hong, Hee Sun, and Fred J. Hickernell. "Algorithm 823: Implementing
213 scrambled digital sequences." ACM Transactions on Mathematical Software
214 (TOMS) 29, no. 2: 95-109, 2003.
215.. [22] Dick, Josef. "Higher order scrambled digital nets achieve the optimal
216 rate of the root mean square error for smooth integrands." The Annals of
217 Statistics 39, no. 3: 1372-1398, 2011.
218.. [23] Niederreiter, Harald. "Multidimensional numerical integration using
219 pseudorandom numbers." In Stochastic Programming 84 Part I, pp. 17-38.
220 Springer, Berlin, Heidelberg, 1986.
221.. [24] Hickernell, Fred J. "Obtaining O (N-2+e) Convergence for Lattice
222 Quadrature Rules." In Monte Carlo and Quasi-Monte Carlo Methods 2000,
223 pp. 274-289. Springer, Berlin, Heidelberg, 2002.
224.. [25] Owen, Art B., and Seth D. Tribble. "A quasi-Monte Carlo Metropolis
225 algorithm." Proceedings of the National Academy of Sciences 102,
226 no. 25: 8844-8849, 2005.
227.. [26] Chen, Su. "Consistency and convergence rate of Markov chain quasi Monte
228 Carlo with examples." PhD diss., Stanford University, 2011.
229.. [27] Joe, Stephen, and Frances Y. Kuo. "Constructing Sobol sequences with
230 better two-dimensional projections." SIAM Journal on Scientific Computing
231 30, no. 5: 2635-2654, 2008.
233"""
234from ._qmc import *