that can be generated, which is ``2**bits``. Maximal value is 64.
    It does not correspond to the return type, which is always
    ``np.float64`` to prevent points from repeating themselves.
    Default is None, which for backward compatibility, corresponds to 30.

    .. versionadded:: 1.9.0
optimization : {None, "random-cd", "lloyd"}, optional
    Whether to use an optimization scheme to improve the quality after
    sampling. Note that this is a post-processing step that does not
    guarantee that all properties of the sample will be conserved.
    Default is None.

    * ``random-cd``: random permutations of coordinates to lower the
      centered discrepancy. The best sample based on the centered
      discrepancy is constantly updated. Centered discrepancy-based
      sampling shows better space-filling robustness toward 2D and 3D
      subprojections compared to using other discrepancy measures.
    * ``lloyd``: Perturb samples using a modified Lloyd-Max algorithm.
      The process converges to equally spaced samples.

    .. versionadded:: 1.10.0

rng : `numpy.random.Generator`, optional
    Pseudorandom number generator state. When `rng` is None, a new
    `numpy.random.Generator` is created using entropy from the
    operating system. Types other than `numpy.random.Generator` are
    passed to `numpy.random.default_rng` to instantiate a ``Generator``.

    .. versionchanged:: 1.15.0

        As part of the `SPEC-007 <https://scientific-python.org/specs/spec-0007/>`_
        transition from use of `numpy.random.RandomState` to
        `numpy.random.Generator`, this keyword was changed from `seed` to
        `rng`. For an interim period, both keywords will continue to work, although
        only one may be specified at a time. After the interim period, function
        calls using the `seed` keyword will emit warnings. Following a
        deprecation period, the `seed` keyword will be removed.

Notes
-----
Sobol' sequences [1]_ provide :math:`n=2^m` low discrepancy points in
:math:`[0,1)^{d}`. Scrambling them [3]_ makes them suitable for singular
integrands, provides a means of error estimation, and can improve their
rate of convergence. The scrambling strategy which is implemented is a
(left) linear matrix scramble (LMS) followed by a digital random shift
(LMS+shift) [2]_.

There are many versions of Sobol' sequences depending on their
'direction numbers'. This code uses direction numbers from [4]_. 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 [5]_.

   If :math:`n=2^m` points are not enough then one should take :math:`2^M`
   points for :math:`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 :math:`B` of bits.
   After :math:`2^B` points have been generated, the sequence would
   repeat. Hence, an error is raised.
   The number of bits can be controlled with the parameter `bits`.

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] J. Matousek, "On the L2-discrepancy for anchored boxes."
   J. of Complexity 14, 527-556, 1998.
.. [3] Art B. Owen, "Scrambling Sobol and Niederreiter-Xing points."
   Journal of Complexity, 14(4):466-489, December 1998.
.. [4] S. Joe and F. Y. Kuo, "Constructing sobol sequences with better
   two-dimensional projections." SIAM Journal on Scientific Computing,
   30(5):2635-2654, 2008.
.. [5] 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]])

Ú