scipy.stats.qmc.Sobol¶
- class scipy.stats.qmc.Sobol(d, *, scramble=True, seed=None)[source]¶
Engine for generating (scrambled) Sobol’ sequences.
Sobol’ sequences are low-discrepancy, quasi-random numbers. Points can be drawn using two methods:
random_base2
: safely draw \(n=2^m\) points. This method guarantees the balance properties of the sequence.random
: draw an arbitrary number of points from the sequence. See warning below.
- Parameters
- dint
Dimensionality of the sequence. Max dimensionality is 21201.
- scramblebool, optional
If True, use Owen scrambling. Otherwise no scrambling is done. Default is True.
- seed{None, int,
numpy.random.Generator
}, optional If seed is None the
numpy.random.Generator
singleton is used. If seed is an int, a newGenerator
instance is used, seeded with seed. If seed is already aGenerator
instance then that instance is used.
Notes
Sobol’ sequences [1] provide \(n=2^m\) low discrepancy points in \([0,1)^{d}\). Scrambling them [2] makes them suitable for singular integrands, provides a means of error estimation, and can improve their rate of convergence.
There are many versions of Sobol’ sequences depending on their ‘direction numbers’. This code uses direction numbers from [3]. Hence, the maximum number of dimension is 21201. The direction numbers have been precomputed with search criterion 6 and can be retrieved at https://web.maths.unsw.edu.au/~fkuo/sobol/.
Warning
Sobol’ sequences are a quadrature rule and they lose their balance properties if one uses a sample size that is not a power of 2, or skips the first point, or thins the sequence [4].
If \(n=2^m\) points are not enough then one should take \(2^M\) points for \(M>m\). When scrambling, the number R of independent replicates does not have to be a power of 2.
Sobol’ sequences are generated to some number \(B\) of bits. After \(2^B\) points have been generated, the sequence will repeat. Currently \(B=30\).
References
- 1
I. M. Sobol. The distribution of points in a cube and the accurate evaluation of integrals. Zh. Vychisl. Mat. i Mat. Phys., 7:784-802, 1967.
- 2
Art B. Owen. Scrambling Sobol and Niederreiter-Xing points. Journal of Complexity, 14(4):466-489, December 1998.
- 3
S. Joe and F. Y. Kuo. Constructing sobol sequences with better two-dimensional projections. SIAM Journal on Scientific Computing, 30(5):2635-2654, 2008.
- 4
Art B. Owen. On dropping the first Sobol’ point. arXiv 2008.08051, 2020.
Examples
Generate samples from a low discrepancy sequence of Sobol’.
>>> from scipy.stats import qmc >>> sampler = qmc.Sobol(d=2, scramble=False) >>> sample = sampler.random_base2(m=3) >>> sample array([[0. , 0. ], [0.5 , 0.5 ], [0.75 , 0.25 ], [0.25 , 0.75 ], [0.375, 0.375], [0.875, 0.875], [0.625, 0.125], [0.125, 0.625]])
Compute the quality of the sample using the discrepancy criterion.
>>> qmc.discrepancy(sample) 0.013882107204860938
To continue an existing design, extra points can be obtained by calling again
random_base2
. Alternatively, you can skip some points like:>>> _ = sampler.reset() >>> _ = sampler.fast_forward(4) >>> sample_continued = sampler.random_base2(m=2) >>> sample_continued array([[0.375, 0.375], [0.875, 0.875], [0.625, 0.125], [0.125, 0.625]])
Finally, samples can be scaled to bounds.
>>> l_bounds = [0, 2] >>> u_bounds = [10, 5] >>> qmc.scale(sample_continued, l_bounds, u_bounds) array([[3.75 , 3.125], [8.75 , 4.625], [6.25 , 2.375], [1.25 , 3.875]])
Methods
fast_forward
(n)Fast-forward the sequence by n positions.
random
([n])Draw next point(s) in the Sobol' sequence.
random_base2
(m)Draw point(s) from the Sobol' sequence.
reset
()Reset the engine to base state.