FIX: global variables in cython for Sobol and tests

scipy/stats/_sobol.pyx

@@ -23,7 +23,7 @@ cdef bint is_initialized = False
 
 # Load direction numbers (taken from https://web.maths.unsw.edu.au/~fkuo/sobol/)
 def initialize_direction_numbers():
-    global is_initialized
+    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"]

scipy/stats/tests/test_qmc.py

@@ -24,9 +24,6 @@ def test_scale(self):
         scaled_back_space = qmc.scale(scaled_space, bounds=corners, reverse=True)
         assert_allclose(scaled_back_space, space)
 
-    def test_discrepancy_star_L2(self):
-
-        qmc.discrepancy_star_L2
 
     def test_discrepancy(self):
         space_0 = [[0.1, 0.5], [0.2, 0.4], [0.3, 0.3], [0.4, 0.2], [0.5, 0.1]]
@@ -251,6 +248,7 @@ def test_other_engine(self):
         assert_array_equal(qmc.multinomial_qmc(100, p, engine=engine, seed=12345),
                            expected)
 
+
 class TestNormalQMC:
     def test_NormalQMC(self):
         # d = 1
@@ -299,7 +297,7 @@ def test_NormalQMCSeeded(self):
         samples_expected = np.array(
             [
                 [1.83169884, -1.40473647, 0.24334828],
-                [0.36596099, 1.2987395, -1.47556275],
+                [0.36596099, 1.2987395, -0.330971],
             ]
         )
         assert_array_almost_equal(samples, samples_expected)
@@ -318,7 +316,7 @@ def test_NormalQMCSeededInvTransform(self):
         samples_expected = np.array(
             [
                 [-1.40525266, 1.37652443, -0.8519666],
-                [-0.166497, -2.3153681, -0.15975676],
+                [-0.166497, -2.3153681, 0.840378],
             ]
         )
         assert_array_almost_equal(samples, samples_expected)
@@ -453,7 +451,7 @@ def test_MultivariateNormalQMCSeeded(self):
         samples_expected = np.array(
             [
                 [2.05178452, -6.35744194, 0.67944512],
-                [0.40993262, 2.60517697, -1.69415825],
+                [0.40993262, 2.60517697, -0.43475],
             ]
         )
         assert_array_almost_equal(samples, samples_expected)
@@ -483,7 +481,7 @@ def test_MultivariateNormalQMCSeededInvTransform(self):
         samples_expected = np.array(
             [
                 [-1.5740992, 5.61057598, -1.28218525],
-                [-0.18650226, -5.41662685, 0.023199],
+                [-0.18650226, -5.41662685, 1.12366],
             ]
         )
         assert_array_almost_equal(samples, samples_expected)
@@ -567,7 +565,7 @@ def test_MultivariateNormalQMCDegenerate(self):
         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.9)
+            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)

scipy/stats/tests/test_sobol_.py

@@ -7,6 +7,7 @@
 from numpy.testing import (assert_, assert_equal, assert_array_almost_equal,
                            assert_array_equal)
 from pytest import raises as assert_raises
+from stats import qmc
 
 
 class TestSobol:
@@ -37,18 +38,8 @@ def test_Unscrambled1DSobol(self):
 
     def test_Unscrambled3DSobol(self):
         self.setUp()
-        expected_dim3 = [
-            0.5,
-            0.75,
-            0.25,
-            0.625,
-            0.125,
-            0.375,
-            0.875,
-            0.3125,
-            0.8125,
-            0.5625,
-        ]
+        expected_dim3 = [0.5, 0.25, 0.75, 0.625, 0.125, 0.875, 0.375, 0.9375,
+                         0.4375, 0.6875]
         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))
@@ -86,8 +77,8 @@ def test_UnscrambledHighDimSobol(self):
         count2 = Counter(engine.random().flatten().tolist())
         count3 = Counter(engine.random().flatten().tolist())
         assert_equal(count1, Counter({0.5: 1111}))
-        assert_equal(count2, Counter({0.25: 580, 0.75: 531}))
-        assert_equal(count3, Counter({0.25: 531, 0.75: 580}))
+        assert_equal(count2, Counter({0.25: 557, 0.75: 554}))
+        assert_equal(count3, Counter({0.25: 554, 0.75: 557}))
 
     def test_UnscrambledSobolBounds(self):
         engine = Sobol(1111)
@@ -131,18 +122,8 @@ def test_Scrambled1DSobol(self):
 
     def test_Scrambled3DSobol(self):
         self.setUp()
-        expected_dim3 = [
-            0.19711632,
-            0.43653634,
-            0.79965184,
-            0.08670237,
-            0.70811484,
-            0.90994149,
-            0.29499525,
-            0.83833538,
-            0.46057166,
-            0.15769824,
-        ]
+        expected_dim3 = [0.19712, 0.79965, 0.43654, 0.0867, 0.70811, 0.295,
+                         0.90994, 0.31903, 0.94863, 0.04729]
         assert_equal(self.draws_scrambled_3d.shape[0], 10)
         assert_equal(self.draws_scrambled_3d.shape[1], 3)
         assert_array_almost_equal(
@@ -196,3 +177,17 @@ def test_0Dim(self):
         engine = Sobol(0)
         draws = engine.random(5)
         assert_array_equal(np.empty((5, 0)), draws)
+
+    def test_discrepancy(self):
+        engine_sobol = Sobol(10)
+        sample_sobol = engine_sobol.random(100)
+
+        engine_olhs = qmc.OrthogonalLatinHypercube(10)
+        sample_olhs = engine_olhs.random(100)
+
+        assert qmc.discrepancy(sample_sobol) < qmc.discrepancy(sample_olhs)
+
+        engine_halton = qmc.Halton(10)
+        sample_halton = engine_halton.random(100)
+
+        assert qmc.discrepancy(sample_sobol) < qmc.discrepancy(sample_halton)