Search code examples
pythonrandomscipy

Recommendations for Low Discrepancy (e.g. Sobol) quasi-random sequences in Python/SciPy?

I would like to use a quasi-random sequence, specifically Sobol, within a SciPy based simulation. Any recommendations on existing, efficient packages?

Solution

  • I would use OpenTURNS, which provides several low discrepancy sequences:

    • Faure sequence,
    • Halton sequence,
    • Reverse Halton sequence,
    • Haselgrove sequence,
    • Sobol sequence.

    Moreover, the sequence can be generated so that the marginals have arbitrary distribution. This is done with an probabilistic transformation, based on the inverse distribution function.

    In the following example, I generate a Sobol' sequence in 2 dimensions, based on the LowDiscrepancyExperiment class. The marginals are uniform in the [-1, 1] interval (which is the default Uniform distribution in OT). I suggest to use a sample size equal to a power of 2, because Sobol' sequence is based on base-2 integer decomposition. The generate method returns a ot.Sample.

    import openturns as ot
dim = 2
distribution = ot.ComposedDistribution([ot.Uniform()]*dim)
bounds = distribution.getRange()
sequence = ot.SobolSequence(dim)
samplesize = 2**5 # Sobol' sequences are in base 2
experiment = ot.LowDiscrepancyExperiment(sequence, distribution,
                                         samplesize, False)
sample = experiment.generate()
print(samplesize[:5])

    The previous sample has size 32. The first 5 elements are:

        y0      y1
0    0       0
1    0.5    -0.5
2   -0.5     0.5
3   -0.25   -0.25
4    0.75    0.75

    The Sobol' sequence in OT can generate an arbitrary sample size, in dimensions up to 1111.

    With a little more work, we may plot the design.

    import openturns.viewer as otv
fig = otv.PlotDesign(sample, bounds, 2**2, 2**1);
fig.set_size_inches(6, 6)

    which produces:

    A Sobol' sequence in 2 dimensions

    See how there is exactly 4 points in each elementary interval.

    If required, the sample can be easily converted into a Numpy array, which may better fits into your Scipy requirement:

    import numpy as np
array = np.array(sample)

    Other examples are provided at: http://openturns.github.io/openturns/latest/auto_reliability_sensitivity/index.html#design-of-experiments