Merge pull request #721 from DarthTrevis/permutations

Permutations

Author: Sebastian Raschka

docs/sources/CHANGELOG.md

@@ -23,7 +23,7 @@ The CHANGELOG for the current development version is available at
 
 ##### Changes
 
-- -
+- `permutation_test` (`mlxtend.evaluate.permutation`) ìs corrected to give the proportion of permutations whose statistic is *at least as extreme* as the one observed. ([#721](https://github.com/rasbt/mlxtend/pull/721) via [Florian Charlier](https://github.com/DarthTrevis))
 
 ##### Bug Fixes
 

docs/sources/user_guide/evaluate/permutation_test.ipynb

@@ -43,13 +43,13 @@
     "4. Divide the permuted dataset into two datasets x' and y' of size *n* and *m*, respectively\n",
     "5. Compute the difference (here: mean) of sample x' and sample y' and record this difference\n",
     "6. Repeat steps 3-5 until all permutations are evaluated\n",
-    "7. Return the p-value as the number of times the recorded differences were more extreme than the original difference from 1. and divide this number by the total number of permutations\n",
+    "7. Return the p-value as the number of times the recorded differences were at least as extreme as the original difference from 1. and divide this number by the total number of permutations\n",
     "\n",
     "Here, the p-value is defined as the probability, given the null hypothesis (no difference between the samples) is true, that we obtain results that are at least as extreme as the results we observed (i.e., the sample difference from 1.).\n",
     "\n",
-    "More formally, we can express the computation of the p-value as follows ([2]):\n",
+    "More formally, we can express the computation of the p-value as follows (adapted from [2]):\n",
     "\n",
-    "$$p(t > t_0) = \\frac{1}{(n+m)!} \\sum^{(n+m)!}_{j=1} I(t_j > t_0),$$\n",
+    "$$p(t \\geq t_0) = \\frac{1}{(n+m)!} \\sum^{(n+m)!}_{j=1} I(t_j \\geq t_0),$$\n",
     "\n",
     "where $t_0$ is the observed value of the test statistic (1. in the list above), and $t$ is the t-value, the statistic computed from the resamples (5.) $t(x'_1, x'_2, ..., x'_n, y'_1, y'_2, ..., x'_m) = |\\bar{x'} - \\bar{y'}|$, and *I* is the indicator function.\n",
     "\n",
@@ -114,7 +114,7 @@
      "name": "stdout",
      "output_type": "stream",
      "text": [
-      "0.0066\n"
+      "0.0066993300669933005\n"
      ]
     }
    ],
@@ -159,7 +159,7 @@
      "output_type": "stream",
      "text": [
       "Observed pearson R: 0.81\n",
-      "P value: 0.09\n"
+      "P value: 0.10\n"
      ]
     }
    ],
@@ -281,7 +281,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.6.4"
+   "version": "3.8.5"
   },
   "toc": {
    "nav_menu": {},

mlxtend/evaluate/permutation.py

@@ -97,7 +97,7 @@ def func(x, y):
     m, n = len(x), len(y)
     combined = np.hstack((x, y))
 
-    more_extreme = 0.
+    at_least_as_extreme = 0.
     reference_stat = func(x, y)
 
     # Note that whether we compute the combinations or permutations
@@ -120,15 +120,21 @@ def func(x, y):
             indices_y = [i for i in range(m + n) if i not in indices_x]
             diff = func(combined[list(indices_x)], combined[indices_y])
 
-            if diff > reference_stat:
-                more_extreme += 1.
+            if diff > reference_stat or np.isclose(diff, reference_stat):
+                at_least_as_extreme += 1.
 
         num_rounds = factorial(m + n) / (factorial(m)*factorial(n))
 
     else:
         for i in range(num_rounds):
             rng.shuffle(combined)
-            if func(combined[:m], combined[m:]) > reference_stat:
-                more_extreme += 1.
+            diff = func(combined[:m], combined[m:])
 
-    return more_extreme / num_rounds
+            if diff > reference_stat or np.isclose(diff, reference_stat):
+                at_least_as_extreme += 1.
+
+        # To cover the actual experiment results
+        at_least_as_extreme += 1
+        num_rounds += 1
+
+    return at_least_as_extreme / num_rounds

mlxtend/evaluate/tests/test_permutation.py

@@ -16,36 +16,36 @@
 def test_one_sided_x_greater_y():
     p = permutation_test(treatment, control,
                          func=lambda x, y: np.mean(x) - np.mean(y))
-    assert round(p, 4) == 0.0274, p
+    assert round(p, 4) == 0.0301, p
 
     p = permutation_test(treatment, control,
                          func="x_mean > y_mean")
-    assert round(p, 4) == 0.0274, p
+    assert round(p, 4) == 0.0301, p
 
 
 def test_one_sided_y_greater_x():
     p = permutation_test(treatment, control,
                          func=lambda x, y: np.mean(y) - np.mean(x))
-    assert round(p, 3) == 1 - 0.03, p
+    assert round(p, 3) == 0.973, p
 
     p = permutation_test(treatment, control,
                          func="x_mean < y_mean")
-    assert round(p, 3) == 1 - 0.03, p
+    assert round(p, 3) == 0.973, p
 
 
 def test_two_sided():
     p = permutation_test(treatment, control,
                          func=lambda x, y: np.abs(np.mean(x) - np.mean(y)))
-    assert round(p, 3) == 0.055, p
+    assert round(p, 3) == 0.060, p
 
     p = permutation_test(treatment, control,
                          func="x_mean != y_mean")
-    assert round(p, 3) == 0.055, p
+    assert round(p, 3) == 0.060, p
 
 
 def test_default():
     p = permutation_test(treatment, control)
-    assert round(p, 3) == 0.055, p
+    assert round(p, 3) == 0.060, p
 
 
 def test_approximateone_sided_x_greater_y():
@@ -55,7 +55,7 @@ def test_approximateone_sided_x_greater_y():
                          method='approximate',
                          num_rounds=5000,
                          seed=123)
-    assert round(p, 3) == 0.028, p
+    assert round(p, 3) == 0.031, p
 
 
 def test_invalid_method():

mlxtend/feature_extraction/tests/test_principal_component_analysis.py

@@ -108,7 +108,7 @@ def test_fail_array_dimension_2():
 def test_variance_explained_ratio():
     pca = PCA()
     pca.fit(X_std)
-    assert np.sum(pca.e_vals_normalized_) == 1.
+    assert_almost_equal(np.sum(pca.e_vals_normalized_), 1.)
     assert np.sum(pca.e_vals_normalized_ < 0.) == 0