Allow conditional independence tests to handle zero variance inputs

.circleci/config.yml

@@ -2,7 +2,6 @@ version: 2.1
 
 orbs:
   python: circleci/python@1.2
-  codecov: codecov/codecov@3.1.1
 
 jobs:
   build-and-test:
@@ -25,8 +24,11 @@ jobs:
           path: test-results
       - store_artifacts:
           path: test-results
-      - codecov/upload:
-          file: "coverage.xml"
+      - run:
+          name: Upload to Codecov
+          command: |
+            pip install codecov-cli
+            codecovcli upload-coverage --file coverage.xml
 
 workflows:
   run-tests:

codecov.yml

@@ -19,9 +19,9 @@ coverage:
 # --------------
 # which folders/files to ignore
 ignore:
-  - *tests*
-  - *__init__.py
-  - *hyppo/_utils.py
+  - "*tests*"
+  - "*__init__.py"
+  - "*hyppo/_utils.py"
 
 # Pull request comments:
 # ----------------------

hyppo/conditional/FCIT.py

@@ -3,10 +3,12 @@
 import joblib
 import numpy as np
 from scipy.stats import ttest_1samp
+from sklearn.base import clone
 from sklearn.metrics import mean_squared_error as mse
 from sklearn.model_selection import GridSearchCV, ShuffleSplit
 from sklearn.preprocessing import StandardScaler
 from sklearn.tree import DecisionTreeRegressor
+from sklearn.utils import check_random_state
 
 from .base import ConditionalIndependenceTest, ConditionalIndependenceTestOutput
 
@@ -30,6 +32,16 @@ class FCIT(ConditionalIndependenceTest):
         Proportion of data to evaluate test stat on.
     discrete: tuple of string
         Whether :math:`X` or :math:`Y` are discrete
+    random_state: int, RandomState instance, or None
+        Controls every source of randomness used internally (the data
+        permutations, the cross-validation splits, and, when ``model``
+        doesn't already pin its own ``random_state``, the regressor
+        fit on each permutation). Pass an int for reproducible results.
+        When ``None`` (the default), results are not reproducible and,
+        because the underlying regressor then draws from NumPy's global
+        random state, the exact statistic/p-value can also drift across
+        scikit-learn versions whenever they change how many times that
+        global state gets consumed internally.
 
     Notes
     -----
@@ -62,13 +74,15 @@ def __init__(
         num_perm=8,
         prop_test=0.1,
         discrete=(False, False),
+        random_state=None,
     ):
 
         self.model = model
         self.cv_grid = cv_grid
         self.num_perm = num_perm
         self.prop_test = prop_test
         self.discrete = discrete
+        self.random_state = random_state
         ConditionalIndependenceTest.__init__(self)
 
     def statistic(self, x, y, z=None):
@@ -91,16 +105,29 @@ def statistic(self, x, y, z=None):
         n_samples = x.shape[0]
         n_test = int(n_samples * self.prop_test)
 
+        rng = check_random_state(self.random_state)
+
         data_permutations = [
-            np.random.permutation(x.shape[0]) for i in range(self.num_perm)
+            rng.permutation(x.shape[0]) for i in range(self.num_perm)
         ]
+        reshuffle_seeds = rng.randint(0, np.iinfo(np.int32).max, size=self.num_perm)
 
-        clf = _cross_val(x, y, z, self.cv_grid, self.model, prop_test=self.prop_test)
+        cv_seed = rng.randint(0, np.iinfo(np.int32).max)
+        clf = _cross_val(
+            x,
+            y,
+            z,
+            self.cv_grid,
+            self.model,
+            prop_test=self.prop_test,
+            random_state=cv_seed,
+        )
         datadict = {
             "x": x,
             "y": y,
             "z": z,
             "data_permutation": data_permutations,
+            "reshuffle_seeds": reshuffle_seeds,
             "n_test": n_test,
             "reshuffle": False,
             "clf": clf,
@@ -113,12 +140,19 @@ def statistic(self, x, y, z=None):
         )
 
         if z.shape[1] == 0:
-            x_indep_y = x[np.random.permutation(n_samples)]
+            x_indep_y = x[rng.permutation(n_samples)]
         else:
             x_indep_y = np.empty([x.shape[0], 0])
 
+        cv_seed = rng.randint(0, np.iinfo(np.int32).max)
         clf = _cross_val(
-            x_indep_y, y, z, self.cv_grid, self.model, prop_test=self.prop_test
+            x_indep_y,
+            y,
+            z,
+            self.cv_grid,
+            self.model,
+            prop_test=self.prop_test,
+            random_state=cv_seed,
         )
 
         datadict["reshuffle"] = True
