Monte Carlo Simulations In Python Using Quasi Random Standard Normal Numbers Using Sobol Sequences Gives Erroneous Values

I'm trying to perform Monte Carlo Simulations using quasi-random standard normal numbers. I understand that we can use Sobol sequences to generate uniform numbers, and then use pro

Solution 1:

It appears there is a bug in sobol_seq. Anaconda, python 3.7, 64bit, Windows 10 x64, installed sobol_seq via pip

pip install sobol_seq

# Name                    Version                   Build  Channel
sobol-seq                 0.1.2                    pypi_0    pypi

Simple code

print(sobol_seq.i4_sobol_generate(1, 1, 0))
print(sobol_seq.i4_sobol_generate(1, 1, 1))
print(sobol_seq.i4_sobol_generate(1, 1, 2))
print(sobol_seq.i4_sobol_generate(1, 1, 3))

produced output

[[0.5]][[0.5]][[0.5]][[0.5]]

Code from http://people.sc.fsu.edu/~jburkardt/py_src/sobol/sobol.html, sobol_lib.py behaves reasonable (well, except first point).

Well, enclosed code looks like it might work, keeping seed together with sampled array. Slow, though...

import sobol_seq
import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import norm

defi4_sobol_generate_std_normal(dim_num, seed, size=None):
    """
    Generates multivariate standard normal quasi-random variables.
    Parameters:
      Input, integer dim_num, the spatial dimension.
      Input, integer n, the number of points to generate.
      Input, integer seed, initial seed
      Output, real np array of shape (n, dim_num).
    """if size isNone:
        q, seed = sobol_seq.i4_sobol(dim_num, seed)
        normals = norm.ppf(q)
        return (normals, seed)

    ifisinstance(size, int) orisinstance(size, np.int32) orisinstance(size, np.int64) orisinstance(size, np.int16):
        rc = np.empty((dim_num, size))

        for k inrange(size):
            q, seed = sobol_seq.i4_sobol(dim_num, seed)
            rc[:,k] = norm.ppf(q)

        return (rc, seed)
    else:
        raise ValueError("Size type is not recognized")

    returnNone


seed = 1
x, seed = i4_sobol_generate_std_normal(1, seed)
print(x)

x, seed = i4_sobol_generate_std_normal(1, seed)
print(x)

seed = 1
x, seed = i4_sobol_generate_std_normal(1, seed, size=10)
print(x)

x, seed = i4_sobol_generate_std_normal(1, seed, size=1000)
print(x)

hist, bins = np.histogram(x, bins=20, range=(-2.5, 2.5), density=True)
plt.bar(bins[:-1], hist, width = 0.22, align='edge')

plt.show()

Here is the picture

Solution 2:

For future reference, we added Sobol' sequence in SciPy 1.7: scipy.stats.qmc.Sobol