Repository: scipy/scipy
Pull request: ENH: add stats.qmc module with quasi Monte Carlo functionality
Add `scipy.stats.qmc`.

#### Reference issue
Closes https://github.com/scipy/scipy/issues/9695
#### What does this implement/fix?
Provide a set of functions to create and assess quasi-Monte Carlo Design of Experiments.

It provides a generic class `scipy.stats.qmc.QMCEngine` which defines a QMC engine/sampler. An engine is state aware: it can be continued, advanced and reseted. 4 base samplers are available:

* `scipy.stats.qmc.Sobol` the well known Sobol' low discrepancy sequence. Several warnings have been added to guide the user into properly using this sampler. The sequence is scrambled by default.
* `scipy.stats.qmc.Halton`: Halton low discrepancy sequence. The sequence is scrambled by default.
* `scipy.stats.qmc.LatinHypercube`: plain LHS design.
* `scipy.stats.qmc.OrthogonalLatinHypercube`: Orthogonal version of LHS which is more uniform than plain LHS.

And 2 special samplers are available:

* `scipy.stats.qmc.MultinomialQMC`: sampling from a multinomial distribution using any of the base `scipy.stats.qmc.QMCEngine`.
* `scipy.stats.qmc.MultivariateNormalQMC`: sampling from a multivariate Normal using any of the base `scipy.stats.qmc.QMCEngine`.

The module also provide the following helpers:

* `scipy.stats.qmc.discrepancy`: assess the quality of a set of points in terms of space coverage.
* `scipy.stats.qmc.update_discrepancy`: can be used in an optimization loop to construct a good set of points.
* `scipy.stats.qmc.scale`: easily scale a set of points from (to) the unit interval to (from) a given range.

The module quality has been assessed by specialists from the domain (see the long discussion bellow for more details). Following are some convergence plots  of Sobol' and Halton.

* gist for Halton: https://gist.github.com/tupui/ee6e685219c06c429338569e7ee7ff2a
* gist for Sobol': https://gist.github.com/tupui/fb6e219b1dd2316b7498ebce231bfff5

