Repository: jakobrunge/tigramite (language: python)

Issue:
Tom branch 1
Implemented simple pairwise method and renamed parameters and changed documentation accordingly


Discussion / hints:

Fix patch:
diff --git a/tutorials/causal_discovery/tigramite_tutorial_pairwise_mult_ci.ipynb b/tutorials/causal_discovery/tigramite_tutorial_pairwise_mult_ci.ipynb
index 8466d094..b3c4af7b 100644
--- a/tutorials/causal_discovery/tigramite_tutorial_pairwise_mult_ci.ipynb
+++ b/tutorials/causal_discovery/tigramite_tutorial_pairwise_mult_ci.ipynb
@@ -45,10 +45,10 @@
     "Suppose one has random vectors $X, Y, Z$ (each potentially multivariate) and suppose that one  wants to decide whether $X$ is independent of $Y$ given $Z$ (that is, the conditional densities (assuming they exist) factorize: $\\forall x,y,z:\\;p_{X,Y|Z=z}(x,y) = p_{X|Z=z}(x)\\cdot p_{Y|Z=z}(y)$.)\n",
     "\n",
     "One approach to test independence between $X$ and $Y$ given $Z$ is to test whether all components $X_i$ and $Y_j$ are independent given $Z$.\n",
-    "The conditional independence test PairwiseMultCI builds on this pairwise idea. PairwiseMultCI tests independence between each pair of components $X_i$ and $Y_j$ given $(Z,S_{ij}).$ Here, $S_{ij}$ is a set that consists of components of either $X$ or $Y$ that are independent of $Y_j$ or $X_i$. One can show that including these components increases the effect size of the corresponding test.\n",
+    "The conditional independence test PairwiseMultCI builds on this pairwise idea. PairwiseMultCI tests independence between each pair of components $X_i$ and $Y_j$ given $(Z,S_{ij}).$ Here, $S_{ij}$ is a set that consists of components of either $X$ or $Y$ that are independent of $Y_j$ or $X_i$. One can show that including these components increases the effect size of the corresponding test. To choose between the pairwise independence testing with and without increased conditioning sets, there is the boolean parameter ``learn_augmented_cond_sets``. Setting it to False will lead to pairwise independence testing without increased conditioning sets, setting it to true will lead to increased conditioning sets.\n",
     "\n",
     "In practice, these sets $S_{ij}$ are usually not known a priori and hence need to be learned.\n",
-    "PairwiseMultCI first learns these conditional independencies on one part of the sample and then uses these conditional independencies to do the testing on the second part of the sample. The sample splitting ratio (``pre_step_sample_fraction``) is a hyperparameter, as is the significance level ($\\alpha_{pre}$) for learning independencies in the first step. (We think that only a small part of the sample should be used for learning independencies in the first step, moreover, $\\alpha_{pre}$ should be rather large to only mark something as independent if there is strong evidence).\n",
+    "PairwiseMultCI first learns these conditional independencies on one part of the sample and then uses these conditional independencies to do the testing on the second part of the sample. The sample splitting ratio (``pre_step_sample_fraction``) is a hyperparameter, as is the significance level ($\\alpha_{pre}$) for learning independencies in the first step. (We think that only a small part of the sample should be used for learning independencies in the first step, moreover, $\\alpha_{pre}$ should be rather large to only mark something as independent if there is strong evidence). As a small side remark, setting ``pre_step_sample_fraction`` to $0$ is the same as setting ``learn_augmented_cond_sets`` to False.\n",
     "\n",
     "We remark that PairwiseMultCI is a \"meta-algorithm\", that is, it is a framework that works with several univariate test statistics. Because of that, PairwiseMultCI also has the parameter cond_ind_test, which can be any other conditional independence test from tigramite/independence_tests"
    ]
@@ -92,8 +92,8 @@
      "name": "stdout",
      "output_type": "stream",
      "text": [
-      "test_stat 0.8330221744959563\n",
-      "pval 6.30177446825862e-21\n"
+      "test_stat 0.8011087470057225\n",
+      "pval 2.584051117524395e-18\n"
      ]
     }
    ],
@@ -101,7 +101,7 @@
     "# alpha_pre: significance level for the first step\n",
     "# pre_step_sample_fraction: relative size of the sample for the first step\n",
     "# cond_ind_test: respective univariate conditional independence test\n",