@@ -192,23 +226,30 @@ def test(self, x, y, z=None):
         return ConditionalIndependenceTestOutput(stat, pvalue)
 
 
-def _cross_val(x, y, z, cv_grid, model, prop_test):
+def _cross_val(x, y, z, cv_grid, model, prop_test, random_state=None):
     """
     Choose the regression hyperparameters by
     cross-validation.
     """
 
-    splitter = ShuffleSplit(n_splits=3, test_size=prop_test)
+    splitter = ShuffleSplit(n_splits=3, test_size=prop_test, random_state=random_state)
     cv = GridSearchCV(estimator=model, cv=splitter, param_grid=cv_grid, n_jobs=-1)
     cv.fit(_interleave(x, z), y)
 
-    return type(model)(**cv.best_params_)
+    best_model = clone(model)
+    best_model.set_params(**cv.best_params_)
+
+    params = best_model.get_params()
+    if params.get("random_state") is None and random_state is not None:
+        best_model.set_params(random_state=random_state)
+
+    return best_model
 
 
 def _interleave(x, z, seed=None):
     """Interleave x and z dimension-wise."""
     state = np.random.get_state()
-    np.random.seed(seed or int(time.time()))
+    np.random.seed(seed if seed is not None else int(time.time()))
     total_ids = np.random.permutation(x.shape[1] + z.shape[1])
     np.random.set_state(state)
     out = np.zeros([x.shape[0], x.shape[1] + z.shape[1]])
@@ -227,7 +268,12 @@ def _obtain_error(data_and_i):
     y = data["y"]
     z = data["z"]
     if data["reshuffle"]:
+        reshuffle_seeds = data.get("reshuffle_seeds")
+        seed = None if reshuffle_seeds is None else int(reshuffle_seeds[i])
+        state = np.random.get_state()
+        np.random.seed(seed)
         perm_ids = np.random.permutation(x.shape[0])
+        np.random.set_state(state)
     else:
         perm_ids = np.arange(x.shape[0])
     data_permutation = data["data_permutation"][i]

hyppo/conditional/_utils.py

@@ -6,13 +6,14 @@
 class _CheckInputs:
     """Checks inputs for all independence tests"""
 
-    def __init__(self, x, y, z, reps=None, max_dims=None, ignore_z_var=False):
+    def __init__(self, x, y, z=None, reps=None, max_dims=None, ignore_z_var=False):
         self.x = x
         self.y = y
         self.z = z
         self.reps = reps
         self.max_dims = max_dims
         self.ignore_z_var = ignore_z_var # to allow for constant z input
+        self.is_zero_variance = False
 
     def __call__(self):
         check_ndarray_xyz(self.x, self.y, self.z)
@@ -59,8 +60,8 @@ def check_dim_xyz(self, max_dims):
                 raise ValueError(
                     f"x, y, z must have be univariate and have shape [n,{max_dims}]"
                 )
-
         self._check_nd_indeptest()
+        self._check_variance()
 
         return self.x, self.y, self.z
 
@@ -86,10 +87,9 @@ def _check_min_samples(self):
 
     def _check_variance(self):
         if np.var(self.x) == 0:
-        # or np.var(self.y) == 0 or np.var(self.z) == 0:
-            raise ValueError("Test cannot be run. Input array x has 0 variance.")
+            self.is_zero_variance = True
         if np.var(self.y) == 0:
-            raise ValueError("Test cannot be run. Input array y has 0 variance")
-        if not self.ignore_z_var:
+            self.is_zero_variance = True
+        if not self.ignore_z_var and self.z is not None:
             if np.var(self.z) == 0:
-                raise ValueError("Test cannot be run. Input array z has 0 variance")
+                self.is_zero_variance = True

hyppo/conditional/cdcorr.py