![halton_art2](https://user-images.githubusercontent.com/23188539/106391045-a4e70500-63eb-11eb-92d7-78c588f7e383.png)
![sobol_typeB](https://user-images.githubusercontent.com/23188539/106391046-a6183200-63eb-11eb-8f15-57037115ead1.png)




Problem statement:
[Question] Desire for qMC library in scipy?
I have a fast Cython implementation of a Sobol low-discrepancy quasi-random number generator using Owen scrambling.  On top of that I have some transformations (inverse transform, Box-Muller) that make it easy to draw quasi-MC samples from, say, a multivariate normal distribution. This code has been used successfully in an industry setting.

Before I go through the effort of creating a PR, I wanted to see whether this is something that people would like to see in scipy. 

The api could look something like the following:
```
from scipy.qmc import Sobol

sobol = Sobol(dimen=d, seed=1234)
samples = sobol.rvs(n)
```

```
from scipy.qmc import MultivariateNormalQMC

mvn = MultivariateNormalQMC(mean=np.zeros(2), cov=np.array([[1, 0.1], [0.1, 1]]))
samples = mvn.rvs(n)
```

Hints from issue discussion:
Yes!
Plus there is already Sobol' QMC in `from scipy.optimize._shgo_lib import sobol_seq`. There is some current work on statsmodels to implement some sampling methods: https://github.com/statsmodels/statsmodels/pull/4104 and WIP: https://github.com/statsmodels/statsmodels/pull/4113
There definitely is enough interest, a PR would be very welcome @Balandat.

> `from scipy.qmc import MultivariateNormalQMC`

A separate top-level module does not seem warranted. This seems to fit well in `scipy.stats`. Depending on how many functions you want to add (sounds like quite a few?), we could create `scipy.stats.qmc`.
Sorry I have been swamped with stuff, hopefully I can get to this soon.

FWIW, we upstreamed a modified implementation into pytorch: https://pytorch.org/docs/stable/torch.html#quasi-random-sampling
@rgommers Is there any interest for other QMC like Halton (which is slightly better in dim <~ 10)? I put a version in statsmodels. I could also propose some LHS. 
> @rgommers Is there any interest for other QMC like Halton (which is slightly better in dim <~ 10)? I put a version in statsmodels. I could also propose some LHS.

In principle yes I think, but we try hard to not duplicate things between `statsmodels` and `scipy.stats`. The Halton function and `statsmodels.tools.sequences` seem to not have made it into the Statsmodels docs though (and your WIP PRs actually moves it to `doe.sequences`), so it's unclear to me if that's public or not - my impression is it's private still.

If we're going to have multiple, either in SciPy or in Statsmodels, it would be good if they have a uniform API as much as possible. The `halton`, `van_der_corput` and `latin_hypercube` signatures are all slightly different and that doesn't look necessary.

I think a self-contained module would be better rather than putting them into a "design of experiments" (sub)module. Other than that, no strong preference between SciPy and Statsmodels.

Cc @josef-pkt, @bashtage  


IMO both of these are better fits for scipy than statsmodels.  I think those functions can be treated as private, and so could be quietly and quickly deprecated (with a pointer to scipy.MOD) if SciPy goes forward with these. 
FWIW there is a sobol squence generator wrapped in Python with support for 21201 dimensions in [pyscenarios](https://github.com/crusaderky/pyscenarios).  Uses numba and so is fast. 
As for a design, it seems like it would make sense to split the transformations from the generation. Transformation is a generic problem that can take any form fo quasi-random-like generates. 
@bashtage I agree with you that these methods are more general than "just" statistically relevant. For lots of optimization/integration problems we want to have better number generators. But I am not sure about keeping these private. I would want to have an alternative to `np.random.random()` which is accessible.

Commit message: Refactor test setup and some raises

Patch to review:
diff --git a/.gitignore b/.gitignore
index 03242cc661dc..45d377db1ce8 100644
--- a/.gitignore
+++ b/.gitignore
@@ -270,6 +270,7 @@ scipy/stats/mvnmodule.c
 scipy/stats/statlibmodule.c
 scipy/stats/vonmises_cython.c
 scipy/stats/_stats.c
+scipy/stats/_sobol.c
 scipy/version.py
 scipy/special/_exprel.c
 scipy/optimize/_group_columns.c
diff --git a/doc/source/stats.qmc.rst b/doc/source/stats.qmc.rst
new file mode 100644
index 000000000000..cb9a544e1e01
--- /dev/null
+++ b/doc/source/stats.qmc.rst
@@ -0,0 +1,4 @@
+.. automodule:: scipy.stats.qmc
+   :no-members:
+   :no-inherited-members:
+   :no-special-members:
diff --git a/doc/source/stats.rst b/doc/source/stats.rst
index 1a4259a679d1..c92dba4a3fda 100644
--- a/doc/source/stats.rst
+++ b/doc/source/stats.rst
@@ -7,3 +7,4 @@
    :hidden:
 
    stats.mstats
+   stats.qmc
diff --git a/scipy/stats/__init__.py b/scipy/stats/__init__.py
index cf769523aef7..6d5646681b9f 100644
--- a/scipy/stats/__init__.py
+++ b/scipy/stats/__init__.py
@@ -373,6 +373,16 @@
 
    gaussian_kde
 
+
+Quasi-Monte Carlo
+=================
+
+.. toctree::
+   :maxdepth: 4
+
+   stats.qmc
+
+
 Warnings used in :mod:`scipy.stats`
 ===================================
 
diff --git a/scipy/stats/_qmc.py b/scipy/stats/_qmc.py
new file mode 100644
index 000000000000..f0075d782af0
--- /dev/null
+++ b/scipy/stats/_qmc.py
@@ -0,0 +1,1528 @@
+"""Quasi-Monte Carlo engines and helpers."""
+from abc import ABC, abstractmethod
+import math
+import warnings
+
+import numpy as np
+
+from scipy.optimize import brute
+from scipy._lib._util import check_random_state
+from scipy.optimize import basinhopping
+from scipy.stats import norm
+from scipy.stats._sobol import (
+    initialize_v, _cscramble, _fill_p_cumulative, _draw, _fast_forward,
+    _categorize, initialize_direction_numbers, _MAXDIM, _MAXBIT
+)
+
+__all__ = ['scale', 'discrepancy', 'QMCEngine', 'Sobol', 'Halton',
+           'OrthogonalLatinHypercube', 'LatinHypercube', 'OptimalDesign',
+           'multinomial_qmc', 'NormalQMC', 'MultivariateNormalQMC']
+
+
+def scale(sample, bounds, reverse=False):
+    r"""Sample scaling from unit hypercube to bounds range.
+
+    To convert a sample from :math:`[0, 1)` to :math:`[a, b), b>a`, the
+    following transformation is used:
+
+    .. math::
+
+        (b - a) \cdot \text{sample} + a
+
+    Parameters
+    ----------
+    sample : array_like (n, d)
+        Sample to scale.
+    bounds : tuple or array_like ([min, d], [max, d])
+        Desired range of transformed data. If `reverse` is True, range of the
+        original data to transform to the unit hypercube.
+    reverse : bool
+        Reverse the transformation from `bounds` to the unit hypercube.
+        Default is False.
+
+    Returns
+    -------
+    sample : array_like (n, d)
+        Scaled sample.
+
+    Examples
+    --------
+    >>> from scipy.stats import qmc
+    >>> bounds = [[-2, 0],
+    ...           [6, 5]]
+    >>> sample = [[0.5 , 0.5 ],
+    ...           [0.75, 0.25]]
+    >>> qmc.scale(sample, bounds)
+    array([[2.  , 2.5 ],
+           [4.  , 1.25]])
+
+    """
+    bounds = np.asarray(bounds)
+    min_ = np.min(bounds, axis=0)
+    max_ = np.max(bounds, axis=0)
+    if not reverse:
+        return sample * (max_ - min_) + min_
+    else:
+        return (sample - min_) / (max_ - min_)
+
+
+def discrepancy(sample, iterative=False, method='CD'):
+    """Discrepancy on a given sample.
+
+    Parameters
+    ----------
+    sample : array_like (n, d)
+        The sample to compute the discrepancy from.
+    iterative : bool
+        Must be False if not using it for updating the discrepancy.
+        Default is False. Refer to the notes for more details.
+    method : str
+        Type of discrepancy, can be ['CD', 'WD', 'MD', 'star']. Refer to
+        the notes for more details. Default is ``CD``.
+
+    Returns
+    -------
+    discrepancy : float
+        Discrepancy.
+
+    Notes
+    -----
+    The discrepancy is a uniformity criterion used to assess the space filling
+    of a number of samples in a hypercube.
+    The discrepancy measures how the spread of the points deviates from a
+    uniform distribution.
+    The lower the value is, the better the coverage of the parameter space is.
+
+    Taking a subspace of the parameter space, we count the number of points in
+    the subpace and compare it to the total number of the points of the
+    sample. This value is substracted by the volume of the subspace. This
+    process is done over all subspaces and the highest value is kept.
+
+    Four methods are available:
+
+    * ``CD``: Centered Discrepancy - subspace involves a corner of the
+      hypercube
+    * ``WD``: Wrap-around Discrepancy - subspace can wrap around bounds
+    * ``MD``: Mixture Discrepancy - mix between CD/WD covering more criteria
+    * ``star``: Star L2-discrepancy - like CD BUT variant to rotation
+
+    Lastly, there using ``iterative=True``, it is possible to compute the
+    discrepancy as if we had :math:`n+1` samples. This is useful if we want
+    to add a point to a sampling and check the candidate which would give the
+    lowest discrepancy. Then you could just update the discrepancy with
+    each candidate using `_update_discrepancy`. This method is faster than
+    computing the discrepancy for a large number of candidates.
+
+    References
+    ----------
+    .. [1] Fang et al. Design and modeling for computer experiments,
+       Computer Science and Data Analysis Series, 2006.
+    .. [2] Zhou Y.-D. et al. Mixture discrepancy for quasi-random point sets
+       Journal of Complexity, 29 (3-4) , pp. 283-301, 2013.
+    .. [3] T. T. Warnock. Computational investigations of low discrepancy point
+       sets, Applications of Number Theory to Numerical
+       Analysis, Academic Press, pp. 319-343, 1972.
+
+    Examples
+    --------
+    Calculate the quality of the sample using the discrepancy:
+
+    >>> from scipy.stats import qmc
+    >>> space = np.array([[1, 3], [2, 6], [3, 2], [4, 5], [5, 1], [6, 4]])
+    >>> bounds = np.array([[0.5, 0.5], [6.5, 6.5]])
+    >>> space = qmc.scale(space, bounds, reverse=True)
+    >>> space
+    array([[0.08333333, 0.41666667],
+           [0.25      , 0.91666667],
+           [0.41666667, 0.25      ],
+           [0.58333333, 0.75      ],
+           [0.75      , 0.08333333],
+           [0.91666667, 0.58333333]])
+    >>> qmc.discrepancy(space)
+    0.008142039609053464
+
+    We can also compute iteratively the discrepancy by using
+    ``iterative=True``.
+
+    >>> disc_init = qmc.discrepancy(space[:-1], iterative=True)
+    >>> disc_init
+    0.04769081147119336
+    >>> _update_discrepancy(space[-1], space[:-1], disc_init) # doctest: +SKIP
+    0.008142039609053513
+
+    """
+    sample = np.asarray(sample)
+
+    n, d = sample.shape
+
+    if iterative:
+        n += 1
+
+    if method == 'CD':
+        abs_ = abs(sample - 0.5)
+        disc1 = np.sum(np.prod(1 + 0.5 * abs_ - 0.5 * abs_ ** 2, axis=1))
+
+        prod_arr = 1
+        for i in range(d):
+            s0 = sample[:, i]
+            prod_arr *= (1 +
+                         0.5 * abs(s0[:, None] - 0.5) + 0.5 * abs(s0 - 0.5) -
+                         0.5 * abs(s0[:, None] - s0))
+        disc2 = prod_arr.sum()
+
+        return ((13.0 / 12.0) ** d - 2.0 / n * disc1 +
+                1.0 / (n ** 2) * disc2)
+    elif method == 'WD':
+        prod_arr = 1
+        for i in range(d):
+            s0 = sample[:, i]
+            x_kikj = abs(s0[:, None] - s0)
+            prod_arr *= 3.0 / 2.0 - x_kikj + x_kikj ** 2
+
+        return - (4.0 / 3.0) ** d + 1.0 / (n ** 2) * prod_arr.sum()
+    elif method == 'MD':
+        abs_ = abs(sample - 0.5)
+        disc1 = np.sum(np.prod(5.0 / 3.0 - 0.25 * abs_ - 0.25 * abs_ ** 2,
+                               axis=1))
+
+        prod_arr = 1
+        for i in range(d):
+            s0 = sample[:, i]
+            prod_arr *= (15.0 / 8.0 -
+                         0.25 * abs(s0[:, None] - 0.5) - 0.25 * abs(s0 - 0.5) -
+                         3.0 / 4.0 * abs(s0[:, None] - s0) +
+                         0.5 * abs(s0[:, None] - s0) ** 2)
+        disc2 = prod_arr.sum()
+
+        disc = (19.0 / 12.0) ** d
+        disc1 = 2.0 / n * disc1
+        disc2 = 1.0 / (n ** 2) * disc2
+
+        return disc - disc1 + disc2
+    elif method == 'star':
+        return np.sqrt(
+            3 ** (-d) - 2 ** (1 - d) / n
+            * np.sum(np.prod(1 - sample ** 2, axis=1))
+            + np.sum([
+                np.prod(1 - np.maximum(sample[k, :], sample[j, :]))
+                for k in range(n) for j in range(n)
+            ]) / n ** 2
+        )
+    else:
+        raise ValueError('{} is not a valid method. Options are '
+                         'CD, WD, MD, star.'.format(method))
+
+
+def _update_discrepancy(x_new, sample, initial_disc):
+    """Update the centered discrepancy with a new sample.
+
+    Parameters
+    ----------
+    x_new : array_like (1, d)
+        The new sample to add in `sample`.
+    sample : array_like (n, d)
+        The initial sample.
+    initial_disc : float
+        Centered discrepancy of the `sample`.
+
+    Returns
+    -------
+    discrepancy : float
+        Centered discrepancy of the sample composed of `x_new` and `sample`.
+
+    """
+    sample = np.asarray(sample)
+    x_new = np.asarray(x_new)
+
+    n = len(sample) + 1
+    abs_ = abs(x_new - 0.5)
+
+    disc1 = - 2 / n * np.prod(1 + 1 / 2 * abs_ - 1 / 2 * abs_ ** 2)
+    disc2 = 2 / (n ** 2) * np.sum(np.prod(1 + 1 / 2 * abs_ +
+                                          1 / 2 * abs(sample - 0.5) -
+                                          1 / 2 * abs(x_new - sample),
+                                          axis=1))
+    disc3 = 1 / (n ** 2) * np.prod(1 + abs_)
+
+    return initial_disc + disc1 + disc2 + disc3
+
+
+def _perturb_discrepancy(sample, i1, i2, k, disc):
+    """Centered discrepancy after and elementary perturbation on a LHS.
+
+    An elementary perturbation consists of an exchange of coordinates between
+    two points: ``sample[i1, k] <-> sample[i2, k]``. By construction,
+    this operation conserves the LHS properties.
+
+    Parameters
+    ----------
+    sample : array_like (n, d)
+        The sample (before permutation) to compute the discrepancy from.
+    i1 : int
+        The first line of the elementary permutation.
+    i2 : int
+        The second line of the elementary permutation.
+    k : int
+        The column of the elementary permutation.
+    disc : float
+        Centered discrepancy of the design before permutation.
+
+    Returns
+    -------
+    discrepancy : float
+        Centered discrepancy.
+
+    References
+    ----------
+    .. [1] Jin et al. "An efficient algorithm for constructing optimal design
+       of computer experiments", Journal of Statistical Planning and
+       Inference, 2005.
+
+    """
+    sample = np.asarray(sample)
+    n = sample.shape[0]
+
+    z_ij = sample - 0.5
+
+    # Eq (19)
+    c_i1j = 1. / n ** 2. * np.prod(0.5 * (2. + abs(z_ij[i1, :]) +
+                                          abs(z_ij) -
+                                          abs(z_ij[i1, :] - z_ij)),
+                                   axis=1)
+    c_i2j = 1. / n ** 2. * np.prod(0.5 * (2. + abs(z_ij[i2, :]) +
+                                          abs(z_ij) -
+                                          abs(z_ij[i2, :] - z_ij)),
+                                   axis=1)
+
+    # Eq (20)
+    c_i1i1 = (1. / n ** 2 * np.prod(1 + abs(z_ij[i1, :])) -
+              2. / n * np.prod(1. + 0.5 * abs(z_ij[i1, :]) -
+                               0.5 * z_ij[i1, :] ** 2))
+    c_i2i2 = (1. / n ** 2 * np.prod(1 + abs(z_ij[i2, :])) -
+              2. / n * np.prod(1. + 0.5 * abs(z_ij[i2, :]) -
+                               0.5 * z_ij[i2, :] ** 2))
+
+    # Eq (22), typo in the article in the denominator i2 -> i1
+    num = (2 + abs(z_ij[i2, k]) + abs(z_ij[:, k]) -
+           abs(z_ij[i2, k] - z_ij[:, k]))
+    denum = (2 + abs(z_ij[i1, k]) + abs(z_ij[:, k]) -
+             abs(z_ij[i1, k] - z_ij[:, k]))
+    gamma = num / denum
+
+    # Eq (23)
+    c_p_i1j = gamma * c_i1j
+    # Eq (24)
+    c_p_i2j = c_i2j / gamma
+
+    alpha = (1 + abs(z_ij[i2, k])) / (1 + abs(z_ij[i1, k]))
+    beta = (2 - abs(z_ij[i2, k])) / (2 - abs(z_ij[i1, k]))
+
+    g_i1 = np.prod(1. + abs(z_ij[i1, :]))
+    g_i2 = np.prod(1. + abs(z_ij[i2, :]))
+    h_i1 = np.prod(1. + 0.5 * abs(z_ij[i1, :]) - 0.5 * (z_ij[i1, :] ** 2))
+    h_i2 = np.prod(1. + 0.5 * abs(z_ij[i2, :]) - 0.5 * (z_ij[i2, :] ** 2))
+
+    # Eq (25), typo in the article g is missing
+    c_p_i1i1 = ((g_i1 * alpha) / (n ** 2) - 2. * alpha * beta * h_i1 / n)
+    # Eq (26), typo in the article n ** 2
+    c_p_i2i2 = ((g_i2 / ((n ** 2) * alpha)) - (2. * h_i2 / (n * alpha * beta)))
+
+    # Eq (26)
+    sum_ = c_p_i1j - c_i1j + c_p_i2j - c_i2j
+
+    mask = np.ones(n, dtype=bool)
+    mask[[i1, i2]] = False
+    sum_ = sum(sum_[mask])
+
+    disc_ep = (disc + c_p_i1i1 - c_i1i1 + c_p_i2i2 - c_i2i2 + 2 * sum_)
+
+    return disc_ep
+
+
+def primes_from_2_to(n):
+    """Prime numbers from 2 to *n*.
+
+    Parameters
+    ----------
+    n : int
+        Sup bound with ``n >= 6``.
+
+    Returns
+    -------
+    primes : list(int)
+        Primes in ``2 <= p < n``.
+
+    References
+    ----------
+    .. [1] `StackOverflow <https://stackoverflow.com/questions/2068372>`_.
+
+    """
+    sieve = np.ones(n // 3 + (n % 6 == 2), dtype=bool)
+    for i in range(1, int(n ** 0.5) // 3 + 1):
+        k = 3 * i + 1 | 1
+        sieve[k * k // 3::2 * k] = False
+        sieve[k * (k - 2 * (i & 1) + 4) // 3::2 * k] = False
+    return np.r_[2, 3, ((3 * np.nonzero(sieve)[0][1:] + 1) | 1)]
+
+
+def n_primes(n):
+    """List of the n-first prime numbers.
+
+    Parameters
+    ----------
+    n : int
+        Number of prime numbers wanted.
+
+    Returns
+    -------
+    primes : list(int)
+        List of primes.
+
+    """
+    primes = [2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59,
+              61, 67, 71, 73, 79, 83, 89, 97, 101, 103, 107, 109, 113, 127,
+              131, 137, 139, 149, 151, 157, 163, 167, 173, 179, 181, 191, 193,
+              197, 199, 211, 223, 227, 229, 233, 239, 241, 251, 257, 263, 269,
+              271, 277, 281, 283, 293, 307, 311, 313, 317, 331, 337, 347, 349,
+              353, 359, 367, 373, 379, 383, 389, 397, 401, 409, 419, 421, 431,
+              433, 439, 443, 449, 457, 461, 463, 467, 479, 487, 491, 499, 503,
+              509, 521, 523, 541, 547, 557, 563, 569, 571, 577, 587, 593, 599,
+              601, 607, 613, 617, 619, 631, 641, 643, 647, 653, 659, 661, 673,
+              677, 683, 691, 701, 709, 719, 727, 733, 739, 743, 751, 757, 761,
+              769, 773, 787, 797, 809, 811, 821, 823, 827, 829, 839, 853, 857,
+              859, 863, 877, 881, 883, 887, 907, 911, 919, 929, 937, 941, 947,
+              953, 967, 971, 977, 983, 991, 997][:n]
+
+    if len(primes) < n:
+        big_number = 2000
+        while 'Not enough primes':
+            primes = primes_from_2_to(big_number)[:n]
+            if len(primes) == n:
+                break
+            big_number += 1000
+
+    return primes
+
+
+def van_der_corput(n, base=2, start_index=0, scramble=False, seed=None):
+    """Van der Corput sequence.
+
+    Pseudo-random number generator based on a b-adic expansion.
+
+    Scrambling uses permutations of the remainders (see [1]_). Multiple
+    permutations are applied to construct a point. The sequence of
+    permutations has to be the same for all points of the sequence.
+
+    Parameters
+    ----------
+    n : int
+        Number of element of the sequence.
+    base : int
+        Base of the sequence.
+    start_index : int
+        Index to start the sequence from.
+    scramble: bool
+        If True, use Owen scrambling. Otherwise no scrambling is done.
+        Default is True.
+    seed : {int or `numpy.random.RandomState` instance}, optional
+        If `seed` is not specified the `numpy.random.RandomState`
+        singleton is used.
+        If `seed` is an int, a new ``RandomState`` instance is used,
+        seeded with `seed`.
+        If `seed` is already a ``RandomState`` instance, then that
+        instance is used.
+
+    Returns
+    -------
+    sequence : list (n,)
+        Sequence of Van der Corput.
+
+    References
+    ----------
+    .. [1] A. B. Owen. "A randomized Halton algorithm in R",
+       arXiv:1706.02808, 2017.
+
+    """
+    rng = check_random_state(seed)
+    sequence = np.zeros(n)
+
+    quotient = np.arange(start_index, start_index + n)
+    b2r = 1 / base
+
+    while (1 - b2r) < 1:
+        remainder = quotient % base
+
+        if scramble:
+            # permutation must be the same for all points of the sequence
+            perm = rng.permutation(base)
+            remainder = perm[np.array(remainder).astype(int)]
+
+        sequence += remainder * b2r
+        b2r /= base
+        quotient = (quotient - remainder) / base
+
+    return sequence
+
+
+class QMCEngine(ABC):
+    """A generic Quasi-Monte Carlo sampler class meant for subclassing.
+
+    QMCEngine is a base class to construct a specific Quasi-Monte Carlo sampler.
+    It cannot be used directly as a sampler.
+
+    Parameters
+    ----------
+    d : int
+        Dimension of the parameter space.
+    seed : {int or `numpy.random.RandomState` instance}, optional
+        If `seed` is not specified, the `numpy.random.RandomState`
+        singleton is used.
+        If `seed` is an int, a new ``RandomState`` instance is used,
+        seeded with `seed`.
+        If `seed` is already a ``RandomState`` instance, then that
+        instance is used.
+
+    Notes
+    -----
+    By convention samples are distributed over the half-open interval
+    ``[0, 1)``. Instances of the class can access the attributes: ``d`` for
+    the dimension; and ``rng`` for the random number generator (used for the
+    ``seed``).
+
+    **Subclassing**
+
+    When subclassing `QMCEngine` to create a new sampler,  ``__init__`` and
+    ``random`` must be redefined.
+
+    * ``__init__(d, seed=None)``: at least fix the dimension. If the sampler
+      does not take advantage of a ``seed`` (deterministic methods like
+      Halton), this parameter can be omitted.
+    * ``random(n)``: draw ``n`` from the engine.
+
+    Optionally, two other methods can be overwritten by subclasses:
+
+    * ``reset``: Reset the engine to it's original state.
+    * ``fast_forward``: If the sequence is deterministic (like Halton
+      sequence), then ``fast_forward(n)`` is skipping the ``n`` first draw.
+
+    Examples
+    --------
+    To create a random sampler based on ``np.random.random``, we would do the
+    following:
+
+    >>> from scipy.stats import qmc
+    >>> class RandomEngine(qmc.QMCEngine):
+    ...     def __init__(self, d, seed):
+    ...         super().__init__(d=d, seed=seed)
+    ...         self.rng_seed_state = self.rng.get_state()
+    ...
+    ...
+    ...     def random(self, n=1):
+    ...         return self.rng.random((n, self.d))
+    ...
+    ...
+    ...     def reset(self):
+    ...         self.rng.set_state(self.rng_seed_state)
+    ...         return self
+    ...
+    ...
+    ...     def fast_forward(self, n):
+    ...         self.rng.random((n, self.d))
+    ...         return self
+
+    After subclassing `QMCEngine` to define the sampling strategy we want to use,
+    we can create an instance to sample from.
+
+    >>> engine = RandomEngine(2, seed=12345)
+    >>> engine.random(5)
+    array([[0.92961609, 0.31637555],
+           [0.18391881, 0.20456028],
+           [0.56772503, 0.5955447 ],
+           [0.96451452, 0.6531771 ],
+           [0.74890664, 0.65356987]])
+
+    We can also reset the state of the generator and resample again.
+
+    >>> _ = engine.reset()
+    >>> engine.random(5)
+    array([[0.92961609, 0.31637555],
+           [0.18391881, 0.20456028],
+           [0.56772503, 0.5955447 ],
+           [0.96451452, 0.6531771 ],
+           [0.74890664, 0.65356987]])
+
+    """
+
+    @abstractmethod
+    def __init__(self, d, seed=None):
+        self.d = d
+        self.rng = check_random_state(seed)
+        self.num_generated = 0
+
+    @abstractmethod
+    def random(self, n=1):
+        """Draw `n` in the half-open interval ``[0, 1)``.
+
+        Parameters
+        ----------
+        n : int
+            Number of samples to generate in the parameter space.
+            Default is 1.
+
+        Returns
+        -------
+        sample : array_like (n, d)
+            QMC sample.
+        """
+        # self.num_generated += n
+
+    def reset(self):
+        """Reset the engine to base state.
+
+        Returns
+        -------
+        engine: QMCEngine
+            Engine reset to its base state.
+
+        """
+        self.num_generated = 0
+        return self
+
+    def fast_forward(self, n):
+        """Fast-forward the sequence by `n` positions.
+
+        Parameters
+        ----------
+        n: int
+            Number of points to skip in the sequence.
+        """
+        self.num_generated += n
+        return self
+
+
+class Halton(QMCEngine):
+    """Halton sequence.
+
+    Pseudo-random number generator that generalize the Van der Corput sequence
+    for multiple dimensions. The Halton sequence uses the base-two Van der Corput
+    sequence for the first dimension, base-three for its second and base-:math:`n` for
+    its n-dimension.
+
+    Parameters
+    ----------
+    d : int
+        Dimension of the parameter space.
+    scramble: bool
+        If True, use Owen scrambling. Otherwise no scrambling is done.
+        Default is True.
+    seed : {int or `numpy.random.RandomState` instance}, optional
+        If `seed` is not specified, the `numpy.random.RandomState`
+        singleton is used.
+        If `seed` is an int, a new ``RandomState`` instance is used,
+        seeded with `seed`.
+        If `seed` is already a ``RandomState`` instance, then that
+        instance is used.
+
+    Notes
+    -----
+    The Halton sequence has severe striping artifacts for even modestly
+    large dimensions. These can be ameliorated by scrambling. Scrambling
+    also supports replication-based error estimates and extends
+    applicabiltiy to unbounded integrands.
+
+    References
+    ----------
+    .. [1] Halton, "On the efficiency of certain quasi-random sequences of
+       points in evaluating multi-dimensional integrals", Numerische
+       Mathematik, 1960.
+    .. [2] A. B. Owen. "A randomized Halton algorithm in R",
+       arXiv:1706.02808, 2017.
+
+    Examples
+    --------
+    Generate samples from a low discrepancy sequence of Halton.
+
+    >>> from scipy.stats import qmc
+    >>> sampler = qmc.Halton(d=2, scramble=False)
+    >>> sample = sampler.random(n=5)
+    >>> sample
+    array([[0.        , 0.        ],
+           [0.5       , 0.33333333],
+           [0.25      , 0.66666667],
+           [0.75      , 0.11111111],
+           [0.125     , 0.44444444]])
+
+    Compute the quality of the sample using the discrepancy criterion.
+
+    >>> qmc.discrepancy(sample)
+    0.088893711419753
+
+    If some wants to continue an existing design, extra points can be obtained
+    by calling again `random`. Alternatively, you can skip some points like:
+
+    >>> sampler.fast_forward(5)
+    >>> sample_continued = sampler.random(n=5)
+    >>> sample_continued
+    array([[0.3125    , 0.37037037],
+           [0.8125    , 0.7037037 ],
+           [0.1875    , 0.14814815],
+           [0.6875    , 0.48148148],
+           [0.4375    , 0.81481481]])
+
+    Finally, samples can be scaled to bounds.
+
+    >>> bounds = [[0, 2], [10, 5]]
+    >>> qmc.scale(sample_continued, bounds)
+    array([[3.125     , 3.11111111],
+           [8.125     , 4.11111111],
+           [1.875     , 2.44444444],
+           [6.875     , 3.44444444],
+           [4.375     , 4.44444444]])
+
+    """
+
+    def __init__(self, d, scramble=True, seed=None):
+        super().__init__(d=d)
+        self.seed = seed
+        self.base = n_primes(d)
+        self.scramble = scramble
+
+    def random(self, n=1):
+        """Draw `n` in the half-open interval ``[0, 1)``.
+
+        Parameters
+        ----------
+        n : int
+            Number of samples to generate in the parameter space. Default is 1.
+
+        Returns
+        -------
+        sample : array_like (n, d)
+            QMC sample.
+
+        """
+        # Generate a sample using a Van der Corput sequence per dimension.
+        # important to have ``type(bdim) == int`` for performance reason
+        sample = [van_der_corput(n, int(bdim), self.num_generated,
+                                 scramble=self.scramble, seed=self.seed)
+                  for bdim in self.base]
+
+        self.num_generated += n
+        return np.array(sample).T
+
+
+class OrthogonalLatinHypercube(QMCEngine):
+    """Orthogonal array-based Latin hypercube sampling (OA-LHS).
+
+    In addition to the constraints from the Latin Hypercube, an orthogonal
+    array of size `n` is defined and only one point is allowed per subspace.
+
+    Parameters
+    ----------
+    d : int
+        Dimension of the parameter space.
+    seed : {int or `numpy.random.RandomState` instance}, optional
+        If `seed` is not specified the `numpy.random.RandomState`
+        singleton is used.
+        If `seed` is an int, a new ``RandomState`` instance is used,
+        seeded with `seed`.
+        If `seed` is already a ``RandomState`` instance, then that
+        instance is used.
+
+    References
+    ----------
+    .. [1] Art B. Owen, "Orthogonal arrays for computer experiments,
+       integration and visualization", Statistica Sinica, 1992.
+
+    Examples
+    --------
+    Generate samples from an orthogonal latin hypercube generator.
+
+    >>> from scipy.stats import qmc
+    >>> sampler = qmc.OrthogonalLatinHypercube(d=2, seed=12345)
+    >>> sample = sampler.random(n=5)
+    >>> sample
+    array([[0.18592322, 0.77846875],
+           [0.64091206, 0.18474763],
+           [0.91354501, 0.80535794],
+           [0.43678376, 0.27478424],
+           [0.26327511, 0.43099467]])
+
+    Compute the quality of the sample using the discrepancy criterion.
+
+    >>> qmc.discrepancy(sample)
+    0.02004864477993462
+
+    Finally, samples can be scaled to bounds.
+
+    >>> bounds = [[0, 2], [10, 5]]
+    >>> qmc.scale(sample, bounds)
+    array([[1.85923219, 4.33540624],
+           [6.40912056, 2.55424289],
+           [9.13545006, 4.41607383],
+           [4.36783762, 2.82435271],
+           [2.63275111, 3.29298401]])
+
+    """
+
+    def __init__(self, d, seed=None):
+        super().__init__(d=d, seed=seed)
+
+    def random(self, n=1):
+        """Draw `n` in the half-open interval ``[0, 1)``.
+
+        Parameters
+        ----------
+        n : int
+            Number of samples to generate in the parameter space. Default is 1.
+
+        Returns
+        -------
+        sample : array_like (n, d)
+            OLHS sample.
+
+        """
+        sample = []
+        step = 1.0 / n
+
+        for _ in range(self.d):
+            # Enforce a unique point per grid
+            j = np.arange(n) * step
+            temp = j + self.rng.uniform(low=0, high=step, size=n)
+            self.rng.shuffle(temp)
+
+            sample.append(temp)
+
+        return np.array(sample).T
+
+
+class LatinHypercube(QMCEngine):
+    """Latin hypercube sampling (LHS).
+
+    A Latin hypercube sample [1]_ generates :math:`n` points in
+    :math:`[0,1)^{d}`. Each univariate marginal distribution is stratified,
+    placing exactly one point in :math:`[j/n, (j+1)/n)` for
+    :math:`j=0,1,...,n-1`. They are still applicable when :math:`n << d`.
+    LHS is extremely effective on integrands that are nearly additive [2]_.
+    LHS on :math:`n` points never has more variance than plain MC on
+    :math:`n-1` points [3]_. There is a central limit theorem for plain
+    LHS [4]_, but not necessarily for optimized LHS.
+
+    Parameters
+    ----------
+    d : int
+        Dimension of the parameter space.
+    centered : bool
+        Center the point within the multi-dimensional grid. Default is False.
+    seed : {int or `numpy.random.RandomState` instance}, optional
+        If `seed` is not specified the `numpy.random.RandomState`
+        singleton is used.
+        If `seed` is an int, a new ``RandomState`` instance is used,
+        seeded with `seed`.
+        If `seed` is already a ``RandomState`` instance, then that
+        instance is used.
+
+    References
+    ----------
+    .. [1] Mckay et al., "A Comparison of Three Methods for Selecting Values
+       of Input Variables in the Analysis of Output from a Computer Code",
+       Technometrics, 1979.
+    .. [2] M. Stein, "Large sample properties of simulations using Latin
+       hypercube sampling." Technometrics 29, no. 2: 143-151, 1987.
+    .. [3] A. B. Owen, "Monte Carlo variance of scrambled net quadrature."
+       SIAM Journal on Numerical Analysis 34, no. 5: 1884-1910, 1997
+    .. [4]  Loh, W.-L. "On Latin hypercube sampling." The annals of statistics
+       24, no. 5: 2058-2080, 1996.
+
+    Examples
+    --------
+    Generate samples from a Latin hypercube generator.
+
+    >>> from scipy.stats import qmc
+    >>> sampler = qmc.LatinHypercube(d=2, seed=12345)
+    >>> sample = sampler.random(n=5)
+    >>> sample
+    array([[0.01407678, 0.53672489],
+           [0.36321624, 0.75908794],
+           [0.28645499, 0.48089106],
+           [0.2070971 , 0.46936458],
+           [0.45021867, 0.66928603]])
+
+    Compute the quality of the sample using the discrepancy criterion.
+
+    >>> qmc.discrepancy(sample)
+    0.1271335273223828
+
+    Finally, samples can be scaled to bounds.
+
+    >>> bounds = [[0, 2], [10, 5]]
+    >>> qmc.scale(sample, bounds)
+    array([[0.14076781, 3.61017467],
+           [3.63216238, 4.27726383],
+           [2.86454994, 3.44267318],
+           [2.07097096, 3.40809374],
+           [4.50218672, 4.00785808]])
+
+    """
+
+    def __init__(self, d, centered=False, seed=None):
+        super().__init__(d=d, seed=seed)
+        self.centered = centered
+
+    def random(self, n=1):
+        """Draw `n` in the half-open interval ``[0, 1)``.
+
+        Parameters
+        ----------
+        n : int
+            Number of samples to generate in the parameter space. Default is 1.
+
+        Returns
+        -------
+        sample : array_like (n, d)
+            LHS sample.
+
+        """
+        if self.centered:
+            r = 0.5
+        else:
+            r = self.rng.random_sample((n, self.d))
+
+        try:
+            rg_integers = self.rng.integers
+        except AttributeError:
+            rg_integers = self.rng.randint
+        q = rg_integers(low=1, high=n, size=(n, self.d))
+
+        return 1. / n * (q - r)
+
+
+class OptimalDesign(QMCEngine):
+    """Optimal design.
+
+    Optimize the design by doing random permutations to lower the centered
+    discrepancy. If `optimization` is False, `niter` design are generated and
+    the one with lowest centered discrepancy is return. This option is faster.
+    Otherwise,
+
+    Centered discrepancy-based design shows better space filling robustness
+    toward 2D and 3D subprojections. Distance-based design shows better space
+    filling but less robustness to subprojections.
+
+    Parameters
+    ----------
+    d : int
+        Dimension of the parameter space.
+    start_design : array_like (n, d), optional
+        Initial design of experiment to optimize.
+    niter : int
+        Number of iterations to perform. Default is 1.
+    force : bool
+        If `optimization`, force *basinhopping* optimization. Otherwise
+        grid search is used. Default is False.
+    optimization : bool
+        Optimal design using global optimization or random generation of
+        `niter` samples. Default is True.
+    seed : {int or `numpy.random.RandomState` instance}, optional
+        If `seed` is not specified the `numpy.random.RandomState`
+        singleton is used.
+        If `seed` is an int, a new ``RandomState`` instance is used,
+        seeded with `seed`.
+        If `seed` is already a ``RandomState`` instance, then that
+        instance is used.
+
+    References
+    ----------
+    .. [1] Fang et al. Design and modeling for computer experiments,
+       Computer Science and Data Analysis Series, 2006.
+    .. [2] Damblin et al., "Numerical studies of space filling designs:
+       optimization of Latin Hypercube Samples and subprojection properties",
+       Journal of Simulation, 2013.
+
+    Examples
+    --------
+    Generate samples from an optimal design.
+
+    >>> from scipy.stats import qmc
+    >>> sampler = qmc.OptimalDesign(d=2, seed=12345)
+    >>> sample = sampler.random(n=5)
+    >>> sample
+    array([[0.64091206, 0.77846875],
+           [0.43678376, 0.18474763],
+           [0.18592322, 0.80535794],
+           [0.91354501, 0.27478424],
+           [0.26327511, 0.43099467]])
+
+    Compute the quality of the sample using the discrepancy criterion.
+
+    >>> qmc.discrepancy(sample)
+    0.019688311022535432
+
+    You can possibly improve the quality of the sample by performing more
+    optimization iterations by using `niter`:
+
+    >>> sampler_2 = qmc.OptimalDesign(d=2, niter=2, seed=12345)
+    >>> sample_2 = sampler_2.random(n=5)
+    >>> qmc.discrepancy(sample_2)
+    0.019607673478802434
+
+    Finally, samples can be scaled to bounds.
+
+    >>> bounds = [[0, 2], [10, 5]]
+    >>> qmc.scale(sample, bounds)
+    array([[6.40912056, 4.33540624],
+           [4.36783762, 2.55424289],
+           [1.85923219, 4.41607383],
+           [9.13545006, 2.82435271],
+           [2.63275111, 3.29298401]])
+
+    """
+
+    def __init__(self, d, start_design=None, niter=1, force=False,
+                 optimization=True, seed=None):
+        super().__init__(d=d, seed=seed)
+        self.start_design = start_design
+        self.niter = niter
+        self.force = force
+        self.optimization = optimization
+
+        self.best_doe = self.start_design
+        if self.start_design is not None:
+            self.best_disc = discrepancy(self.start_design)
+        else:
+            self.best_disc = np.inf
+
+        self.olhs = OrthogonalLatinHypercube(self.d, seed=self.rng)
+
+    def random(self, n=1):
+        """Draw `n` in the half-open interval ``[0, 1)``.
+
+        Parameters
+        ----------
+        n : int
+            Number of samples to generate in the parameter space. Default is 1.
+
+        Returns
+        -------
+        sample : array_like (n, d)
+            Optimal sample.
+
+        """
+        if self.optimization:
+            if self.best_doe is None:
+                self.best_doe = self.olhs.random(n)
+                self.best_disc = discrepancy(self.best_doe)
+
+            def _perturb_best_doe(x):
+                """Perturb the DoE and keep track of the best DoE.
+
+                Parameters
+                ----------
+                x : list of int
+                    It is a list of:
+                        idx : int
+                            Index value of the components to compute
+
+                Returns
+                -------
+                discrepancy : float
+                    Centered discrepancy.
+
+                """
+                # Perturb the DoE
+                doe = self.best_doe.copy()
+                col, row_1, row_2 = np.round(x).astype(int)
+                doe[row_1, col], doe[row_2, col] = doe[row_2, col], \
+                                                   doe[row_1, col]
+
+                disc = _perturb_discrepancy(self.best_doe, row_1, row_2, col,
+                                            self.best_disc)
+
+                if disc < self.best_disc:
+                    self.best_disc = disc
+                    self.best_doe = doe
+
+                return disc
+
+            # Total number of possible design
+            complexity = self.d * n ** 2
+
+            if (complexity > 1e6) or self.force:
+                bounds_optim = ([0, self.d - 1],
+                                [0, n - 1],
+                                [0, n - 1])
+            else:
+                bounds_optim = (slice(0, self.d - 1, 1),
+                                slice(0, n - 1, 1),
+                                slice(0, n - 1, 1))
+
+            for _ in range(self.niter):
+                if (complexity > 1e6) or self.force:
+                    minimizer_kwargs = {"method": "L-BFGS-B",
+                                        "bounds": bounds_optim}
+                    _ = basinhopping(_perturb_best_doe, [0, 0, 0], niter=100,
+                                     minimizer_kwargs=minimizer_kwargs)
+                else:
+                    _ = brute(_perturb_best_doe, ranges=bounds_optim,
+                              finish=None)
+        else:
+            for _ in range(self.niter):
+                doe = self.olhs.random(n)
+                disc = discrepancy(doe)
+                if disc < self.best_disc:
+                    self.best_disc = disc
+                    self.best_doe = doe
+
+        return self.best_doe
+
+
+class Sobol(QMCEngine):
+    """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 :math:`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
+    ----------
+    d : int
+        Dimensionality of the sequence. Max dimensionality is 21201.
+    scramble : bool
+        If True, use Owen scrambling. Otherwise no scrambling is done.
+        Default is True.
+    seed : {int or `numpy.random.RandomState` instance}, optional
+        If `seed` is not specified the `numpy.random.RandomState`
+        singleton is used.
+        If `seed` is an int, a new ``RandomState`` instance is used,
+        seeded with `seed`.
+        If `seed` is already a ``RandomState`` instance, then that
+        instance is used.
+
+    Notes
+    -----
+    Sobol' sequences [1]_ provide :math:`n=2^m` low discrepancy points in
+    :math:`[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 :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 will repeat.
+       Currently :math:`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.
+
+    >>> bounds = [[0, 2], [10, 5]]
+    >>> qmc.scale(sample_continued, bounds)
+    array([[3.75 , 3.125],
+           [8.75 , 4.625],
+           [6.25 , 2.375],
+           [1.25 , 3.875]])
+
+    """
+
+    MAXDIM = _MAXDIM
+    MAXBIT = _MAXBIT
+
+    def __init__(self, d, scramble=True, seed=None):
+        if d > self.MAXDIM:
+            raise ValueError(
+                "Maximum supported dimensionality is {}.".format(self.MAXDIM)
+            )
+        super().__init__(d=d, seed=seed)
+
+        # initialize direction numbers
+        initialize_direction_numbers()
+
+        # v is d x MAXBIT matrix
+        self._sv = np.zeros((d, self.MAXBIT), dtype=int)
+        initialize_v(self._sv, d)
+
+        if not scramble:
+            self._shift = np.zeros(d, dtype=int)
+        else:
+            self._scramble()
+
+        self._quasi = self._shift.copy()
+        self._first_point = (self._quasi / 2 ** self.MAXBIT).reshape(1, -1)
+
+    def _scramble(self):
+        """Scramble the sequence."""
+        try:
+            rg_integers = self.rng.integers
+        except AttributeError:
+            rg_integers = self.rng.randint
+        # Generate shift vector
+        self._shift = np.dot(
+            rg_integers(2, size=(self.d, self.MAXBIT)),
+            2 ** np.arange(self.MAXBIT),
+        )
+        self._quasi = self._shift.copy()
+        # Generate lower triangular matrices (stacked across dimensions)
+        ltm = np.tril(rg_integers(2, size=(self.d, self.MAXBIT, self.MAXBIT)))
+        _cscramble(self.d, ltm, self._sv)
+        self.num_generated = 0
+
+    def random(self, n=1):
+        """Draw next point(s) in the Sobol' sequence.
+
+        Parameters
+        ----------
+        n : int
+            Number of samples to generate in the parameter space. Default is 1.
+
+        Returns
+        -------
+        sample : array_like (n, d)
+            Sobol' sample.
+
+        """
+        sample = np.empty((n, self.d), dtype=float)
+
+        if self.num_generated == 0:
+            # verify n is 2**n
+            if not (n & (n - 1) == 0):
+                warnings.warn("The balance properties of Sobol' points require"
+                              " n to be a power of 2.")
+
+            if n == 1:
+                sample = self._first_point
+            else:
+                _draw(n - 1, self.num_generated, self.d, self._sv,
+                      self._quasi, sample)
+                sample = np.concatenate([self._first_point, sample])[:n]
+        else:
+            _draw(n, self.num_generated - 1, self.d, self._sv,
+                  self._quasi, sample)
+
+        self.num_generated += n
+        return sample
+
+    def random_base2(self, m=1):
+        """Draw point(s) from the Sobol' sequence.
+
+        This function draws :math:`n=2^m` points in the parameter space
+        ensuring the balance properties of the sequence.
+
+        Parameters
+        ----------
+        m : int
+            Logarithm in base 2 of the number of samples; i.e., n = 2^m.
+            Default is 1.
+
+        Returns
+        -------
+        sample : array_like (n, d)
+            Sobol' sample.
+
+        """
+        n = 2 ** m
+
+        total_n = self.num_generated + n
+        if not (total_n & (total_n - 1) == 0):
+            raise ValueError("The balance properties of Sobol' points require "
+                             "n to be a power of 2. {0} points have been "
+                             "previously generated, then: n={0}+2**{1}={2}. "
+                             "If you still want to do this, the function "
+                             "'Sobol.random()' can be used."
+                             .format(self.num_generated, m, total_n))
+
+        return self.random(n)
+
+    def reset(self):
+        """Reset the engine to base state.
+
+        Returns
+        -------
+        engine: Sobol
+            Engine reset to its base state.
+
+        """
+        self._quasi = self._shift.copy()
+        self.num_generated = 0
+        return self
+
+    def fast_forward(self, n):
+        """Fast-forward the sequence by `n` positions.
+
+        Parameters
+        ----------
+        n: int
+            Number of points to skip in the sequence.
+
+        Returns
+        -------
+        engine: Sobol
+            The fast-forwarded engine.
+
+        """
+        if self.num_generated == 0:
+            _fast_forward(n - 1, self.num_generated, self.d,
+                          self._sv, self._quasi)
+        else:
+            _fast_forward(n, self.num_generated - 1, self.d,
+                          self._sv, self._quasi)
+        self.num_generated += n
+        return self
+
+
+def multinomial_qmc(n, pvals, engine=None, seed=None):
+    """Draw low-discreancy quasi-random samples from multinomial distribution.
+
+    Parameters
+    ----------
+    n : int
+        Number of experiments.
+    pvals: Iterable[float]
+        float vector of probabilities of size ``p``. Elements must be
+        non-negative and sum to 1.
+    engine: QMCEngine, optional
+        Quasi-Monte Carlo engine sampler. If None, Sobol' is used.
+    seed : {int or `numpy.random.RandomState` instance}, optional
+        If `seed` is not specified the `numpy.random.RandomState`
+        singleton is used.
+        If `seed` is an int, a new ``RandomState`` instance is used,
+        seeded with `seed`.
+        If `seed` is already a ``RandomState`` instance, then that
+        instance is used.
+
+    Returns
+    -------
+    samples: array_like (pvals,)
+        int vector of size ``p`` summing to `n`.
+
+    """
+    if np.min(pvals) < 0:
+        raise ValueError('Elements of pvals must be non-negative.')
+    if not np.isclose(np.sum(pvals), 1):
+        raise ValueError('Elements of pvals must sum to 1.')
+
+    if engine is None:
+        engine = Sobol(1, scramble=True, seed=seed)
+    draws = engine.random(n).ravel()
+    p_cumulative = np.empty_like(pvals, dtype=float)
+    _fill_p_cumulative(np.array(pvals, dtype=float), p_cumulative)
+    sample = np.zeros_like(pvals, dtype=int)
+    _categorize(draws, p_cumulative, sample)
+    return sample
+
+
+class NormalQMC(QMCEngine):
+    """Engine for QMC sampling from a multivariate normal :math:`N(0, I_d)`.
+
+    Parameters
+    ----------
+    d: int
+        The dimension of the samples.
+    inv_transform: bool
+        If True, use inverse transform instead of Box-Muller. Default is True.
+    engine: QMCEngine, optional
+        Quasi-Monte Carlo engine sampler. If None, Sobol' is used.
+    seed : {int or `numpy.random.RandomState` instance}, optional
+        If `seed` is not specified the `numpy.random.RandomState`
+        singleton is used.
+        If `seed` is an int, a new ``RandomState`` instance is used,
+        seeded with `seed`.
+        If `seed` is already a ``RandomState`` instance, then that
+        instance is used.
+
+    Examples
+    --------
+    >>> import matplotlib.pyplot as plt
+    >>> from scipy.stats import qmc
+    >>> engine = qmc.NormalQMC(2)
+    >>> sample = engine.random(512)
+    >>> plt.scatter(sample[:, 0], sample[:, 1])
+    >>> plt.show()
+
+    """
+
+    def __init__(self, d, inv_transform=True, engine=None, seed=None):
+        super().__init__(d=d, seed=seed)
+        self._inv_transform = inv_transform
+        if not inv_transform:
+            # to apply Box-Muller, we need an even number of dimensions
+            engine_dim = 2 * math.ceil(d / 2)
+        else:
+            engine_dim = d
+
+        if engine is None:
+            self.engine = Sobol(d=engine_dim, scramble=True, seed=seed)
+        else:
+            self.engine = engine
+
+    def random(self, n=1):
+        """Draw `n` QMC samples from the standard Normal.
+
+        Parameters
+        ----------
+        n : int
+            Number of samples to generate in the parameter space. Default is 1.
+
+        Returns
+        -------
+        sample : array_like (n, d)
+            Sample.
+
+        """
+        # get base samples
+        samples = self.engine.random(n)
+        if self._inv_transform:
+            # apply inverse transform
+            # (values to close to 0/1 result in inf values)
+            return norm.ppf(0.5 + (1 - 1e-10) * (samples - 0.5))
+        else:
+            # apply Box-Muller transform (note: indexes starting from 1)
+            even = np.arange(0, samples.shape[-1], 2)
+            Rs = np.sqrt(-2 * np.log(samples[:, even]))
+            thetas = 2 * math.pi * samples[:, 1 + even]
+            cos = np.cos(thetas)
+            sin = np.sin(thetas)
+            transf_samples = np.stack([Rs * cos, Rs * sin],
+                                      -1).reshape(n, -1)
+            # make sure we only return the number of dimension requested
+            return transf_samples[:, : self.d]
+
+
+class MultivariateNormalQMC(QMCEngine):
+    r"""QMC sampling from a multivariate Normal :math:`N(\mu, \Sigma)`.
+
+    Parameters
+    ----------
+    mean: array_like (d,)
+        The mean vector.
+    cov: array_like (d, d)
+        The covariance matrix.
+    inv_transform: bool
+        If True, use inverse transform instead of Box-Muller. Default is True.
+    engine: QMCEngine, optional
+        Quasi-Monte Carlo engine sampler. If None, Sobol' is used.
+    seed : {int or `numpy.random.RandomState` instance}, optional
+        If `seed` is not specified the `numpy.random.RandomState`
+        singleton is used.
+        If `seed` is an int, a new ``RandomState`` instance is used,
+        seeded with `seed`.
+        If `seed` is already a ``RandomState`` instance, then that
+        instance is used.
+
+    Examples
+    --------
+    >>> import matplotlib.pyplot as plt
+    >>> from scipy.stats import qmc
+    >>> engine = qmc.MultivariateNormalQMC(mean=[0, 5], cov=[[1, 0], [0, 1]])
+    >>> sample = engine.random(512)
+    >>> plt.scatter(sample[:, 0], sample[:, 1])
+    >>> plt.show()
+
+    """
+
+    def __init__(self, mean, cov, inv_transform=True, engine=None, seed=None):
+        # check for square/symmetric cov matrix and mean vector has the same d
+        mean = np.array(mean, copy=False, ndmin=1)
+        cov = np.array(cov, copy=False, ndmin=2)
+        if not mean.shape[0] == cov.shape[0]:
+            raise ValueError("Dimension mismatch between mean and covariance.")
+        if not np.allclose(cov, cov.transpose()):
+            raise ValueError("Covariance matrix is not symmetric.")
+
+        super().__init__(d=mean.shape[0])
+        self._mean = mean
+        self._normal_engine = NormalQMC(d=self.d, inv_transform=inv_transform,
+                                        engine=engine, seed=seed)
+        # compute Cholesky decomp; if it fails, do the eigen decomposition
+        try:
+            self._corr_matrix = np.linalg.cholesky(cov).transpose()
+        except np.linalg.LinAlgError:
+            eigval, eigvec = np.linalg.eigh(cov)
+            if not np.all(eigval >= -1.0e-8):
+                raise ValueError("Covariance matrix not PSD.")
+            eigval = np.clip(eigval, 0.0, None)
+            self._corr_matrix = (eigvec * np.sqrt(eigval)).transpose()
+
+    def random(self, n=1) -> np.ndarray:
+        """Draw `n` QMC samples from the multivariate Normal.
+
+        Parameters
+        ----------
+        n : int
+            Number of samples to generate in the parameter space. Default is 1.
+
+        Returns
+        -------
+        sample : array_like (n, d)
+            Sample.
+
+        """
+        base_samples = self._normal_engine.random(n)
+        qmc_samples = base_samples @ self._corr_matrix + self._mean
+        return qmc_samples
diff --git a/scipy/stats/_sobol.pyx b/scipy/stats/_sobol.pyx
new file mode 100644
index 000000000000..9235e15b815a
--- /dev/null
+++ b/scipy/stats/_sobol.pyx
@@ -0,0 +1,357 @@
+from __future__ import division, absolute_import
+
+cimport cython
+cimport numpy as cnp
+
+import os
+import numpy as np
+
+# Parameters are linked to the direction numbers list.
+# See `initialize_direction_numbers` for more details.
+# Declared using DEF to be known at compilation time for ``poly`` et ``vinit``
+DEF MAXDIM = 21201  # max number of dimensions
+DEF MAXDEG = 18  # max polynomial degree
+DEF MAXBIT = 30  # max number of bits
+
+
+# Needed to be accessed with python
+cdef extern from *:
+    """
+    int MAXDIM_DEFINE = 21201;
+    int MAXDEG_DEFINE = 18;
+    int MAXBIT_DEFINE = 30;
+    """
+    int MAXDIM_DEFINE  # max number of dimensions
+    int MAXDEG_DEFINE  # max polynomial degree
+    int MAXBIT_DEFINE  # max number of bits
+
+_MAXDIM = MAXDIM_DEFINE
+_MAXDEG = MAXDEG_DEFINE
+_MAXBIT = MAXBIT_DEFINE
+
+cdef int poly[MAXDIM]
+cdef int vinit[MAXDIM][MAXDEG]
+
+cdef int LARGEST_NUMBER = 2 ** MAXBIT  # largest possible integer
+cdef float RECIPD = 1.0 / LARGEST_NUMBER  # normalization constant
+
+cdef bint is_initialized = False
+
+
+def initialize_direction_numbers():
+    """Load direction numbers.
+
+    Direction numbers obtained using the search criterion D(6)
+    up to the dimension 21201. This is the recommended choice by the authors.
+
+    Original data can be found at https://web.maths.unsw.edu.au/~fkuo/sobol/.
+    For additional details on the quantities involved, see [1].
+
+    [1] S. Joe and F. Y. Kuo. Remark on algorithm 659: Implementing sobol's
+        quasirandom sequence generator. ACM Trans. Math. Softw., 29(1):49-57,
+        Mar. 2003.
+
+    The C-code generated from putting the numbers in as literals is obscenely
+    large/inefficient. The data file was thus packaged and save as an .npz data
+    file for fast loading using the following code (this assumes that the file
+    https://web.maths.unsw.edu.au/~fkuo/sobol/new-joe-kuo-6.21201 is present in
+    the working directory):
+
+        import pandas as pd
+        import numpy as np
+
+        # read in file content
+        with open("./new-joe-kuo-6.21201", "r") as f:
+            lines = f.readlines()
+
+        rows = []
+
+        # parse data from file line by line
+        for l in lines[1:]:
+            nums = [int(n) for n in l.replace(" \n", "").split()]
+            d, s, a = nums[:3]
+            vs = {f"v{i}": int(v) for i,v in enumerate(nums[3:])}
+            rows.append({"d": d, "s": s, "a": a, **vs})
+
+
+        # read in as dataframe, explicitly use zero values
+        df = pd.DataFrame(rows).fillna(0).astype(int)
+
+        # peform conversion
+        df["poly"] = 2 * df["a"] + 2 ** df["s"] + 1
+
+        # ensure columns are properly ordered
+        vs = df[[f"v{i}" for i in range(18)]].values
+
+        # add the degenerate d=1 column (not included in the data file)
+        vs = np.vstack([vs[0][np.newaxis, :], vs])
+        poly = np.concatenate([[1], df["poly"].values])
+
+        # save as compressed .npz file to minimize size of distribution
+        np.savez_compressed("./_sobol_direction_numbers", vinit=vs, poly=poly)
+
+    """
+    global is_initialized, poly, vinit
+    if not is_initialized:
+        dns = np.load(os.path.join(os.path.dirname(__file__), "_sobol_direction_numbers.npz"))
+        poly = dns["poly"]
+        vinit = dns["vinit"]
+        is_initialized = True
+
+
+@cython.boundscheck(False)
+@cython.wraparound(False)
+cdef int bit_length(const int n):
+    cdef int bits = 0
+    cdef int nloc = n
+    while nloc != 0:
+        nloc >>= 1
+        bits += 1
+    return bits
+
+
+@cython.boundscheck(False)
+@cython.wraparound(False)
+cdef int low_0_bit(const int x) nogil:
+    """Get the position of the right-most 0 bit for an integer.
+
+    Examples:
+    >>> low_0_bit(0)
+    1
+    >>> low_0_bit(1)
+    2
+    >>> low_0_bit(2)
+    1
+    >>> low_0_bit(5)
+    2
+    >>> low_0_bit(7)
+    4
+
+    Parameters
+    ----------
+    x: int
+        An integer.
+
+    Returns
+    -------
+    position: int
+        Position of the right-most 0 bit.
+
+    """
+    cdef int z = x
+    cdef int i = 0
+    while True:
+        i += 1
+        if z % 2 == 0:
+            break
+        z = z // 2
+    return i
+
+
+@cython.boundscheck(False)
+@cython.wraparound(False)
+cdef int ibits(const int x, const int pos, const int length) nogil:
+    """Extract a sequence of bits from the bit representation of an integer.
+
+    Extract the sequence from position `pos` (inclusive) to ``pos + length``
+    (not inclusive), leftwise.
+
+    Examples:
+    >>> ibits(1, 0, 1)
+    1
+    >>> ibits(1, 1, 1)
+    0
+    >>> ibits(2, 0, 1)
+    0
+    >>> ibits(2, 0, 2)
+    2
+    >>> ibits(25, 1, 5)
+    12
+
+
+    Parameters
+    ----------
+    x: int
+        Integer to convert to bit representation.
+    pos: int
+        Starting position of sequence in bit representation of integer.
+    length: int
+        Length of sequence (number of bits).
+
+    Returns
+    -------
+    ibits: int
+        Integer value corresponding to bit sequence.
+
+    """
+    return (x >> pos) & ((1 << length) - 1)
+
+
+@cython.boundscheck(False)
+@cython.wraparound(False)
+cpdef void initialize_v(cnp.int_t[:, :] v, const int dim):
+    cdef int d, i, j, k, m, p, newv, pow2
+
+    if dim == 0:
+        return
+
+    # first row of v is all 1s
+    for i in range(MAXBIT):
+        v[0, i] = 1
+
+    # Remaining rows of v (row 2 through dim, indexed by [1:dim])
+    for d in range(1, dim):
+        p = poly[d]
+        m = bit_length(p) - 1
+
+        # First m elements of row d comes from vinit
+        for j in range(m):
+            v[d, j] = vinit[d][j]
+
+        # Fill in remaining elements of v as in Section 2 (top of pg. 90) of:
+        #
+        # P. Bratley and B. L. Fox. Algorithm 659: Implementing sobol's
+        # quasirandom sequence generator. ACM Trans.
+        # Math. Softw., 14(1):88-100, Mar. 1988.
+        #
+        for j in range(m, MAXBIT):
+            newv = v[d, j - m]
+            pow2 = 1
+            for k in range(m):
+                pow2 = pow2 << 1
+                if (p >> (m - 1 - k)) & 1:
+                    newv = newv ^ (pow2 * v[d, j - k - 1])
+            v[d, j] = newv
+
+    # Multiply each column of v by power of 2:
+    # v * [2^(maxbit-1), 2^(maxbit-2),..., 2, 1]
+    pow2 = 1
+    for d in range(MAXBIT):
+        for i in range(dim):
+            v[i, MAXBIT - 1 - d] *= pow2
+        pow2 = pow2 << 1
+
+
+@cython.boundscheck(False)
+@cython.wraparound(False)
+cpdef void _draw(const int n,
+                 const int num_gen,
+                 const int dim,
+                 cnp.int_t[:, :] sv,
+                 cnp.int_t[:] quasi,
+                 cnp.float_t[:, :] result) nogil:
+    cdef int i, j, l, qtmp
+    cdef int num_gen_loc = num_gen
+    for i in range(n):
+        l = low_0_bit(num_gen_loc)
+        for j in range(dim):
+            qtmp = quasi[j] ^ sv[j, l - 1]
+            quasi[j] = qtmp
+            result[i, j] = qtmp * RECIPD
+        num_gen_loc += 1
+
+
+@cython.boundscheck(False)
+@cython.wraparound(False)
+cpdef void _fast_forward(const int n,
+                         const int num_gen,
+                         const int dim,
+                         cnp.int_t[:, :] sv,
+                         cnp.int_t[:] quasi) nogil:
+    cdef int i, j, l
+    cdef int num_gen_loc = num_gen
+    for i in range(n):
+        l = low_0_bit(num_gen_loc)
+        for j in range(dim):
+            quasi[j] = quasi[j] ^ sv[j, l - 1]
+        num_gen_loc += 1
+
+
+@cython.boundscheck(False)
+@cython.wraparound(False)
+cdef int cdot_pow2(cnp.int_t[:] a) nogil:
+    cdef int i
+    cdef int size = a.shape[0]
+    cdef int z = 0
+    cdef int pow2 = 1
+    for i in range(size):
+        z += a[size - 1 - i] * pow2
+        pow2 *= 2
+    return z
+
+
+@cython.boundscheck(False)
+@cython.wraparound(False)
+cpdef void _cscramble(const int dim,
+                      cnp.int_t[:, :, :] ltm,
+                      cnp.int_t[:, :] sv) nogil:
+    cdef int d, i, j, k, l, lsm, lsmdp, p, t1, t2, vdj
+
+    # Set diagonals of maxbit x maxbit arrays to 1
+    for d in range(dim):
+        for i in range(MAXBIT):
+            ltm[d, i, i] = 1
+
+    for d in range(dim):
+        for j in range(MAXBIT):
+            vdj = sv[d, j]
+            l = 1
+            t2 = 0
+            for p in range(MAXBIT - 1, -1, -1):
+                lsmdp = cdot_pow2(ltm[d, p, :])
+                t1 = 0
+                for k in range(MAXBIT):
+                    t1 += ibits(lsmdp, k, 1) * ibits(vdj, k, 1)
+                t1 = t1 % 2
+                t2 = t2 + t1 * l
+                l = 2 * l
+            sv[d, j] = t2
+
+
+@cython.boundscheck(False)
+@cython.wraparound(False)
+cpdef void _fill_p_cumulative(cnp.float_t[:] p,
+                              cnp.float_t[:] p_cumulative) nogil:
+    cdef int i
+    cdef int len_p = p.shape[0]
+    cdef float tot = 0
+    cdef float t
+    for i in range(len_p):
+        t = tot + p[i]
+        p_cumulative[i] = t
+        tot = t
+
+
+@cython.boundscheck(False)
+@cython.wraparound(False)
+cpdef void _categorize(cnp.float_t[:] draws,
+                       cnp.float_t[:] p_cumulative,
+                       cnp.int_t[:] result) nogil:
+    cdef int i
+    cdef int n_p = p_cumulative.shape[0]
+    for i in range(draws.shape[0]):
+        j = _find_index(p_cumulative, n_p, draws[i])
+        result[j] = result[j] + 1
+
+
+@cython.boundscheck(False)
+@cython.wraparound(False)
+cdef int _find_index(cnp.float_t[:] p_cumulative,
+                     const int size,
+                     const float value) nogil:
+    cdef int l = 0
+    cdef int r = size - 1
+    cdef int m
+    while r > l:
+        m = (l + r) // 2
+        if value > p_cumulative[m]:
+            l = m + 1
+        else:
+            r = m
+    return r
+
+
+def _test_find_index(p_cumulative, size, value):
+    # type: (np.ndarray, int, float) -> int
+    """Wrapper for testing in python"""
+    return _find_index(p_cumulative, size, value)
diff --git a/scipy/stats/_sobol_direction_numbers.npz b/scipy/stats/_sobol_direction_numbers.npz
new file mode 100644
index 000000000000..da96d4517d80
Binary files /dev/null and b/scipy/stats/_sobol_direction_numbers.npz differ
diff --git a/scipy/stats/qmc.py b/scipy/stats/qmc.py
new file mode 100644
index 000000000000..27a0082ee862
--- /dev/null
+++ b/scipy/stats/qmc.py
@@ -0,0 +1,229 @@
+# -*- coding: utf-8 -*-
+r"""
+====================================================
+Quasi-Monte Carlo submodule (:mod:`scipy.stats.qmc`)
+====================================================
+
+.. currentmodule:: scipy.stats.qmc
+
+This module provides Quasi-Monte Carlo generators and associated helper
+functions.
+
+
+Quasi-Monte Carlo
+=================
+
+Engines
+-------
+
+.. autosummary::
+   :toctree: generated/
+
+   QMCEngine
+   Sobol
+   Halton
+   LatinHypercube
+   OrthogonalLatinHypercube
+   OptimalDesign
+   NormalQMC
+   MultivariateNormalQMC
+
+Helpers
+-------
+
+.. autosummary::
+   :toctree: generated/
+
+   discrepancy
+   scale
+   multinomial_qmc
+
+
+Introduction to Quasi-Monte Carlo
+=================================
+
+Quasi-Monte Carlo (QMC) methods [1]_, [2]_, [3]_ provide an
+:math:`n \times d` array of numbers in :math:`[0,1]`. They can be used in
+place of n points from the :math:`U[0,1]^{d}` distribution. Compared to
+random points, QMC points are designed to have fewer gaps and clumps. This is
+quantified by discrepancy measures [4]_. From the Koksma-Hlawka
+inequality [5]_ we know that low discrepancy reduces a bound on
+integration error. Averaging a function :math:`f` over :math:`n` QMC points
+can achieve an integration error close to :math:`O(n^{-1})` for well
+behaved functions [2]_.
+
+Most QMC constructions are designed for special values of :math:`n`
+such as powers of 2 or large primes. Changing the sample
+size by even one can degrade their performance, even their
+rate of convergence [6]_. For instance :math:`n=100` points may give less
+accuracy than :math:`n=64` if the method was designed for :math:`n=2^m`.
+
+Some QMC constructions are extensible in :math:`n`: we can find
+another special sample size :math:`n' > n` and often an infinite
+sequence of increasing special sample sizes. Some QMC
+constructions are extensible in :math:`d`: we can increase the dimension,
+possibly to some upper bound, and typically without requiring
+special values of :math:`d`. Some QMC methods are extensible in
+both :math:`n` and :math:`d`.
+
+QMC points are deterministic. That makes it hard to estimate
+the accuracy of averages. Randomized QMC (RQMC) [7]_
+points are constructed so that each point is individually :math:`U[0,1]^{d}`
+while collectively the :math:`n` points retain their low discrepancy.
+One can make :math:`R` independent replications of RQMC points to
+see how stable a computation is. From :math:`R` independent values,
+a t-test (or bootstrap t-test [8]_) then gives approximate confidence
+intervals on the mean value.  Some RQMC methods produce a
+root mean squared error that is actually :math:`o(1/n)` and smaller than
+the rate seen in unrandomized QMC.  An intuitive explanation is
+that the error is a sum of many small ones and random errors
+cancel in a way that deterministic ones do not. RQMC also
+has advantages on integrands that are singular or, for other
+reasons, fail to be Riemann integrable.
+
+(R)QMC cannot beat Bahkvalov's curse of dimension (see [9]_). For
+any random or deterministic method, there are worst case functions
+that will give it poor performance in high dimensions. A worst
+case function for QMC might be 0 at all n points but very
+large elsewhere. Worst case analyses get very pessimistic
+in high dimensions. (R)QMC can bring a great improvement over
+MC when the functions on which it is used are not worst case.
+For instance (R)QMC can be especially effective on integrands
+that are well approximated by sums of functions of
+some small number of their input variables at a time [10]_, [11]_.
+That property is often a surprising finding about those functions.
+
+Scrambled nets are a kind of RQMC that have some valuable robustness
+properties [12]_. If the integrand is square integrable, they give variance
+:math:`var_{SNET} = o(1/n)`. There is a finite upper bound on
+:math:`var_{SNET} / var_{MC}` that holds simultaneously for every square
+integrable integrand. Scrambled nets satisfy a strong law of large numbers
+for :math:`f` in :math:`L^p` when :math:`p>1`. In some
+special cases there is a central limit theorem [13]_. For smooth enough
+integrands they can achieve RMSE nearly :math:`O(n^{-3})`. See [12]_
+for references about these properties.
+
+The main kinds of QMC methods are lattice rules [14]_ and digital
+nets and sequences [2]_, [15]_. The theories meet up in polynomial
+lattice rules [16]_ which can produce digital nets. Lattice rules
+require some form of search for good constructions. For digital
+nets there are widely used default constructions.
+
+The most widely used QMC methods are Sobol' sequences [17]_.
+These are digital nets. They are extensible in both :math:`n` and :math:`d`.
+They can be scrambled.  The special sample sizes are powers
+of 2. Another popular method are Halton sequences [18]_.
+The constructions resemble those of digital nets. The earlier
+dimensions have much better equidistribution properties than
+later ones. There are essentially no special sample sizes.
+They are not thought to be as accurate as Sobol' sequences.
+They can be scrambled. The nets of Faure [19]_ are also widely
+used. All dimensions are equally good, but the special sample
+sizes grow rapidly with dimension :math:`d`. They can be scrambled.
+The nets of Niederreiter and Xing [20]_ have the best asymptotic
+properties but have not shown good empirical performance [21]_.
+
+Higher order digital nets are formed by a digit interleaving process
+in the digits of the constructed points. They can achieve higher
+levels of asymptotic accuracy given higher smoothness conditions on :math:`f`
+and they can be scrambled [22]_. There is little or no empirical work
+showing the improved rate to be attained.
+
+Using QMC is like using the entire period of a small random
+number generator. The constructions are similar and so
+therefore are the computational costs [23]_.
+
+(R)QMC is sometimes improved by passing the points through
+a baker's transformation (tent function) prior to using them.
+That function has the form :math:`1-2|x-1/2|`.  As :math:`x` goes from 0 to
+1, this function goes from 0 to 1 and then back.  It is very
+useful to produce a periodic function for lattice rules [14]_,
+and sometimes it improves the convergence rate [24]_.
+
+It is not straightforward to apply QMC methods to Markov
+chain Monte Carlo (MCMC).  We can think of MCMC as using
+:math:`n=1` point in :math:`[0,1]^{d}` for very large :math:`d`, with
+ergodic results corresponding to :math:`d \to \infty`.  One proposal is
+in [25]_ and under strong conditions an improved rate of convergence
+has been shown [26]_.
+
+Returning to Sobol' points: there are many versions depending
+on what are called direction numbers. Those are the result of
+searches and are tabulated. A very widely used set of direction
+numbers come from [27]_. It is extensible in dimension up to
+:math:`d=21201`.
+
+References
+----------
+.. [1] Owen, Art B. "Monte Carlo Book: the Quasi-Monte Carlo parts." (2019).
+.. [2] Niederreiter, Harald. Random number generation and quasi-Monte Carlo
+   methods. Society for Industrial and Applied Mathematics, 1992.
+.. [3] Dick, Josef, Frances Y. Kuo, and Ian H. Sloan. "High-dimensional
+   integration: the quasi-Monte Carlo way." Acta Numerica 22 (2013): 133.
+.. [4] Aho, A. V., C. Aistleitner, T. Anderson, K. Appel, V. Arnol'd, N.
+   Aronszajn, D. Asotsky et al. "W. Chen et al.(eds.), A Panorama of
+   Discrepancy Theory (2014): 679. Sringer International Publishing,
+   Switzerland.
+.. [5] Hickernell, Fred J. "Koksma-Hlawka Inequality." Wiley StatsRef:
+   Statistics Reference Online (2014).
+.. [6] Owen, Art B. "On dropping the first Sobol' point." arXiv preprint
+   arXiv:2008.08051 (2020).
+.. [7] L'Ecuyer, Pierre, and Christiane Lemieux. "Recent advances in randomized
+   quasi-Monte Carlo methods." In Modeling uncertainty, pp. 419-474. Springer,
+   New York, NY, 2002.
+.. [8] DiCiccio, Thomas J., and Bradley Efron. "Bootstrap confidence
+   intervals." Statistical science (1996): 189-212.
+.. [9] Dimov, Ivan T. Monte Carlo methods for applied scientists. World
+   Scientific, 2008.
+.. [10] Caflisch, Russel E., William J. Morokoff, and Art B. Owen. Valuation of
+   mortgage backed securities using Brownian bridges to reduce effective
+   dimension. Journal of Computational Finance, (1997): 1, no. 1 27-46.
+.. [11] Sloan, Ian H., and Henryk Wozniakowski. "When are quasi-Monte Carlo
+   algorithms efficient for high dimensional integrals?." Journal of Complexity
+   14, no. 1 (1998): 1-33.
+.. [12] Owen, Art B., and Daniel Rudolf "A strong law of large numbers for
+   scrambled net integration." SIAM Review, to appear.
+.. [13] Loh, Wei-Liem. "On the asymptotic distribution of scrambled net
+   quadrature." The Annals of Statistics 31, no. 4 (2003): 1282-1324.
+.. [14] Sloan, Ian H. and S. Joe. Lattice methods for multiple integration.
+   Oxford University Press, 1994.
+.. [15] Dick, Josef, and Friedrich Pillichshammer. Digital nets and sequences:
+   discrepancy theory and quasi-Monte Carlo integration. Cambridge University
+   Press, 2010.
+.. [16] Dick, Josef, F. Kuo, Friedrich Pillichshammer, and I. Sloan.
+   "Construction algorithms for polynomial lattice rules for multivariate
+   integration." Mathematics of computation 74, no. 252 (2005): 1895-1921.
+.. [17] Sobol', Il'ya Meerovich. "On the distribution of points in a cube and
+   the approximate evaluation of integrals." Zhurnal Vychislitel'noi Matematiki
+   i Matematicheskoi Fiziki 7, no. 4 (1967): 784-802.
+.. [18] Halton, John H. "On the efficiency of certain quasi-random sequences of
+   points in evaluating multi-dimensional integrals." Numerische Mathematik 2,
+   no. 1 (1960): 84-90.
+.. [19] Faure, Henri. "Discrepance de suites associees a un systeme de
+   numeration (en dimension s)." Acta arithmetica 41, no. 4 (1982): 337-351.
+.. [20] Niederreiter, Harold, and Chaoping Xing. "Low-discrepancy sequences and
+   global function fields with many rational places." Finite Fields and their
+   applications 2, no. 3 (1996): 241-273.
+.. [21] Hong, Hee Sun, and Fred J. Hickernell. "Algorithm 823: Implementing
+   scrambled digital sequences." ACM Transactions on Mathematical Software
+   (TOMS) 29, no. 2 (2003): 95-109.
+.. [22] Dick, Josef. "Higher order scrambled digital nets achieve the optimal
+   rate of the root mean square error for smooth integrands." The Annals of
+   Statistics 39, no. 3 (2011): 1372-1398.
+.. [23] Niederreiter, Harald. "Multidimensional numerical integration using
+   pseudorandom numbers." In Stochastic Programming 84 Part I, pp. 17-38.
+   Springer, Berlin, Heidelberg, 1986.
+.. [24] Hickernell, Fred J. "Obtaining O (N-2+e) Convergence for Lattice
+   Quadrature Rules." In Monte Carlo and Quasi-Monte Carlo Methods 2000,
+   pp. 274-289. Springer, Berlin, Heidelberg, 2002.
+.. [25] Owen, Art B., and Seth D. Tribble. "A quasi-Monte Carlo Metropolis
+   algorithm." Proceedings of the National Academy of Sciences 102, no. 25
+   (2005): 8844-8849.
+.. [26] Chen, Su. "Consistency and convergence rate of Markov chain quasi Monte
+   Carlo with examples." PhD diss., Stanford University, 2011.
+.. [27] Joe, Stephen, and Frances Y. Kuo. "Constructing Sobol sequences with
+   better two-dimensional projections." SIAM Journal on Scientific Computing
+   30, no. 5 (2008): 2635-2654.
+
+"""
+from ._qmc import *
diff --git a/scipy/stats/setup.py b/scipy/stats/setup.py
index 94dc721bbc73..fcffed566a92 100644
--- a/scipy/stats/setup.py
+++ b/scipy/stats/setup.py
@@ -25,9 +25,15 @@ def configuration(parent_package='',top_path=None):
 
     # add mvn module
     config.add_extension('mvn',
-        sources=['mvn.pyf','mvndst.f'],
+        sources=['mvn.pyf', 'mvndst.f'],
     )
 
+    # add _sobol module
+    config.add_extension('_sobol',
+        sources=['_sobol.c', ],
+    )
+    config.add_data_files('_sobol_direction_numbers.npz')
+
     return config
 
 
diff --git a/scipy/stats/tests/test_qmc.py b/scipy/stats/tests/test_qmc.py
new file mode 100644
index 000000000000..2121aa11c86c
--- /dev/null
+++ b/scipy/stats/tests/test_qmc.py
@@ -0,0 +1,624 @@
+import copy
+
+import pytest
+import numpy as np
+from numpy.testing import (assert_allclose, assert_almost_equal, assert_,
+                           assert_equal, assert_array_almost_equal,
+                           assert_array_equal)
+from scipy.stats import shapiro
+
+from scipy.stats._sobol import _test_find_index
+from scipy.stats import qmc
+from scipy.stats._qmc import (van_der_corput, n_primes, primes_from_2_to,
+                              _perturb_discrepancy, _update_discrepancy)
+
+
+class TestUtils(object):
+    def test_scale(self):
+        space = [[0, 0], [1, 1], [0.5, 0.5]]
+        corners = np.array([[-2, 0], [6, 5]])
+        out = [[-2, 0], [6, 5], [2, 2.5]]
+
+        scaled_space = qmc.scale(space, bounds=corners)
+
+        assert_allclose(scaled_space, out)
+
+        scaled_back_space = qmc.scale(scaled_space, bounds=corners,
+                                      reverse=True)
+        assert_allclose(scaled_back_space, space)
+
+    def test_discrepancy(self):
+        space_1 = np.array([[1, 3], [2, 6], [3, 2], [4, 5], [5, 1], [6, 4]])
+        space_1 = (2.0 * space_1 - 1.0) / (2.0 * 6.0)
+        space_2 = np.array([[1, 5], [2, 4], [3, 3], [4, 2], [5, 1], [6, 6]])
+        space_2 = (2.0 * space_2 - 1.0) / (2.0 * 6.0)
+
+        # From Fang et al. Design and modeling for computer experiments, 2006
+        assert_allclose(qmc.discrepancy(space_1), 0.0081, atol=1e-4)
+        assert_allclose(qmc.discrepancy(space_2), 0.0105, atol=1e-4)
+
+        # From Zhou Y.-D. et al. Mixture discrepancy for quasi-random point
+        # sets. Journal of Complexity, 29 (3-4), pp. 283-301, 2013.
+        sample = np.array([[2, 1, 1, 2, 2, 2],
+                           [1, 2, 2, 2, 2, 2],
+                           [2, 1, 1, 1, 1, 1],
+                           [1, 1, 1, 1, 2, 2],
+                           [1, 2, 2, 2, 1, 1],
+                           [2, 2, 2, 2, 1, 1],
+                           [2, 2, 2, 1, 2, 2]])
+        sample = (2.0 * sample - 1.0) / (2.0 * 2.0)
+
+        assert_allclose(qmc.discrepancy(sample, method='MD'), 2.5000,
+                        atol=1e-4)
+        assert_allclose(qmc.discrepancy(sample, method='WD'), 1.3680,
+                        atol=1e-4)
+        assert_allclose(qmc.discrepancy(sample, method='CD'), 0.3172,
+                        atol=1e-4)
+        assert_allclose(qmc.discrepancy(sample, method='star'), 0.037451,
+                        atol=1e-4)
+
+        with pytest.raises(ValueError, match=r"toto is not a valid method."):
+            qmc.discrepancy(sample, method='toto')
+
+    def test_update_discrepancy(self):
+        space_1 = np.array([[1, 3], [2, 6], [3, 2], [4, 5], [5, 1], [6, 4]])
+        space_1 = (2.0 * space_1 - 1.0) / (2.0 * 6.0)
+
+        disc_init = qmc.discrepancy(space_1[:-1], iterative=True)
+        disc_iter = _update_discrepancy(space_1[-1], space_1[:-1],
+                                        disc_init)
+
+        assert_allclose(disc_iter, 0.0081, atol=1e-4)
+
+    def test_perm_discrepancy(self):
+        doe_init = np.array([[1, 3], [2, 6], [3, 2], [4, 5], [5, 1], [6, 4]])
+        doe_init = (2.0 * doe_init - 1.0) / (2.0 * 6.0)
+
+        disc_init = qmc.discrepancy(doe_init)
+
+        row_1, row_2, col = 5, 2, 1
+
+        doe = copy.deepcopy(doe_init)
+        doe[row_1, col], doe[row_2, col] = doe[row_2, col], doe[row_1, col]
+
+        disc_valid = qmc.discrepancy(doe)
+        disc_perm = _perturb_discrepancy(doe_init, row_1, row_2, col,
+                                         disc_init)
+
+        assert_allclose(disc_valid, disc_perm)
+
+    def test_n_primes(self):
+        primes = n_primes(10)
+        assert primes[-1] == 29
+
+        primes = n_primes(168)
+        assert primes[-1] == 997
+
+        primes = n_primes(350)
+        assert primes[-1] == 2357
+
+    def test_primes(self):
+        primes = primes_from_2_to(50)
+        out = [2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47]
+        assert_allclose(primes, out)
+
+
+class TestQMC(object):
+    def test_van_der_corput(self):
+        sample = van_der_corput(10)
+        out = [0., 0.5, 0.25, 0.75, 0.125, 0.625, 0.375, 0.875, 0.0625, 0.5625]
+        assert_almost_equal(sample, out)
+
+        sample = van_der_corput(7, start_index=3)
+        assert_almost_equal(sample, out[3:])
+
+    def test_van_der_corput_scramble(self):
+        out = van_der_corput(10, scramble=True, seed=123456)
+
+        sample = van_der_corput(7, start_index=3, scramble=True,
+                                seed=123456)
+        assert_almost_equal(sample, out[3:])
+
+    def test_halton(self):
+        engine = qmc.Halton(d=2, scramble=False)
+        sample = engine.random(n=3)
+
+        out = np.array([[0, 0], [1 / 2, 1 / 3], [1 / 4, 2 / 3], [3 / 4, 1 / 9],
+                        [1 / 8, 4 / 9], [5 / 8, 7 / 9]])
+        assert_almost_equal(sample, out[:3], decimal=1)
+
+        assert engine.num_generated == 3
+
+        # continuing
+        sample = engine.random(n=3)
+        assert_almost_equal(sample, out[3:], decimal=1)
+
+        # reset
+        engine.reset()
+        sample = engine.random(n=6)
+        assert_almost_equal(sample, out, decimal=1)
+
+        # continuing with fast_forward
+        engine = qmc.Halton(d=2, scramble=False)
+        engine.fast_forward(2)
+        sample = engine.random(n=4)
+
+        assert_almost_equal(sample, out[2:], decimal=1)
+
+    def test_halton_scramble(self):
+        engine = qmc.Halton(d=2, scramble=True, seed=123456)
+        out = engine.random(n=10)
+
+        engine = qmc.Halton(d=2, scramble=True, seed=123456)
+        sample = engine.random(n=3)
+        assert_almost_equal(sample, out[:3], decimal=1)
+
+        # continuing
+        sample = engine.random(n=7)
+        assert_almost_equal(sample, out[3:], decimal=1)
+
+        # reset
+        engine.reset()
+        sample = engine.random(n=10)
+        assert_almost_equal(sample, out, decimal=1)
+
+        # continuing with fast_forward
+        engine = qmc.Halton(d=2, scramble=True, seed=123456)
+        engine.fast_forward(2)
+        sample = engine.random(n=8)
+
+        assert_almost_equal(sample, out[2:], decimal=1)
+
+
+class TestLHS(object):
+    def test_lhs(self):
+        np.random.seed(123456)
+
+        corners = np.array([[0, 2], [10, 5]])
+
+        lhs = qmc.LatinHypercube(d=2)
+        sample = lhs.random(n=5)
+        sample = qmc.scale(sample, bounds=corners)
+        out = np.array([[5.7, 3.2], [5.5, 3.9], [5.2, 3.6],
+                        [5.1, 3.3], [5.8, 4.1]])
+        assert_almost_equal(sample, out, decimal=1)
+
+        lhs = qmc.LatinHypercube(d=2, centered=True)
+        sample = lhs.random(n=5)
+        out = np.array([[0.1, 0.5], [0.3, 0.1], [0.7, 0.1],
+                        [0.1, 0.1], [0.3, 0.7]])
+        assert_almost_equal(sample, out, decimal=1)
+
+    def test_orthogonal_lhs(self):
+        np.random.seed(123456)
+
+        corners = np.array([[0, 2], [10, 5]])
+
+        olhs = qmc.OrthogonalLatinHypercube(d=2)
+        sample = olhs.random(n=5)
+        sample = qmc.scale(sample, bounds=corners)
+        out = np.array([[3.933, 2.670], [7.794, 4.031], [4.520, 2.129],
+                        [0.253, 4.976], [8.753, 3.249]])
+        assert_almost_equal(sample, out, decimal=1)
+
+        # Checking independency of the random numbers generated
+        n_samples = 500
+        sample = olhs.random(n=n_samples)
+        min_b = 50  # number of bins
+        bins = np.linspace(0, 1, min(min_b, n_samples) + 1)
+        hist = np.histogram(sample[:, 0], bins=bins)
+        out = np.array([n_samples / min_b] * min_b)
+        assert_equal(hist[0], out)
+
+        hist = np.histogram(sample[:, 1], bins=bins)
+        assert_equal(hist[0], out)
+
+    def test_optimal_design(self):
+        seed = 123456
+
+        olhs = qmc.OrthogonalLatinHypercube(d=2, seed=seed)
+        sample_ref = olhs.random(n=20)
+        disc_ref = qmc.discrepancy(sample_ref)
+
+        optimal_1 = qmc.OptimalDesign(d=2, start_design=sample_ref,
+                                      seed=seed)
+        sample_1 = optimal_1.random(n=20)
+        disc_1 = qmc.discrepancy(sample_1)
+
+        assert disc_1 < disc_ref
+
+        optimal_ = qmc.OptimalDesign(d=2, seed=seed)
+        sample_ = optimal_.random(n=20)
+        disc_ = qmc.discrepancy(sample_)
+
+        assert disc_ < disc_ref
+
+        optimal_2 = qmc.OptimalDesign(d=2, start_design=sample_ref, niter=2,
+                                      seed=seed)
+        sample_2 = optimal_2.random(n=20)
+        disc_2 = qmc.discrepancy(sample_2)
+        assert disc_2 < disc_1
+
+        # resample is like doing another iteration
+        sample_1 = optimal_1.random(n=20)
+        assert_allclose(sample_1, sample_2)
+
+        optimal_3 = qmc.OptimalDesign(d=2, start_design=sample_ref,
+                                      force=True)
+        sample_3 = optimal_3.random(n=20)
+        disc_3 = qmc.discrepancy(sample_3)
+        assert disc_3 < disc_ref
+
+        # no optimization
+        optimal_4 = qmc.OptimalDesign(d=2, start_design=sample_ref,
+                                      optimization=False, niter=100, seed=seed)
+        sample_4 = optimal_4.random(n=20)
+        disc_4 = qmc.discrepancy(sample_4)
+
+        assert disc_4 < disc_ref
+
+
+class TestMultinomialQMC:
+    def test_MultinomialNegativePs(self):
+        p = np.array([0.12, 0.26, -0.05, 0.35, 0.22])
+        with pytest.raises(ValueError, match=r"Elements of pvals must "
+                                             r"be non-negative."):
+            qmc.multinomial_qmc(10, p)
+
+    def test_MultinomialSumOfPTooLarge(self):
+        p = np.array([0.12, 0.26, 0.1, 0.35, 0.22])
+        with pytest.raises(ValueError, match=r"Elements of pvals must sum "
+                                             r"to 1."):
+            qmc.multinomial_qmc(10, p)
+
+    @pytest.mark.filterwarnings('ignore::UserWarning')
+    def test_MultinomialBasicDraw(self):
+        p = np.array([0.12, 0.26, 0.05, 0.35, 0.22])
+        expected = np.array([12, 25, 6, 35, 22])
+        assert_array_equal(qmc.multinomial_qmc(100, p, seed=12345), expected)
+
+    def test_MultinomialDistribution(self):
+        p = np.array([0.12, 0.26, 0.05, 0.35, 0.22])
+        draws = qmc.multinomial_qmc(8192, p, seed=12345)
+        assert_array_almost_equal(draws / np.sum(draws), p, decimal=4)
+
+    def test_FindIndex(self):
+        p_cumulative = np.array([0.1, 0.4, 0.45, 0.6, 0.75, 0.9, 0.99, 1.0])
+        size = len(p_cumulative)
+        assert_equal(_test_find_index(p_cumulative, size, 0.0), 0)
+        assert_equal(_test_find_index(p_cumulative, size, 0.4), 2)
+        assert_equal(_test_find_index(p_cumulative, size, 0.44999), 2)
+        assert_equal(_test_find_index(p_cumulative, size, 0.45001), 3)
+        assert_equal(_test_find_index(p_cumulative, size, 1.0), size - 1)
+
+    @pytest.mark.filterwarnings('ignore::UserWarning')
+    def test_other_engine(self):
+        # same as test_MultinomialBasicDraw with different engine
+        p = np.array([0.12, 0.26, 0.05, 0.35, 0.22])
+        expected = np.array([12, 25, 6, 35, 22])
+        engine = qmc.Sobol(1, scramble=True, seed=12345)
+        assert_array_equal(qmc.multinomial_qmc(100, p, engine=engine,
+                                               seed=12345),
+                           expected)
+
+
+class TestNormalQMC:
+    def test_NormalQMC(self):
+        # d = 1
+        engine = qmc.NormalQMC(d=1)
+        samples = engine.random()
+        assert_equal(samples.shape, (1, 1))
+        samples = engine.random(n=5)
+        assert_equal(samples.shape, (5, 1))
+        # d = 2
+        engine = qmc.NormalQMC(d=2)
+        samples = engine.random()
+        assert_equal(samples.shape, (1, 2))
+        samples = engine.random(n=5)
+        assert_equal(samples.shape, (5, 2))
+
+    def test_NormalQMCInvTransform(self):
+        # d = 1
+        engine = qmc.NormalQMC(d=1, inv_transform=True)
+        samples = engine.random()
+        assert_equal(samples.shape, (1, 1))
+        samples = engine.random(n=5)
+        assert_equal(samples.shape, (5, 1))
+        # d = 2
+        engine = qmc.NormalQMC(d=2, inv_transform=True)
+        samples = engine.random()
+        assert_equal(samples.shape, (1, 2))
+        samples = engine.random(n=5)
+        assert_equal(samples.shape, (5, 2))
+
+    def test_other_engine(self):
+        engine = qmc.NormalQMC(d=2, engine=qmc.Sobol(d=2, scramble=False),
+                               inv_transform=True)
+        samples = engine.random()
+        assert_equal(samples.shape, (1, 2))
+
+    def test_NormalQMCSeeded(self):
+        # test even dimension
+        engine = qmc.NormalQMC(d=2, inv_transform=False, seed=12345)
+        samples = engine.random(n=2)
+        samples_expected = np.array(
+            [[-0.943472, 0.405116], [-0.63099602, -1.32950772]]
+        )
+        assert_array_almost_equal(samples, samples_expected)
+        # test odd dimension
+        engine = qmc.NormalQMC(d=3, inv_transform=False, seed=12345)
+        samples = engine.random(n=2)
+        samples_expected = np.array(
+            [
+                [-0.943472, 0.405116, 0.268828],
+                [1.83169884, -1.40473647, 0.24334828],
+            ]
+        )
+        assert_array_almost_equal(samples, samples_expected)
+
+    def test_NormalQMCSeededInvTransform(self):
+        # test even dimension
+        engine = qmc.NormalQMC(d=2, seed=12345, inv_transform=True)
+        samples = engine.random(n=2)
+        samples_expected = np.array(
+            [[0.228309, -0.162516], [-0.41622922, 0.46622792]]
+        )
+        assert_array_almost_equal(samples, samples_expected)
+        # test odd dimension
+        engine = qmc.NormalQMC(d=3, seed=12345, inv_transform=True)
+        samples = engine.random(n=2)
+        samples_expected = np.array(
+            [
+                [0.228309, -0.162516, 0.167352],
+                [-1.40525266, 1.37652443, -0.8519666],
+            ]
+        )
+        assert_array_almost_equal(samples, samples_expected)
+
+    def test_NormalQMCShapiro(self):
+        engine = qmc.NormalQMC(d=2, seed=12345)
+        samples = engine.random(n=256)
+        assert_(all(np.abs(samples.mean(axis=0)) < 1e-2))
+        assert_(all(np.abs(samples.std(axis=0) - 1) < 1e-2))
+        # perform Shapiro-Wilk test for normality
+        for i in (0, 1):
+            _, pval = shapiro(samples[:, i])
+            assert_(pval > 0.9)
+        # make sure samples are uncorrelated
+        cov = np.cov(samples.transpose())
+        assert_(np.abs(cov[0, 1]) < 1e-2)
+
+    def test_NormalQMCShapiroInvTransform(self):
+        engine = qmc.NormalQMC(d=2, seed=12345, inv_transform=True)
+        samples = engine.random(n=256)
+        assert_(all(np.abs(samples.mean(axis=0)) < 1e-2))
+        assert_(all(np.abs(samples.std(axis=0) - 1) < 1e-2))
+        # perform Shapiro-Wilk test for normality
+        for i in (0, 1):
+            _, pval = shapiro(samples[:, i])
+            assert_(pval > 0.9)
+        # make sure samples are uncorrelated
+        cov = np.cov(samples.transpose())
+        assert_(np.abs(cov[0, 1]) < 1e-2)
+
+
+class TestMultivariateNormalQMC:
+    def test_MultivariateNormalQMCNonPSD(self):
+        # try with non-psd, non-pd cov and expect an assertion error
+        with pytest.raises(ValueError, match=r"Covariance matrix not PSD."):
+            qmc.MultivariateNormalQMC([0, 0], [[1, 2], [2, 1]])
+
+    def test_MultivariateNormalQMCNonPD(self):
+        # try with non-pd but psd cov; should work
+        engine = qmc.MultivariateNormalQMC(
+            [0, 0, 0], [[1, 0, 1], [0, 1, 1], [1, 1, 2]]
+        )
+        assert_(engine._corr_matrix is not None)
+
+    def test_MultivariateNormalQMCSymmetric(self):
+        # try with non-symmetric cov and expect an error
+        with pytest.raises(ValueError, match=r"Covariance matrix is not "
+                                             r"symmetric."):
+            qmc.MultivariateNormalQMC([0, 0], [[1, 0], [2, 1]])
+
+    def test_MultivariateNormalQMCDim(self):
+        # incompatible dimension of mean/cov
+        with pytest.raises(ValueError, match=r"Dimension mismatch between "
+                                             r"mean and covariance."):
+            qmc.MultivariateNormalQMC([0], [[1, 0], [0, 1]])
+
+    def test_MultivariateNormalQMC(self):
+        # d = 1 scalar
+        engine = qmc.MultivariateNormalQMC(mean=0, cov=5)
+        samples = engine.random()
+        assert_equal(samples.shape, (1, 1))
+        samples = engine.random(n=5)
+        assert_equal(samples.shape, (5, 1))
+
+        # d = 2 list
+        engine = qmc.MultivariateNormalQMC(mean=[0, 1], cov=[[1, 0], [0, 1]])
+        samples = engine.random()
+        assert_equal(samples.shape, (1, 2))
+        samples = engine.random(n=5)
+        assert_equal(samples.shape, (5, 2))
+
+        # d = 3 np.array
+        mean = np.array([0, 1, 2])
+        cov = np.array([[1, 0, 0], [0, 1, 0], [0, 0, 1]])
+        engine = qmc.MultivariateNormalQMC(mean, cov)
+        samples = engine.random()
+        assert_equal(samples.shape, (1, 3))
+        samples = engine.random(n=5)
+        assert_equal(samples.shape, (5, 3))
+
+    def test_MultivariateNormalQMCInvTransform(self):
+        # d = 1 scalar
+        engine = qmc.MultivariateNormalQMC(mean=0, cov=5, inv_transform=True)
+        samples = engine.random()
+        assert_equal(samples.shape, (1, 1))
+        samples = engine.random(n=5)
+        assert_equal(samples.shape, (5, 1))
+
+        # d = 2 list
+        engine = qmc.MultivariateNormalQMC(
+            mean=[0, 1], cov=[[1, 0], [0, 1]], inv_transform=True
+        )
+        samples = engine.random()
+        assert_equal(samples.shape, (1, 2))
+        samples = engine.random(n=5)
+        assert_equal(samples.shape, (5, 2))
+
+        # d = 3 np.array
+        mean = np.array([0, 1, 2])
+        cov = np.array([[1, 0, 0], [0, 1, 0], [0, 0, 1]])
+        engine = qmc.MultivariateNormalQMC(mean, cov, inv_transform=True)
+        samples = engine.random()
+        assert_equal(samples.shape, (1, 3))
+        samples = engine.random(n=5)
+        assert_equal(samples.shape, (5, 3))
+
+    def test_MultivariateNormalQMCSeeded(self):
+        # test even dimension
+        np.random.seed(54321)
+        a = np.random.randn(2, 2)
+        A = a @ a.transpose() + np.diag(np.random.rand(2))
+        engine = qmc.MultivariateNormalQMC(np.array([0, 0]), A,
+                                           inv_transform=False, seed=12345)
+        samples = engine.random(n=2)
+        samples_expected = np.array(
+            [[-1.010703, -0.324223], [-0.67595995, -2.27437872]]
+        )
+        assert_array_almost_equal(samples, samples_expected)
+
+        # test odd dimension
+        np.random.seed(54321)
+        a = np.random.randn(3, 3)
+        A = a @ a.transpose() + np.diag(np.random.rand(3))
+        engine = qmc.MultivariateNormalQMC(np.array([0, 0, 0]), A,
+                                           inv_transform=False, seed=12345)
+        samples = engine.random(n=2)
+        samples_expected = np.array(
+            [
+                [-1.056834, 2.493251, 0.114556],
+                [2.05178452, -6.35744194, 0.67944512],
+            ]
+        )
+        assert_array_almost_equal(samples, samples_expected)
+
+    def test_MultivariateNormalQMCSeededInvTransform(self):
+        # test even dimension
+        np.random.seed(54321)
+        a = np.random.randn(2, 2)
+        A = a @ a.transpose() + np.diag(np.random.rand(2))
+        engine = qmc.MultivariateNormalQMC(
+            np.array([0, 0]), A, seed=12345, inv_transform=True
+        )
+        samples = engine.random(n=2)
+        samples_expected = np.array(
+            [[0.244578, -0.004441], [-0.44588916, 0.22657776]]
+        )
+        assert_array_almost_equal(samples, samples_expected)
+
+        # test odd dimension
+        np.random.seed(54321)
+        a = np.random.randn(3, 3)
+        A = a @ a.transpose() + np.diag(np.random.rand(3))
+        engine = qmc.MultivariateNormalQMC(
+            np.array([0, 0, 0]), A, seed=12345, inv_transform=True
+        )
+        samples = engine.random(n=2)
+        samples_expected = np.array(
+            [
+                [0.255741, -0.761559, 0.234236],
+                [-1.5740992, 5.61057598, -1.28218525],
+            ]
+        )
+        assert_array_almost_equal(samples, samples_expected)
+
+    def test_MultivariateNormalQMCShapiro(self):
+        # test the standard case
+        engine = qmc.MultivariateNormalQMC(
+            mean=[0, 0], cov=[[1, 0], [0, 1]], seed=12345
+        )
+        samples = engine.random(n=256)
+        assert_(all(np.abs(samples.mean(axis=0)) < 1e-2))
+        assert_(all(np.abs(samples.std(axis=0) - 1) < 1e-2))
+        # perform Shapiro-Wilk test for normality
+        for i in (0, 1):
+            _, pval = shapiro(samples[:, i])
+            assert_(pval > 0.9)
+        # make sure samples are uncorrelated
+        cov = np.cov(samples.transpose())
+        assert_(np.abs(cov[0, 1]) < 1e-2)
+
+        # test the correlated, non-zero mean case
+        engine = qmc.MultivariateNormalQMC(
+            mean=[1.0, 2.0], cov=[[1.5, 0.5], [0.5, 1.5]], seed=12345
+        )
+        samples = engine.random(n=256)
+        assert_(all(np.abs(samples.mean(axis=0) - [1, 2]) < 1e-2))
+        assert_(all(np.abs(samples.std(axis=0) - np.sqrt(1.5)) < 1e-2))
+        # perform Shapiro-Wilk test for normality
+        for i in (0, 1):
+            _, pval = shapiro(samples[:, i])
+            assert_(pval > 0.9)
+        # check covariance
+        cov = np.cov(samples.transpose())
+        assert_(np.abs(cov[0, 1] - 0.5) < 1e-2)
+
+    def test_MultivariateNormalQMCShapiroInvTransform(self):
+        # test the standard case
+        engine = qmc.MultivariateNormalQMC(
+            mean=[0, 0], cov=[[1, 0], [0, 1]], seed=12345, inv_transform=True
+        )
+        samples = engine.random(n=256)
+        assert_(all(np.abs(samples.mean(axis=0)) < 1e-2))
+        assert_(all(np.abs(samples.std(axis=0) - 1) < 1e-2))
+        # perform Shapiro-Wilk test for normality
+        for i in (0, 1):
+            _, pval = shapiro(samples[:, i])
+            assert_(pval > 0.9)
+        # make sure samples are uncorrelated
+        cov = np.cov(samples.transpose())
+        assert_(np.abs(cov[0, 1]) < 1e-2)
+
+        # test the correlated, non-zero mean case
+        engine = qmc.MultivariateNormalQMC(
+            mean=[1.0, 2.0],
+            cov=[[1.5, 0.5], [0.5, 1.5]],
+            seed=12345,
+            inv_transform=True,
+        )
+        samples = engine.random(n=256)
+        assert_(all(np.abs(samples.mean(axis=0) - [1, 2]) < 1e-2))
+        assert_(all(np.abs(samples.std(axis=0) - np.sqrt(1.5)) < 1e-2))
+        # perform Shapiro-Wilk test for normality
+        for i in (0, 1):
+            _, pval = shapiro(samples[:, i])
+            assert_(pval > 0.9)
+        # check covariance
+        cov = np.cov(samples.transpose())
+        assert_(np.abs(cov[0, 1] - 0.5) < 1e-2)
+
+    def test_MultivariateNormalQMCDegenerate(self):
+        # X, Y iid standard Normal and Z = X + Y, random vector (X, Y, Z)
+        engine = qmc.MultivariateNormalQMC(
+            mean=[0.0, 0.0, 0.0],
+            cov=[[1.0, 0.0, 1.0], [0.0, 1.0, 1.0], [1.0, 1.0, 2.0]],
+            seed=12345,
+        )
+        samples = engine.random(n=512)
+        assert_(all(np.abs(samples.mean(axis=0)) < 1e-2))
+        assert_(np.abs(np.std(samples[:, 0]) - 1) < 1e-2)
+        assert_(np.abs(np.std(samples[:, 1]) - 1) < 1e-2)
+        assert_(np.abs(np.std(samples[:, 2]) - np.sqrt(2)) < 1e-2)
+        for i in (0, 1, 2):
+            _, pval = shapiro(samples[:, i])
+            assert_(pval > 0.8)
+        cov = np.cov(samples.transpose())
+        assert_(np.abs(cov[0, 1]) < 1e-2)
+        assert_(np.abs(cov[0, 2] - 1) < 1e-2)
+        # check to see if X + Y = Z almost exactly
+        assert_(
+            all(np.abs(samples[:, 0] + samples[:, 1] - samples[:, 2]) < 1e-5)
+        )
diff --git a/scipy/stats/tests/test_sobol_.py b/scipy/stats/tests/test_sobol_.py
new file mode 100644
index 000000000000..2f62180bcfc9
--- /dev/null
+++ b/scipy/stats/tests/test_sobol_.py
@@ -0,0 +1,212 @@
+from collections import Counter
+
+import pytest
+import numpy as np
+from numpy.testing import (assert_, assert_equal, assert_array_almost_equal,
+                           assert_array_equal)
+
+from scipy.stats.qmc import Sobol
+import scipy.stats.qmc as qmc
+
+
+class TestSobol:
+
+    @pytest.fixture(autouse=True)
+    def setup_sequences(self):
+        # 10 samples for testing. Silent the warning as 10 is not 2**n
+        with pytest.warns(UserWarning):
+            engine_unscrambled_1d = Sobol(1, scramble=False)
+            self.draws_unscrambled_1d = engine_unscrambled_1d.random(10)
+            engine_unscrambled_3d = Sobol(3, scramble=False)
+            self.draws_unscrambled_3d = engine_unscrambled_3d.random(10)
+            engine_scrambled_1d = Sobol(1, scramble=True, seed=12345)
+            self.draws_scrambled_1d = engine_scrambled_1d.random(10)
+            engine_scrambled_3d = Sobol(3, scramble=True, seed=12345)
+            self.draws_scrambled_3d = engine_scrambled_3d.random(10)
+
+    def test_warning(self):
+        with pytest.warns(UserWarning):
+            engine = Sobol(1)
+            engine.random(10)
+
+    def test_random_base2(self):
+        engine = Sobol(1, scramble=False)
+        sample = engine.random_base2(2)
+        assert_array_equal(self.draws_unscrambled_1d[:4],
+                           sample)
+
+        # resampling still having N=2**n
+        sample = engine.random_base2(2)
+        assert_array_equal(self.draws_unscrambled_1d[4:8],
+                           sample)
+
+        # resampling again but leading to N!=2**n
+        with pytest.raises(ValueError, match=r"The balance properties of "
+                                             r"Sobol' points"):
+            engine.random_base2(2)
+
+    def test_raise(self):
+        with pytest.raises(ValueError, match=r"Maximum supported "
+                                             r"dimensionality"):
+            Sobol(Sobol.MAXDIM + 1)
+
+    def test_Unscrambled1DSobol(self):
+        expected = [
+            0, 0.5, 0.75, 0.25, 0.375, 0.875, 0.625, 0.125, 0.1875, 0.6875
+        ]
+        assert_equal(self.draws_unscrambled_1d.shape[0], 10)
+        assert_equal(self.draws_unscrambled_1d.shape[1], 1)
+        assert_array_equal(self.draws_unscrambled_1d.flatten(),
+                           np.array(expected))
+
+    def test_Unscrambled3DSobol(self):
+        expected_dim3 = [0, 0.5, 0.25, 0.75, 0.625, 0.125, 0.875, 0.375,
+                         0.9375, 0.4375]
+        assert_equal(self.draws_unscrambled_3d.shape[0], 10)
+        assert_equal(self.draws_unscrambled_3d.shape[1], 3)
+        assert_array_equal(self.draws_unscrambled_3d[:, 2],
+                           np.array(expected_dim3))
+        assert_array_equal(
+            self.draws_unscrambled_3d[:, 0],
+            self.draws_unscrambled_1d.flatten()
+        )
+
+    def test_Unscrambled3DAsyncSobol(self):
+        engine_unscrambled_3d = Sobol(3, scramble=False)
+        draws = np.vstack([engine_unscrambled_3d.random() for _ in range(10)])
+        assert_array_equal(self.draws_unscrambled_3d, draws)
+
+    def test_UnscrambledFastForwardAndResetSobol(self):
+        engine_unscrambled_3d = Sobol(3, scramble=False).fast_forward(5)
+        draws = engine_unscrambled_3d.random(5)
+        assert_array_equal(self.draws_unscrambled_3d[5:10, :], draws)
+
+        engine_unscrambled_3d.reset()
+        even_draws = []
+        for i in range(10):
+            if i % 2 == 0:
+                even_draws.append(engine_unscrambled_3d.random())
+            else:
+                engine_unscrambled_3d.fast_forward(1)
+        assert_array_equal(
+            self.draws_unscrambled_3d[[i for i in range(10) if i % 2 == 0]],
+            np.vstack(even_draws),
+        )
+
+    def test_UnscrambledHighDimSobol(self):
+        engine = Sobol(1111, scramble=False)
+        count1 = Counter(engine.random().flatten().tolist())
+        count2 = Counter(engine.random().flatten().tolist())
+        count3 = Counter(engine.random().flatten().tolist())
+        assert_equal(count1, Counter({0.0: 1111}))
+        assert_equal(count2, Counter({0.5: 1111}))
+        assert_equal(count3, Counter({0.25: 557, 0.75: 554}))
+
+    def test_UnscrambledSobolBounds(self):
+        engine = Sobol(1111, scramble=False)
+        draws = engine.random(1024)
+        assert_(np.all(draws >= 0))
+        assert_(np.all(draws <= 1))
+
+    def test_UnscrambledDistributionSobol(self):
+        engine = Sobol(1111, scramble=False)
+        draws = engine.random(1024)
+        assert_array_almost_equal(
+            np.mean(draws, axis=0), np.repeat(0.5, 1111), decimal=2
+        )
+        assert_array_almost_equal(
+            np.percentile(draws, 25, axis=0), np.repeat(0.25, 1111), decimal=2
+        )
+        assert_array_almost_equal(
+            np.percentile(draws, 75, axis=0), np.repeat(0.75, 1111), decimal=2
+        )
+
+    def test_Scrambled1DSobol(self):
+        expected = [
+            0.590297,
+            0.46784395,
+            0.03562005,
+            0.91319746,
+            0.86014303,
+            0.23796839,
+            0.25856809,
+            0.63636296,
+            0.69455189,
+            0.316758
+        ]
+        print(self.draws_scrambled_1d.flatten())
+        assert_equal(self.draws_scrambled_1d.shape[0], 10)
+        assert_equal(self.draws_scrambled_1d.shape[1], 1)
+        assert_array_almost_equal(
+            self.draws_scrambled_1d.flatten(), np.array(expected)
+        )
+
+    def test_Scrambled3DSobol(self):
+        expected_dim3 = [0.56645, 0.19712, 0.79965, 0.43654, 0.0867, 0.70811,
+                         0.295, 0.90994, 0.31903, 0.94863]
+        assert_equal(self.draws_scrambled_3d.shape[0], 10)
+        assert_equal(self.draws_scrambled_3d.shape[1], 3)
+        assert_array_almost_equal(
+            self.draws_scrambled_3d[:, 2], np.array(expected_dim3), decimal=5
+        )
+
+    def test_Scrambled3DAsyncSobol(self):
+        engine_unscrambled_3d = Sobol(3, scramble=False)
+        draws = np.vstack([engine_unscrambled_3d.random() for _ in range(10)])
+        assert_array_equal(self.draws_unscrambled_3d, draws)
+
+    def test_ScrambledSobolBounds(self):
+        engine = Sobol(100, scramble=True)
+        draws = engine.random(1024)
+        assert_(np.all(draws >= 0))
+        assert_(np.all(draws <= 1))
+
+    def test_ScrambledFastForwardAndResetSobol(self):
+        engine_scrambled_3d = Sobol(3, scramble=True,
+                                    seed=12345).fast_forward(5)
+        draws = engine_scrambled_3d.random(5)
+        assert_array_equal(self.draws_scrambled_3d[5:10], draws)
+
+        engine_scrambled_3d.reset()
+        even_draws = []
+        for i in range(10):
+            if i % 2 == 0:
+                even_draws.append(engine_scrambled_3d.random())
+            else:
+                engine_scrambled_3d.fast_forward(1)
+        assert_array_equal(
+            self.draws_scrambled_3d[[i for i in range(10) if i % 2 == 0]],
+            np.vstack(even_draws),
+        )
+
+    def test_ScrambledDistributionSobol(self):
+        engine = Sobol(10, scramble=True, seed=12345)
+        draws = engine.random(512)
+        assert_array_almost_equal(
+            np.mean(draws, axis=0), np.repeat(0.5, 10), decimal=2
+        )
+        assert_array_almost_equal(
+            np.percentile(draws, 25, axis=0), np.repeat(0.25, 10), decimal=2
+        )
+        assert_array_almost_equal(
+            np.percentile(draws, 75, axis=0), np.repeat(0.75, 10), decimal=2
+        )
+
+    def test_0Dim(self):
+        engine = Sobol(0, scramble=False)
+        draws = engine.random(4)
+        assert_array_equal(np.empty((4, 0)), draws)
+
+    def test_discrepancy(self):
+        engine_sobol = Sobol(10, scramble=False)
+        sample_sobol = engine_sobol.random(128)
+
+        engine_olhs = qmc.OrthogonalLatinHypercube(10)
+        sample_olhs = engine_olhs.random(128)
+
+        assert qmc.discrepancy(sample_sobol) < qmc.discrepancy(sample_olhs)
+
+        engine_halton = qmc.Halton(10)
+        sample_halton = engine_halton.random(128)
+
+        assert qmc.discrepancy(sample_sobol) < qmc.discrepancy(sample_halton)
diff --git a/tools/refguide_check.py b/tools/refguide_check.py
index 5b8c3910b404..8aafd377fd62 100755
--- a/tools/refguide_check.py
+++ b/tools/refguide_check.py
@@ -97,6 +97,7 @@
     'stats',
     'stats.mstats',
     'stats.contingency',
+    'stats.qmc',
 ]
 
 # Docs for these modules are included in the parent module


Code review comments:
diff_hunk: @@ -0,0 +1,212 @@
+from collections import Counter
+
+import pytest
+import numpy as np
+from numpy.testing import (assert_, assert_equal, assert_array_almost_equal,
+                           assert_array_equal)
+
+from scipy.stats.qmc import Sobol
+import scipy.stats.qmc as qmc
+
+
+class TestSobol:
+
+    @pytest.fixture(autouse=True)
+    def setup_sequences(self):
+        # 10 samples for testing. Silent the warning as 10 is not 2**n
+        with pytest.warns(UserWarning):
+            engine_unscrambled_1d = Sobol(1, scramble=False)
+            self.draws_unscrambled_1d = engine_unscrambled_1d.random(10)
+            engine_unscrambled_3d = Sobol(3, scramble=False)
+            self.draws_unscrambled_3d = engine_unscrambled_3d.random(10)
+            engine_scrambled_1d = Sobol(1, scramble=True, seed=12345)
+            self.draws_scrambled_1d = engine_scrambled_1d.random(10)
+            engine_scrambled_3d = Sobol(3, scramble=True, seed=12345)
+            self.draws_scrambled_3d = engine_scrambled_3d.random(10)
+
+    def test_warning(self):
+        with pytest.warns(UserWarning):
+            engine = Sobol(1)
+            engine.random(10)
+
+    def test_random_base2(self):
+        engine = Sobol(1, scramble=False)
+        sample = engine.random_base2(2)
+        assert_array_equal(self.draws_unscrambled_1d[:4],
+                           sample)
+
+        # resampling still having N=2**n
+        sample = engine.random_base2(2)
+        assert_array_equal(self.draws_unscrambled_1d[4:8],
+                           sample)
+
+        # resampling again but leading to N!=2**n
+        with pytest.raises(ValueError, match=r"The balance properties of "
+                                             r"Sobol' points"):
+            engine.random_base2(2)
+
+    def test_raise(self):
+        with pytest.raises(ValueError, match=r"Maximum supported "
+                                             r"dimensionality"):
+            Sobol(Sobol.MAXDIM + 1)
+
+    def test_Unscrambled1DSobol(self):
+        expected = [
+            0, 0.5, 0.75, 0.25, 0.375, 0.875, 0.625, 0.125, 0.1875, 0.6875
+        ]
+        assert_equal(self.draws_unscrambled_1d.shape[0], 10)
+        assert_equal(self.draws_unscrambled_1d.shape[1], 1)
+        assert_array_equal(self.draws_unscrambled_1d.flatten(),
+                           np.array(expected))
+
+    def test_Unscrambled3DSobol(self):
+        expected_dim3 = [0, 0.5, 0.25, 0.75, 0.625, 0.125, 0.875, 0.375,
+                         0.9375, 0.4375]
+        assert_equal(self.draws_unscrambled_3d.shape[0], 10)
+        assert_equal(self.draws_unscrambled_3d.shape[1], 3)
+        assert_array_equal(self.draws_unscrambled_3d[:, 2],
+                           np.array(expected_dim3))
+        assert_array_equal(
+            self.draws_unscrambled_3d[:, 0],
+            self.draws_unscrambled_1d.flatten()
+        )
+
+    def test_Unscrambled3DAsyncSobol(self):
+        engine_unscrambled_3d = Sobol(3, scramble=False)
+        draws = np.vstack([engine_unscrambled_3d.random() for _ in range(10)])
+        assert_array_equal(self.draws_unscrambled_3d, draws)
+
+    def test_UnscrambledFastForwardAndResetSobol(self):
+        engine_unscrambled_3d = Sobol(3, scramble=False).fast_forward(5)
+        draws = engine_unscrambled_3d.random(5)
+        assert_array_equal(self.draws_unscrambled_3d[5:10, :], draws)
+
+        engine_unscrambled_3d.reset()
+        even_draws = []
+        for i in range(10):
+            if i % 2 == 0:
+                even_draws.append(engine_unscrambled_3d.random())
+            else:
+                engine_unscrambled_3d.fast_forward(1)
+        assert_array_equal(
+            self.draws_unscrambled_3d[[i for i in range(10) if i % 2 == 0]],
+            np.vstack(even_draws),
+        )
+
+    def test_UnscrambledHighDimSobol(self):
+        engine = Sobol(1111, scramble=False)
+        count1 = Counter(engine.random().flatten().tolist())
+        count2 = Counter(engine.random().flatten().tolist())
+        count3 = Counter(engine.random().flatten().tolist())
+        assert_equal(count1, Counter({0.0: 1111}))
+        assert_equal(count2, Counter({0.5: 1111}))
+        assert_equal(count3, Counter({0.25: 557, 0.75: 554}))
+
+    def test_UnscrambledSobolBounds(self):
+        engine = Sobol(1111, scramble=False)
+        draws = engine.random(1024)
+        assert_(np.all(draws >= 0))
+        assert_(np.all(draws <= 1))
+
+    def test_UnscrambledDistributionSobol(self):
+        engine = Sobol(1111, scramble=False)
+        draws = engine.random(1024)
+        assert_array_almost_equal(
+            np.mean(draws, axis=0), np.repeat(0.5, 1111), decimal=2
+        )
+        assert_array_almost_equal(
+            np.percentile(draws, 25, axis=0), np.repeat(0.25, 1111), decimal=2
+        )
+        assert_array_almost_equal(
+            np.percentile(draws, 75, axis=0), np.repeat(0.75, 1111), decimal=2
+        )
+
+    def test_Scrambled1DSobol(self):
+        expected = [
+            0.590297,
+            0.46784395,
+            0.03562005,
+            0.91319746,
+            0.86014303,
+            0.23796839,
+            0.25856809,
+            0.63636296,
+            0.69455189,
+            0.316758
+        ]
+        print(self.draws_scrambled_1d.flatten())
+        assert_equal(self.draws_scrambled_1d.shape[0], 10)
+        assert_equal(self.draws_scrambled_1d.shape[1], 1)
+        assert_array_almost_equal(
+            self.draws_scrambled_1d.flatten(), np.array(expected)
+        )
+
+    def test_Scrambled3DSobol(self):
+        expected_dim3 = [0.56645, 0.19712, 0.79965, 0.43654, 0.0867, 0.70811,
+                         0.295, 0.90994, 0.31903, 0.94863]
+        assert_equal(self.draws_scrambled_3d.shape[0], 10)
+        assert_equal(self.draws_scrambled_3d.shape[1], 3)
+        assert_array_almost_equal(
+            self.draws_scrambled_3d[:, 2], np.array(expected_dim3), decimal=5
+        )
+
+    def test_Scrambled3DAsyncSobol(self):
+        engine_unscrambled_3d = Sobol(3, scramble=False)
+        draws = np.vstack([engine_unscrambled_3d.random() for _ in range(10)])
+        assert_array_equal(self.draws_unscrambled_3d, draws)
+
+    def test_ScrambledSobolBounds(self):
path: scipy/stats/tests/test_sobol_.py
text: @user1:
Code reuse with `test_ScrambledSobolBounds` above?

diff_hunk: @@ -0,0 +1,212 @@
+from collections import Counter
+
+import pytest
+import numpy as np
+from numpy.testing import (assert_, assert_equal, assert_array_almost_equal,
+                           assert_array_equal)
+
+from scipy.stats.qmc import Sobol
+import scipy.stats.qmc as qmc
+
+
+class TestSobol:
+
+    @pytest.fixture(autouse=True)
+    def setup_sequences(self):
+        # 10 samples for testing. Silent the warning as 10 is not 2**n
+        with pytest.warns(UserWarning):
+            engine_unscrambled_1d = Sobol(1, scramble=False)
+            self.draws_unscrambled_1d = engine_unscrambled_1d.random(10)
+            engine_unscrambled_3d = Sobol(3, scramble=False)
+            self.draws_unscrambled_3d = engine_unscrambled_3d.random(10)
+            engine_scrambled_1d = Sobol(1, scramble=True, seed=12345)
+            self.draws_scrambled_1d = engine_scrambled_1d.random(10)
+            engine_scrambled_3d = Sobol(3, scramble=True, seed=12345)
+            self.draws_scrambled_3d = engine_scrambled_3d.random(10)
+
+    def test_warning(self):
+        with pytest.warns(UserWarning):
+            engine = Sobol(1)
+            engine.random(10)
+
+    def test_random_base2(self):
+        engine = Sobol(1, scramble=False)
+        sample = engine.random_base2(2)
+        assert_array_equal(self.draws_unscrambled_1d[:4],
+                           sample)
+
+        # resampling still having N=2**n
+        sample = engine.random_base2(2)
+        assert_array_equal(self.draws_unscrambled_1d[4:8],
+                           sample)
+
+        # resampling again but leading to N!=2**n
+        with pytest.raises(ValueError, match=r"The balance properties of "
+                                             r"Sobol' points"):
+            engine.random_base2(2)
+
+    def test_raise(self):
+        with pytest.raises(ValueError, match=r"Maximum supported "
+                                             r"dimensionality"):
+            Sobol(Sobol.MAXDIM + 1)
+
+    def test_Unscrambled1DSobol(self):
+        expected = [
+            0, 0.5, 0.75, 0.25, 0.375, 0.875, 0.625, 0.125, 0.1875, 0.6875
+        ]
+        assert_equal(self.draws_unscrambled_1d.shape[0], 10)
+        assert_equal(self.draws_unscrambled_1d.shape[1], 1)
+        assert_array_equal(self.draws_unscrambled_1d.flatten(),
+                           np.array(expected))
+
+    def test_Unscrambled3DSobol(self):
+        expected_dim3 = [0, 0.5, 0.25, 0.75, 0.625, 0.125, 0.875, 0.375,
+                         0.9375, 0.4375]
+        assert_equal(self.draws_unscrambled_3d.shape[0], 10)
+        assert_equal(self.draws_unscrambled_3d.shape[1], 3)
+        assert_array_equal(self.draws_unscrambled_3d[:, 2],
+                           np.array(expected_dim3))
+        assert_array_equal(
+            self.draws_unscrambled_3d[:, 0],
+            self.draws_unscrambled_1d.flatten()
+        )
+
+    def test_Unscrambled3DAsyncSobol(self):
+        engine_unscrambled_3d = Sobol(3, scramble=False)
+        draws = np.vstack([engine_unscrambled_3d.random() for _ in range(10)])
+        assert_array_equal(self.draws_unscrambled_3d, draws)
+
+    def test_UnscrambledFastForwardAndResetSobol(self):
+        engine_unscrambled_3d = Sobol(3, scramble=False).fast_forward(5)
+        draws = engine_unscrambled_3d.random(5)
+        assert_array_equal(self.draws_unscrambled_3d[5:10, :], draws)
+
+        engine_unscrambled_3d.reset()
+        even_draws = []
+        for i in range(10):
+            if i % 2 == 0:
+                even_draws.append(engine_unscrambled_3d.random())
+            else:
+                engine_unscrambled_3d.fast_forward(1)
+        assert_array_equal(
+            self.draws_unscrambled_3d[[i for i in range(10) if i % 2 == 0]],
+            np.vstack(even_draws),
+        )
+
+    def test_UnscrambledHighDimSobol(self):
+        engine = Sobol(1111, scramble=False)
+        count1 = Counter(engine.random().flatten().tolist())
+        count2 = Counter(engine.random().flatten().tolist())
+        count3 = Counter(engine.random().flatten().tolist())
+        assert_equal(count1, Counter({0.0: 1111}))
+        assert_equal(count2, Counter({0.5: 1111}))
+        assert_equal(count3, Counter({0.25: 557, 0.75: 554}))
+
+    def test_UnscrambledSobolBounds(self):
+        engine = Sobol(1111, scramble=False)
+        draws = engine.random(1024)
+        assert_(np.all(draws >= 0))
+        assert_(np.all(draws <= 1))
+
+    def test_UnscrambledDistributionSobol(self):
+        engine = Sobol(1111, scramble=False)
+        draws = engine.random(1024)
+        assert_array_almost_equal(
+            np.mean(draws, axis=0), np.repeat(0.5, 1111), decimal=2
+        )
+        assert_array_almost_equal(
+            np.percentile(draws, 25, axis=0), np.repeat(0.25, 1111), decimal=2
+        )
+        assert_array_almost_equal(
+            np.percentile(draws, 75, axis=0), np.repeat(0.75, 1111), decimal=2
+        )
+
+    def test_Scrambled1DSobol(self):
+        expected = [
+            0.590297,
+            0.46784395,
+            0.03562005,
+            0.91319746,
+            0.86014303,
+            0.23796839,
+            0.25856809,
+            0.63636296,
+            0.69455189,
+            0.316758
+        ]
+        print(self.draws_scrambled_1d.flatten())
+        assert_equal(self.draws_scrambled_1d.shape[0], 10)
+        assert_equal(self.draws_scrambled_1d.shape[1], 1)
+        assert_array_almost_equal(
+            self.draws_scrambled_1d.flatten(), np.array(expected)
+        )
+
+    def test_Scrambled3DSobol(self):
+        expected_dim3 = [0.56645, 0.19712, 0.79965, 0.43654, 0.0867, 0.70811,
+                         0.295, 0.90994, 0.31903, 0.94863]
+        assert_equal(self.draws_scrambled_3d.shape[0], 10)
+        assert_equal(self.draws_scrambled_3d.shape[1], 3)
+        assert_array_almost_equal(
+            self.draws_scrambled_3d[:, 2], np.array(expected_dim3), decimal=5
+        )
+
+    def test_Scrambled3DAsyncSobol(self):
+        engine_unscrambled_3d = Sobol(3, scramble=False)
+        draws = np.vstack([engine_unscrambled_3d.random() for _ in range(10)])
+        assert_array_equal(self.draws_unscrambled_3d, draws)
+
+    def test_ScrambledSobolBounds(self):
+        engine = Sobol(100, scramble=True)
+        draws = engine.random(1024)
+        assert_(np.all(draws >= 0))
+        assert_(np.all(draws <= 1))
+
+    def test_ScrambledFastForwardAndResetSobol(self):
+        engine_scrambled_3d = Sobol(3, scramble=True,
+                                    seed=12345).fast_forward(5)
+        draws = engine_scrambled_3d.random(5)
+        assert_array_equal(self.draws_scrambled_3d[5:10], draws)
+
+        engine_scrambled_3d.reset()
+        even_draws = []
+        for i in range(10):
+            if i % 2 == 0:
+                even_draws.append(engine_scrambled_3d.random())
+            else:
+                engine_scrambled_3d.fast_forward(1)
+        assert_array_equal(
+            self.draws_scrambled_3d[[i for i in range(10) if i % 2 == 0]],
+            np.vstack(even_draws),
+        )
+
+    def test_ScrambledDistributionSobol(self):
path: scipy/stats/tests/test_sobol_.py
text: @user1:
Code reuse with `test_UnscrambledDistributionSobol` above?

diff_hunk: @@ -0,0 +1,212 @@
+from collections import Counter
+
+import pytest
+import numpy as np
+from numpy.testing import (assert_, assert_equal, assert_array_almost_equal,
+                           assert_array_equal)
+
+from scipy.stats.qmc import Sobol
+import scipy.stats.qmc as qmc
+
+
+class TestSobol:
+
+    @pytest.fixture(autouse=True)
+    def setup_sequences(self):
+        # 10 samples for testing. Silent the warning as 10 is not 2**n
+        with pytest.warns(UserWarning):
+            engine_unscrambled_1d = Sobol(1, scramble=False)
+            self.draws_unscrambled_1d = engine_unscrambled_1d.random(10)
+            engine_unscrambled_3d = Sobol(3, scramble=False)
+            self.draws_unscrambled_3d = engine_unscrambled_3d.random(10)
+            engine_scrambled_1d = Sobol(1, scramble=True, seed=12345)
+            self.draws_scrambled_1d = engine_scrambled_1d.random(10)
+            engine_scrambled_3d = Sobol(3, scramble=True, seed=12345)
+            self.draws_scrambled_3d = engine_scrambled_3d.random(10)
+
+    def test_warning(self):
+        with pytest.warns(UserWarning):
+            engine = Sobol(1)
+            engine.random(10)
+
+    def test_random_base2(self):
+        engine = Sobol(1, scramble=False)
+        sample = engine.random_base2(2)
+        assert_array_equal(self.draws_unscrambled_1d[:4],
+                           sample)
+
+        # resampling still having N=2**n
+        sample = engine.random_base2(2)
+        assert_array_equal(self.draws_unscrambled_1d[4:8],
+                           sample)
+
+        # resampling again but leading to N!=2**n
+        with pytest.raises(ValueError, match=r"The balance properties of "
+                                             r"Sobol' points"):
+            engine.random_base2(2)
+
+    def test_raise(self):
+        with pytest.raises(ValueError, match=r"Maximum supported "
+                                             r"dimensionality"):
+            Sobol(Sobol.MAXDIM + 1)
+
+    def test_Unscrambled1DSobol(self):
+        expected = [
+            0, 0.5, 0.75, 0.25, 0.375, 0.875, 0.625, 0.125, 0.1875, 0.6875
+        ]
+        assert_equal(self.draws_unscrambled_1d.shape[0], 10)
+        assert_equal(self.draws_unscrambled_1d.shape[1], 1)
+        assert_array_equal(self.draws_unscrambled_1d.flatten(),
+                           np.array(expected))
+
+    def test_Unscrambled3DSobol(self):
+        expected_dim3 = [0, 0.5, 0.25, 0.75, 0.625, 0.125, 0.875, 0.375,
+                         0.9375, 0.4375]
+        assert_equal(self.draws_unscrambled_3d.shape[0], 10)
+        assert_equal(self.draws_unscrambled_3d.shape[1], 3)
+        assert_array_equal(self.draws_unscrambled_3d[:, 2],
+                           np.array(expected_dim3))
+        assert_array_equal(
+            self.draws_unscrambled_3d[:, 0],
+            self.draws_unscrambled_1d.flatten()
+        )
+
+    def test_Unscrambled3DAsyncSobol(self):
+        engine_unscrambled_3d = Sobol(3, scramble=False)
+        draws = np.vstack([engine_unscrambled_3d.random() for _ in range(10)])
+        assert_array_equal(self.draws_unscrambled_3d, draws)
+
+    def test_UnscrambledFastForwardAndResetSobol(self):
+        engine_unscrambled_3d = Sobol(3, scramble=False).fast_forward(5)
+        draws = engine_unscrambled_3d.random(5)
+        assert_array_equal(self.draws_unscrambled_3d[5:10, :], draws)
+
+        engine_unscrambled_3d.reset()
+        even_draws = []
+        for i in range(10):
+            if i % 2 == 0:
+                even_draws.append(engine_unscrambled_3d.random())
+            else:
+                engine_unscrambled_3d.fast_forward(1)
+        assert_array_equal(
+            self.draws_unscrambled_3d[[i for i in range(10) if i % 2 == 0]],
+            np.vstack(even_draws),
+        )
+
+    def test_UnscrambledHighDimSobol(self):
+        engine = Sobol(1111, scramble=False)
+        count1 = Counter(engine.random().flatten().tolist())
+        count2 = Counter(engine.random().flatten().tolist())
+        count3 = Counter(engine.random().flatten().tolist())
+        assert_equal(count1, Counter({0.0: 1111}))
+        assert_equal(count2, Counter({0.5: 1111}))
+        assert_equal(count3, Counter({0.25: 557, 0.75: 554}))
+
+    def test_UnscrambledSobolBounds(self):
+        engine = Sobol(1111, scramble=False)
+        draws = engine.random(1024)
+        assert_(np.all(draws >= 0))
+        assert_(np.all(draws <= 1))
+
+    def test_UnscrambledDistributionSobol(self):
+        engine = Sobol(1111, scramble=False)
+        draws = engine.random(1024)
+        assert_array_almost_equal(
+            np.mean(draws, axis=0), np.repeat(0.5, 1111), decimal=2
+        )
+        assert_array_almost_equal(
+            np.percentile(draws, 25, axis=0), np.repeat(0.25, 1111), decimal=2
+        )
+        assert_array_almost_equal(
+            np.percentile(draws, 75, axis=0), np.repeat(0.75, 1111), decimal=2
+        )
+
+    def test_Scrambled1DSobol(self):
+        expected = [
+            0.590297,
+            0.46784395,
+            0.03562005,
+            0.91319746,
+            0.86014303,
+            0.23796839,
+            0.25856809,
+            0.63636296,
+            0.69455189,
+            0.316758
+        ]
+        print(self.draws_scrambled_1d.flatten())
+        assert_equal(self.draws_scrambled_1d.shape[0], 10)
+        assert_equal(self.draws_scrambled_1d.shape[1], 1)
+        assert_array_almost_equal(
+            self.draws_scrambled_1d.flatten(), np.array(expected)
+        )
+
+    def test_Scrambled3DSobol(self):
+        expected_dim3 = [0.56645, 0.19712, 0.79965, 0.43654, 0.0867, 0.70811,
+                         0.295, 0.90994, 0.31903, 0.94863]
+        assert_equal(self.draws_scrambled_3d.shape[0], 10)
+        assert_equal(self.draws_scrambled_3d.shape[1], 3)
+        assert_array_almost_equal(
+            self.draws_scrambled_3d[:, 2], np.array(expected_dim3), decimal=5
+        )
+
+    def test_Scrambled3DAsyncSobol(self):
+        engine_unscrambled_3d = Sobol(3, scramble=False)
+        draws = np.vstack([engine_unscrambled_3d.random() for _ in range(10)])
+        assert_array_equal(self.draws_unscrambled_3d, draws)
+
+    def test_ScrambledSobolBounds(self):
+        engine = Sobol(100, scramble=True)
+        draws = engine.random(1024)
+        assert_(np.all(draws >= 0))
+        assert_(np.all(draws <= 1))
+
+    def test_ScrambledFastForwardAndResetSobol(self):
path: scipy/stats/tests/test_sobol_.py
text: @user1:
Code reuse with `test_UnscrambledFastForwardAndResetSobol` above?

diff_hunk: @@ -0,0 +1,212 @@
+from collections import Counter
+
+import pytest
+import numpy as np
+from numpy.testing import (assert_, assert_equal, assert_array_almost_equal,
+                           assert_array_equal)
+
+from scipy.stats.qmc import Sobol
+import scipy.stats.qmc as qmc
+
+
+class TestSobol:
+
+    @pytest.fixture(autouse=True)
+    def setup_sequences(self):
+        # 10 samples for testing. Silent the warning as 10 is not 2**n
+        with pytest.warns(UserWarning):
+            engine_unscrambled_1d = Sobol(1, scramble=False)
+            self.draws_unscrambled_1d = engine_unscrambled_1d.random(10)
+            engine_unscrambled_3d = Sobol(3, scramble=False)
+            self.draws_unscrambled_3d = engine_unscrambled_3d.random(10)
+            engine_scrambled_1d = Sobol(1, scramble=True, seed=12345)
+            self.draws_scrambled_1d = engine_scrambled_1d.random(10)
+            engine_scrambled_3d = Sobol(3, scramble=True, seed=12345)
+            self.draws_scrambled_3d = engine_scrambled_3d.random(10)
+
+    def test_warning(self):
+        with pytest.warns(UserWarning):
+            engine = Sobol(1)
+            engine.random(10)
+
+    def test_random_base2(self):
+        engine = Sobol(1, scramble=False)
+        sample = engine.random_base2(2)
+        assert_array_equal(self.draws_unscrambled_1d[:4],
+                           sample)
+
+        # resampling still having N=2**n
+        sample = engine.random_base2(2)
+        assert_array_equal(self.draws_unscrambled_1d[4:8],
+                           sample)
+
+        # resampling again but leading to N!=2**n
+        with pytest.raises(ValueError, match=r"The balance properties of "
+                                             r"Sobol' points"):
+            engine.random_base2(2)
+
+    def test_raise(self):
+        with pytest.raises(ValueError, match=r"Maximum supported "
+                                             r"dimensionality"):
+            Sobol(Sobol.MAXDIM + 1)
+
+    def test_Unscrambled1DSobol(self):
+        expected = [
+            0, 0.5, 0.75, 0.25, 0.375, 0.875, 0.625, 0.125, 0.1875, 0.6875
+        ]
+        assert_equal(self.draws_unscrambled_1d.shape[0], 10)
+        assert_equal(self.draws_unscrambled_1d.shape[1], 1)
+        assert_array_equal(self.draws_unscrambled_1d.flatten(),
+                           np.array(expected))
+
+    def test_Unscrambled3DSobol(self):
+        expected_dim3 = [0, 0.5, 0.25, 0.75, 0.625, 0.125, 0.875, 0.375,
+                         0.9375, 0.4375]
+        assert_equal(self.draws_unscrambled_3d.shape[0], 10)
+        assert_equal(self.draws_unscrambled_3d.shape[1], 3)
+        assert_array_equal(self.draws_unscrambled_3d[:, 2],
+                           np.array(expected_dim3))
+        assert_array_equal(
+            self.draws_unscrambled_3d[:, 0],
+            self.draws_unscrambled_1d.flatten()
+        )
+
+    def test_Unscrambled3DAsyncSobol(self):
+        engine_unscrambled_3d = Sobol(3, scramble=False)
+        draws = np.vstack([engine_unscrambled_3d.random() for _ in range(10)])
+        assert_array_equal(self.draws_unscrambled_3d, draws)
+
+    def test_UnscrambledFastForwardAndResetSobol(self):
+        engine_unscrambled_3d = Sobol(3, scramble=False).fast_forward(5)
+        draws = engine_unscrambled_3d.random(5)
+        assert_array_equal(self.draws_unscrambled_3d[5:10, :], draws)
+
+        engine_unscrambled_3d.reset()
+        even_draws = []
+        for i in range(10):
+            if i % 2 == 0:
+                even_draws.append(engine_unscrambled_3d.random())
+            else:
+                engine_unscrambled_3d.fast_forward(1)
+        assert_array_equal(
+            self.draws_unscrambled_3d[[i for i in range(10) if i % 2 == 0]],
+            np.vstack(even_draws),
+        )
+
+    def test_UnscrambledHighDimSobol(self):
+        engine = Sobol(1111, scramble=False)
+        count1 = Counter(engine.random().flatten().tolist())
+        count2 = Counter(engine.random().flatten().tolist())
+        count3 = Counter(engine.random().flatten().tolist())
+        assert_equal(count1, Counter({0.0: 1111}))
+        assert_equal(count2, Counter({0.5: 1111}))
+        assert_equal(count3, Counter({0.25: 557, 0.75: 554}))
+
+    def test_UnscrambledSobolBounds(self):
+        engine = Sobol(1111, scramble=False)
+        draws = engine.random(1024)
+        assert_(np.all(draws >= 0))
+        assert_(np.all(draws <= 1))
+
+    def test_UnscrambledDistributionSobol(self):
+        engine = Sobol(1111, scramble=False)
+        draws = engine.random(1024)
+        assert_array_almost_equal(
+            np.mean(draws, axis=0), np.repeat(0.5, 1111), decimal=2
+        )
+        assert_array_almost_equal(
+            np.percentile(draws, 25, axis=0), np.repeat(0.25, 1111), decimal=2
+        )
+        assert_array_almost_equal(
+            np.percentile(draws, 75, axis=0), np.repeat(0.75, 1111), decimal=2
+        )
+
+    def test_Scrambled1DSobol(self):
+        expected = [
+            0.590297,
+            0.46784395,
+            0.03562005,
+            0.91319746,
+            0.86014303,
+            0.23796839,
+            0.25856809,
+            0.63636296,
+            0.69455189,
+            0.316758
+        ]
+        print(self.draws_scrambled_1d.flatten())
+        assert_equal(self.draws_scrambled_1d.shape[0], 10)
+        assert_equal(self.draws_scrambled_1d.shape[1], 1)
+        assert_array_almost_equal(
+            self.draws_scrambled_1d.flatten(), np.array(expected)
+        )
+
+    def test_Scrambled3DSobol(self):
+        expected_dim3 = [0.56645, 0.19712, 0.79965, 0.43654, 0.0867, 0.70811,
+                         0.295, 0.90994, 0.31903, 0.94863]
+        assert_equal(self.draws_scrambled_3d.shape[0], 10)
+        assert_equal(self.draws_scrambled_3d.shape[1], 3)
+        assert_array_almost_equal(
+            self.draws_scrambled_3d[:, 2], np.array(expected_dim3), decimal=5
+        )
+
+    def test_Scrambled3DAsyncSobol(self):
path: scipy/stats/tests/test_sobol_.py
text: @user1:
This was probably copy-pasted from `test_Unscrambled3DAsyncSobol`, and it actually wasn't modified. (I assume it was supposed to be.) That suggests there should be some code-re-use - especially when you consider that this sort of test really should apply to all `QMCEngine` subclasseses.  With such a large PR, avoiding duplication will really help. 

That, combined with comments about what each tests is looking for, will help reviewers wrap their heads around the test suite as a whole to decide whether it is reasonably complete.

@author:
I fixed this so the current version is complete. In the meantime I am still thinking about doing a test class as you suggested.