-    "ci_test = PairwiseMultCI(alpha_pre = 0.5, pre_step_sample_fraction = 0.2, cond_ind_test = ParCorr())\n",
+    "ci_test = PairwiseMultCI(learn_augmented_cond_sets = True, alpha_pre = 0.5, pre_step_sample_fraction = 0.2, cond_ind_test = ParCorr())\n",
     "test_stat, pval = ci_test.run_test_raw(x = x, y = y, z = z)\n",
     "print(\"test_stat\", test_stat)\n",
     "print(\"pval\", pval)"
@@ -154,16 +154,16 @@
     "    y2 = y1 + 0.3 * np.random.normal(0, 1, T). reshape(T, 1) + 0.1 * x\n",
     "    y = np.hstack((y1, y2))\n",
     "    # PairwiseMultCI\n",
-    "    ci_test = PairwiseMultCI(alpha_pre = 0.5, pre_step_sample_fraction = 0.2, cond_ind_test = par_corr)\n",
-    "    test_stat, pval = ci_test.run_test_raw(x = x, y = y, z = z)\n",
-    "    if (pval <= alpha):\n",
+    "    ci_test = PairwiseMultCI(learn_augmented_cond_sets = True, alpha_pre = 0.5, pre_step_sample_fraction = 0.2, cond_ind_test = par_corr)\n",
+    "    test_stat1, pval1 = ci_test.run_test_raw(x = x, y = y, z = z)\n",
+    "    if (pval1 <= alpha):\n",
     "        rejection_bool[i] = 1  \n",
     "           \n",
     "    ## competitor:\n",
-    "    pval1 = par_corr.run_test_raw(x, y1, z)[1]\n",
-    "    pval2 = par_corr.run_test_raw(x, y2, z)[1]  \n",
-    "    if (np.min((pval1, pval2)) * 2 <= alpha):\n",
-    "        rejection_bool_competitor[i] = 1 \n",
+    "    ci_test_competitor = PairwiseMultCI()\n",
+    "    test_stat2, pval2 = ci_test_competitor.run_test_raw(x = x, y = y, z = z)\n",
+    "    if (pval2 <= alpha):\n",
+    "        rejection_bool_competitor[i] = 1  \n",
     "print(\"Rejection Rate PairwiseMultCI:\", np.mean(rejection_bool))\n",
     "print(\"Rejection Rate Competitor:\",np.mean(rejection_bool_competitor))"
    ]
@@ -213,16 +213,16 @@
     "    y2 = 0.3 * np.random.normal(0, 1, T). reshape(T, 1) + 0.1 * x\n",
     "    y = np.hstack((y1, y2))\n",
     "    # PairwiseMultCI\n",
-    "    ci_test = PairwiseMultCI(alpha_pre = 0.5, pre_step_sample_fraction = 0.2, cond_ind_test = par_corr)\n",
-    "    test_stat, pval = ci_test.run_test_raw(x = x, y = y, z = z)\n",
-    "    if (pval <= alpha):\n",
+    "    ci_test = PairwiseMultCI(learn_augmented_cond_sets = True, alpha_pre = 0.5, pre_step_sample_fraction = 0.2, cond_ind_test = par_corr)\n",
+    "    test_stat1, pval1 = ci_test.run_test_raw(x = x, y = y, z = z)\n",
+    "    if (pval1 <= alpha):\n",
     "        rejection_bool[i] = 1  \n",
     "           \n",
     "    ## competitor:\n",
-    "    pval1 = par_corr.run_test_raw(x, y1, z)[1]\n",
-    "    pval2 = par_corr.run_test_raw(x, y2, z)[1]  \n",
-    "    if (np.min((pval1, pval2)) * 2 <= alpha):\n",
-    "        rejection_bool_competitor[i] = 1 \n",
+    "    ci_test_competitor = PairwiseMultCI()\n",
+    "    test_stat2, pval2 = ci_test_competitor.run_test_raw(x = x, y = y, z = z)\n",
+    "    if (pval2 <= alpha):\n",
+    "        rejection_bool_competitor[i] = 1  \n",
     "print(\"Rejection Rate PairwiseMultCI:\", np.mean(rejection_bool))\n",
     "print(\"Rejection Rate Competitor:\",np.mean(rejection_bool_competitor))"
    ]