@@ -172,6 +172,12 @@ def test(
         check_input = _CheckInputs(x, y, z, reps=reps, ignore_z_var=True)
         x, y, z = check_input()
 
+        if check_input.is_zero_variance:
+            self.stat = 0.0
+            self.pvalue = 1.0
+            self.null_dist = None
+            return ConditionalIndependenceTestOutput(0.0, 1.0)
+
         if not self.is_distance:
             x, y = compute_dist(x, y, metric=self.compute_distance, **self.kwargs)
             z = self._compute_kde(z)

hyppo/conditional/kci.py

@@ -112,6 +112,10 @@ def test(self, x, y):
 
         T = len(y)
 
+        if np.var(x) == 0 or np.var(y) == 0:
+            self.stat = 0.0
+            return ConditionalIndependenceTestOutput(0.0, 1.0)
+
         Kx, Ky = self.compute_kern(x, y)
         stat = self.statistic(x, y)
 

hyppo/conditional/pcorr.py

@@ -148,6 +148,12 @@ def test(
         check_input = _CheckInputs(x, y, z, reps=reps, max_dims=1)
         x, y, z = check_input()
 
+        if check_input.is_zero_variance:
+            self.stat = 0.0
+            self.pvalue = 1.0
+            self.null_dist = None
+            return ConditionalIndependenceTestOutput(0.0, 1.0)
+
         if auto:  # run t-stat
             stat = self.statistic(x, y, z)
 

hyppo/conditional/pdcorr.py

@@ -181,6 +181,12 @@ def test(
         check_input = _CheckInputs(x, y, z, reps=reps)
         x, y, z = check_input()
 
+        if check_input.is_zero_variance:
+            self.stat = 0.0
+            self.pvalue = 1.0
+            self.null_dist = None
+            return ConditionalIndependenceTestOutput(0.0, 1.0)
+
         if not self.is_distance:
             x, y, z = compute_dist(x, y, z, metric=self.compute_distance, **self.kwargs)
             self.is_distance = True

hyppo/conditional/tests/test_FCIT.py

@@ -10,22 +10,20 @@ class TestFCIT:
     @pytest.mark.parametrize(
         "n, obs_stat, obs_pvalue",
         [
-            (2000, 11.677197, 3.8168e-06),
+            (2000, 15.854040, 4.815497e-07),
         ],
     )
     def test_linear_oned(self, n, obs_stat, obs_pvalue):
         np.random.seed(123456789)
         x, y = rot_ksamp("linear", n, 1, k=2)
-        np.random.seed(123456789)
-        stat, pvalue = FCIT().test(x, y)
+        stat, pvalue = FCIT(random_state=0).test(x, y)
 
         assert_almost_equal(stat, obs_stat, decimal=-1)
         assert_almost_equal(pvalue, obs_pvalue, decimal=4)
 
     @pytest.mark.parametrize(
         "dim, n, obs_stat, obs_pvalue",
-        # 0.56139, -0.16024
-        [(1, 100000, -0.06757, 0.52599), (2, 100000, -4.59882, 0.99876)],
+        [(1, 100000, -1.339294, 0.888837), (2, 100000, -6.985438, 0.999893)],
     )
     def test_null(self, dim, n, obs_stat, obs_pvalue):
         np.random.seed(12)
@@ -47,18 +45,18 @@ def test_null(self, dim, n, obs_stat, obs_pvalue):
             ).T
         )
 
-        np.random.seed(122)
-        stat, pvalue = FCIT().test(x1.T, y1.T, z1)
+        # random_state pins every source of randomness FCIT uses internally,
+        # so this is reproducible across scikit-learn versions (gh-427).
+        stat, pvalue = FCIT(random_state=0).test(x1.T, y1.T, z1)
 
         assert_almost_equal(pvalue, obs_pvalue, decimal=4)
         assert_almost_equal(stat, obs_stat, decimal=4)
 
     @pytest.mark.parametrize(
         "dim, n, obs_stat, obs_pvalue",
         [
-            #89.271754, 161.35165
-            (1, 100000, 89.184784, 2.91447597e-12),
-            (2, 100000, 161.35105, 4.63412957e-14),
+            (1, 100000, 100.371119, 1.283942e-12),
+            (2, 100000, 187.502936, 1.619644e-14),
         ],
     )
     def test_alternative(self, dim, n, obs_stat, obs_pvalue):
@@ -83,8 +81,7 @@ def test_alternative(self, dim, n, obs_stat, obs_pvalue):
             ).T
         )
 
-        np.random.seed(122)
-        stat, pvalue = FCIT().test(x2.T, y2.T, z2)
+        stat, pvalue = FCIT(random_state=0).test(x2.T, y2.T, z2)
 
         assert_almost_equal(pvalue, obs_pvalue, decimal=12)
         assert_almost_equal(stat, obs_stat, decimal=4)

hyppo/conditional/tests/test_utils.py

@@ -64,7 +64,9 @@ def test_constant_input(self):
         y = np.arange(20).reshape(-1, 1)
         z = np.ones(20).reshape(-1, 1)
 
-        assert_raises(ValueError, _CheckInputs(x, y, z, ignore_z_var=False))
+        check_input = _CheckInputs(x, y, z, ignore_z_var=False)
+        check_input()
+        assert check_input.is_zero_variance == True
 
         try:
             _CheckInputs(x, y, z, ignore_z_var=True)