Merge f51e028 into 894dce8

mlxtend/evaluate/permutation.py

@@ -7,7 +7,7 @@
 # License: BSD 3 clause
 
 import numpy as np
-from itertools import combinations
+from itertools import chain, combinations
 from math import factorial
 try:
     from nose.tools import nottest
@@ -17,11 +17,29 @@ def nottest(f):
         return f
 
 
+def all_combinations(m, n, paired=False):
+    if paired:
+        if m != n:
+            raise ValueError("Populations must have same size in paired test")
+            powerset_iter = chain.from_iterable(
+                combinations(range(n), r) for r in range(n)
+            )
+            for indices_1 in powerset_iter:
+                indices_2 = [i for i in range(n) if i not in indices_1]
+                indices_x = list(indices_1) + [i + n for i in indices_2]
+                indices_y = indices_2 + [i + n for i in indices_1]
+                yield indices_x, indices_y
+    else:
+        for indices_x in combinations(range(m + n), m):
+            indices_y = [i for i in range(m + n) if i not in indices_x]
+            yield list(indices_x), indices_y
+
+
 # decorator to prevent nose to consider
 # this as a unit test due to "test" in the name
 @nottest
 def permutation_test(x, y, func='x_mean != y_mean', method='exact',
-                     num_rounds=1000, seed=None):
+                     num_rounds=1000, seed=None, paired=False):
     """
     Nonparametric permutation test
 
@@ -50,6 +68,9 @@ def permutation_test(x, y, func='x_mean != y_mean', method='exact',
         given by `num_rounds`.
         Note that 'exact' is typically not feasible unless the dataset
         size is relatively small.
+    paired : bool
+        If True, a paired test is performed by only exchanging each
+        datapoint with its associate.
     num_rounds : int (default: 1000)
         The number of permutation samples if `method='approximate'`.
     seed : int or None (default: None)
@@ -114,21 +135,29 @@ def func(x, y):
     # and using combinations instead of permutations simply saves computational
     # time
 
-    if method == 'exact':
-        for indices_x in combinations(range(m + n), m):
-
-            indices_y = [i for i in range(m + n) if i not in indices_x]
+    if method == "exact":
+        for indices_x, indices_y in all_combinations(m, n, paired=paired):
             diff = func(combined[list(indices_x)], combined[indices_y])
 
             if diff > reference_stat or np.isclose(diff, reference_stat):
-                at_least_as_extreme += 1.
+                at_least_as_extreme += 1.0
 
-        num_rounds = factorial(m + n) / (factorial(m)*factorial(n))
+        if paired:
+            num_rounds = 2 ** n
+        else:
+            num_rounds = factorial(m + n) / (factorial(m) * factorial(n))
 
     else:
         for i in range(num_rounds):
-            rng.shuffle(combined)
-            diff = func(combined[:m], combined[m:])
+            if paired:
+                indices_1 = rng.randn(n) > 0
+                indices_2 = [i for i in range(n) if i not in indices_1]
+                indices_x = list(indices_1) + [i + n for i in indices_2]
+                indices_y = indices_2 + [i + n for i in indices_1]
+                diff = func(combined[indices_x], combined[indices_y])
+            else:
+                rng.shuffle(combined)
+                diff = func(combined[:m], combined[m:])
 
             if diff > reference_stat or np.isclose(diff, reference_stat):
                 at_least_as_extreme += 1.

mlxtend/evaluate/tests/test_permutation.py

@@ -79,3 +79,20 @@ def test_invalid_func():
                   [1, 2, 3],
                   [3, 4, 5],
                   'myfunc')
+
+
+def test_paired_runs():
+    permutation_test([1, 2, 3], [4, 5, 6], paired=True)
+
+
+def test_paired_runs_approximate():
+    permutation_test(
+        [1, 2, 3], [4, 5, 6], paired=True, method="approximate", num_rounds=10
+    )
+
+
+def test_paired_invalid_lengths():
+    msg = "Populations must have same size in paired test"
+    assert_raises(
+        ValueError, msg, permutation_test, [1, 2, 3], [3, 4], paired=True
+    )