Test patch:
diff --git a/tests/test_PairwiseMultCI.py b/tests/test_PairwiseMultCI.py
index 501fb0d9..cc3db3fa 100644
--- a/tests/test_PairwiseMultCI.py
+++ b/tests/test_PairwiseMultCI.py
@@ -15,43 +15,49 @@
 @pytest.fixture(params=[
     # Generate PairwiseMultCI test instances
     # test PairwiseMultCI for different combinations of:
-    # significance (of the univariate tests),  alpha_pre, pre_step_sample_fraction, cond_ind_test, cond_ind_test_thres, cond_ind_test_thres_pre
-    ('fixed_thres', 0.5, 0.2, ParCorr, 1, 0.5),
-    ('fixed_thres', 0.5, 0.5, ParCorr, 1, 0.5),
-    ('fixed_thres', 0.8, 0.2, ParCorr, 1, 0.5),
-    ('fixed_thres', 0.8, 0.5, ParCorr, 1, 0.5),
-    ('analytic', 0.5, 0.2, ParCorr, None, None),
-    ('shuffle_test', 0.5, 0.2, ParCorr, None, None),
-    ('analytic', 0.5, 0.5, ParCorr, None, None),
-    ('shuffle_test', 0.5, 0.5, ParCorr, None, None),
-    ('analytic', 0.8, 0.2, ParCorr, None, None),
-    ('shuffle_test', 0.8, 0.2, ParCorr, None, None),
-    ('analytic', 0.8, 0.5, ParCorr, None, None),
-    ('shuffle_test', 0.8, 0.5, ParCorr, None, None),
-    ('analytic', 0.5, 0.2, RobustParCorr, None, None),
-    ('shuffle_test', 0.5, 0.2, RobustParCorr, None, None),
-    ('analytic', 0.5, 0.5, RobustParCorr, None, None),
-    ('shuffle_test', 0.5, 0.5, RobustParCorr, None, None),
-    ('analytic', 0.8, 0.2, RobustParCorr, None, None),
-    ('shuffle_test', 0.8, 0.2, RobustParCorr, None, None),
-    ('analytic', 0.8, 0.5, RobustParCorr, None, None),
-    ('shuffle_test', 0.8, 0.5, RobustParCorr, None, None),
-    ('shuffle_test', 0.5, 0.2, CMIknn, None, None),
-    ('shuffle_test', 0.5, 0.5, CMIknn, None, None),
-    ('shuffle_test', 0.8, 0.2, CMIknn, None, None),
-    ('shuffle_test', 0.8, 0.5, CMIknn, None, None),
+    # learn_augmented_cond_sets, significance (of the univariate tests),  alpha_pre, pre_step_sample_fraction, cond_ind_test, cond_ind_test_thres, cond_ind_test_thres_pre
+    ('True', 'fixed_thres', 0.5, 0.2, ParCorr, 1, 0.5),
+    ('True', 'fixed_thres', 0.5, 0.5, ParCorr, 1, 0.5),
+    ('True', 'fixed_thres', 0.8, 0.2, ParCorr, 1, 0.5),
+    ('True', 'fixed_thres', 0.8, 0.5, ParCorr, 1, 0.5),
+    ('True', 'analytic', 0.5, 0.2, ParCorr, None, None),
+    ('True', 'shuffle_test', 0.5, 0.2, ParCorr, None, None),
+    ('True', 'analytic', 0.5, 0.5, ParCorr, None, None),
+    ('True', 'shuffle_test', 0.5, 0.5, ParCorr, None, None),
+    ('True', 'analytic', 0.8, 0.2, ParCorr, None, None),
+    ('True', 'shuffle_test', 0.8, 0.2, ParCorr, None, None),
+    ('True', 'analytic', 0.8, 0.5, ParCorr, None, None),
+    ('True', 'shuffle_test', 0.8, 0.5, ParCorr, None, None),
+    ('True', 'analytic', 0.5, 0.2, RobustParCorr, None, None),
+    ('True', 'shuffle_test', 0.5, 0.2, RobustParCorr, None, None),
+    ('True', 'analytic', 0.5, 0.5, RobustParCorr, None, None),
+    ('True', 'shuffle_test', 0.5, 0.5, RobustParCorr, None, None),
+    ('True', 'analytic', 0.8, 0.2, RobustParCorr, None, None),
+    ('True', 'shuffle_test', 0.8, 0.2, RobustParCorr, None, None),
+    ('True', 'analytic', 0.8, 0.5, RobustParCorr, None, None),
+    ('True', 'shuffle_test', 0.8, 0.5, RobustParCorr, None, None),
+    ('True', 'shuffle_test', 0.5, 0.2, CMIknn, None, None),
+    ('True', 'shuffle_test', 0.5, 0.5, CMIknn, None, None),
+    ('True', 'shuffle_test', 0.8, 0.2, CMIknn, None, None),
+    ('True', 'shuffle_test', 0.8, 0.5, CMIknn, None, None),
+    ('False', 'fixed_thres', 0.0, 0.0, ParCorr, 1, 0.5),
+    ('False', 'analytic', 0.0, 0.0, ParCorr, None, None),
+    ('False', 'shuffle_test', 0.0, 0.0, ParCorr, None, None),
+    ('False', 'analytic', 0.0, 0.0, RobustParCorr, None, None),
+    ('False', 'shuffle_test', 0.0, 0.0, RobustParCorr, None, None),
+    ('False', 'shuffle_test', 0.0, 0.0, CMIknn, None, None),
 ])
 
 def pairwise_mult_ci(request):
     # Unpack the parameters
-    sig, alpha_pre, pre_step_sample_fraction, cond_ind_test, cond_ind_test_thres, cond_ind_test_thres_pre = request.param
+    learn_augmented_cond_sets, sig, alpha_pre, pre_step_sample_fraction, cond_ind_test, cond_ind_test_thres, cond_ind_test_thres_pre = request.param
     # Generate the par_corr_wls independence test
     if sig != "fixed_thres":
-        return PairwiseMultCI(cond_ind_test = cond_ind_test(significance = sig),
+        return PairwiseMultCI(learn_augmented_cond_sets = learn_augmented_cond_sets, cond_ind_test = cond_ind_test(significance = sig),
                           alpha_pre = alpha_pre,
                           pre_step_sample_fraction = pre_step_sample_fraction)
     else:
-        return PairwiseMultCI(cond_ind_test = cond_ind_test(significance = sig),
+        return PairwiseMultCI(learn_augmented_cond_sets = learn_augmented_cond_sets, cond_ind_test = cond_ind_test(significance = sig),
                           alpha_pre = None,
                           pre_step_sample_fraction=pre_step_sample_fraction,
                           significance= sig,
diff --git a/tigramite/independence_tests/pairwise_CI.py b/tigramite/independence_tests/pairwise_CI.py
index 222d6b9c..66aac9cf 100644
--- a/tigramite/independence_tests/pairwise_CI.py
+++ b/tigramite/independence_tests/pairwise_CI.py
@@ -38,6 +38,13 @@ class PairwiseMultCI(CondIndTest):
 
     Parameters
     ----------
+
+    learn_augmented_cond_sets: boolean
+       Boolean on whether PairwiseMultCI should be executed with a pre step that learns conditional independencies
+       and then augments conditioning sets.
+        That is explained in the tutorial or the above mentioned paper in more detail.
+        False indicates that conditioning sets should not be increased
+
     alpha_pre: float
         Significance level for the first step of the algorithm.
         If cond_ind_test is instantiated with significance = "fixed_thres",
@@ -70,13 +77,14 @@ def measure(self):
         """
         return self._measure
 
-    def __init__(self, cond_ind_test = ParCorr(), alpha_pre = 0.5, pre_step_sample_fraction = 0.2, fixed_thres_pre = None, **kwargs):
+    def __init__(self, learn_augmented_cond_sets = False, cond_ind_test = ParCorr(), alpha_pre = 0.5, pre_step_sample_fraction = 0.2, fixed_thres_pre = None, **kwargs):
         self._measure = 'pairwise_CI'
         self.cond_ind_test = cond_ind_test
         self.alpha_pre = alpha_pre
         self.pre_step_sample_fraction = pre_step_sample_fraction
         self.fixed_thres_pre = fixed_thres_pre
         self.two_sided = False
+        self.learn_augmented_cond_sets = learn_augmented_cond_sets
         CondIndTest.__init__(self, **kwargs)
 
 
@@ -118,6 +126,9 @@ def calculate_dep_measure_and_significance(self, array, xyz, data_type = None, c
             raise ValueError("If significance == 'fixed_thres', fixed_thres_pre for the"
                              " pre-step needs to be defined in initializing PairwiseMultCI.")
 
+        if self.learn_augmented_cond_sets == False:
+            self.pre_step_sample_fraction = 0
+
         x_indices = np.where(xyz == 0)[0]
         y_indices = np.where(xyz == 1)[0]
         z_indices = np.where(xyz == 2)[0]
@@ -130,6 +141,10 @@ def calculate_dep_measure_and_significance(self, array, xyz, data_type = None, c
         size_first_block = np.floor(self.pre_step_sample_fraction * T)
         size_first_block = size_first_block.astype(int)
 
+        if size_first_block < 10:
+            self.pre_step_sample_fraction = 0
+
+
         # split the sample
         array_s1 = array[:, 0:size_first_block]
         x_s1 = array_s1[x_indices]
@@ -169,31 +184,34 @@ def calculate_dep_measure_and_significance(self, array, xyz, data_type = None, c
             z_type_s2 = None
 
         ## Step 1: estimate conditional independencies
-        if fixed_thres_bool == False:
-            p_vals_pre = np.zeros((dim_x, dim_y))
-        else:
-            vals_pre = np.zeros((dim_x, dim_y))
-        for j in np.arange(0, dim_x):
-            for jj in np.arange(0, dim_y):
-                if fixed_thres_bool == False:
-                    p_vals_pre[j, jj] = self.cond_ind_test.run_test_raw(x_s1[j].reshape(size_first_block, 1),
-                                                                   y_s1[jj].reshape(size_first_block, 1),
-                                                                   z_s1.reshape(size_first_block, dim_z),
-                                                                   x_type = x_type_s1,
-                                                                   y_type = y_type_s1,
-                                                                   z_type = z_type_s1)[1]
-                else:
-                    vals_pre[j, jj] = self.cond_ind_test.run_test_raw(x_s1[j].reshape(size_first_block, 1),
-                                                                   y_s1[jj].reshape(size_first_block, 1),
-                                                                   z_s1.reshape(size_first_block, dim_z),
-                                                                   x_type = x_type_s1,
-                                                                   y_type = y_type_s1,
-                                                                   z_type = z_type_s1,
-                                                                   alpha_or_thres=999.)[0] # just a dummy
-        if fixed_thres_bool == False:
-            indep_set = np.where(p_vals_pre > self.alpha_pre)
+        if self.pre_step_sample_fraction > 0:
+            if fixed_thres_bool == False:
+                p_vals_pre = np.zeros((dim_x, dim_y))
+            else:
+                vals_pre = np.zeros((dim_x, dim_y))
+            for j in np.arange(0, dim_x):
+                for jj in np.arange(0, dim_y):
+                    if fixed_thres_bool == False:
+                        p_vals_pre[j, jj] = self.cond_ind_test.run_test_raw(x_s1[j].reshape(size_first_block, 1),
+                                                                       y_s1[jj].reshape(size_first_block, 1),
+                                                                       z_s1.reshape(size_first_block, dim_z),
+                                                                       x_type = x_type_s1,
+                                                                       y_type = y_type_s1,
+                                                                       z_type = z_type_s1)[1]
+                    else:
+                        vals_pre[j, jj] = self.cond_ind_test.run_test_raw(x_s1[j].reshape(size_first_block, 1),
+                                                                       y_s1[jj].reshape(size_first_block, 1),
+                                                                       z_s1.reshape(size_first_block, dim_z),
+                                                                       x_type = x_type_s1,
+                                                                       y_type = y_type_s1,
+                                                                       z_type = z_type_s1,
+                                                                       alpha_or_thres=999.)[0] # just a dummy
+            if fixed_thres_bool == False:
+                indep_set = np.where(p_vals_pre > self.alpha_pre)
+            else:
+                indep_set = np.where(np.abs(vals_pre) >= self.fixed_thres_pre)
         else:
-            indep_set = np.where(np.abs(vals_pre) >= self.fixed_thres_pre)
+            indep_set = np.where(np.zeros((dim_x, dim_y)) > np.zeros((dim_x, dim_y)))
 
         # Step 2: test conditional independencies with increased effect sizes
         if self.cond_ind_test.significance != "fixed_thres":
@@ -308,7 +326,7 @@ def get_analytic_significance(self, value = None, T = None, dim = None, xyz = No
     from tigramite.independence_tests.robust_parcorr import RobustParCorr
     from tigramite.independence_tests.cmiknn import CMIknn
     alpha = 0.05
-    ci = PairwiseMultCI(fixed_thres_pre = 0.5, cond_ind_test = ParCorr(significance = "fixed_thres"))
+    ci = PairwiseMultCI(learn_augmented_cond_sets = True, fixed_thres_pre = 0.5, cond_ind_test = RobustParCorr(significance = "fixed_thres"))
     # ci = PairwiseMultCI(cond_ind_test=ParCorr(significance="fixed_thres"))
     # ci = PairwiseMultCI(cond_ind_test=ParCorr(significance="analytic"))
     T = 100