[med-svn] [python-mne] 68/376: ENH : adding exact permutation test with all permutations + some tests and cosmits

[Yaroslav Halchenko]

This is an automated email from the git hooks/post-receive script.

yoh pushed a commit to annotated tag v0.1
in repository python-mne.

commit 71702c1366d14e8731e58f1fe874fdebadf93297
Author: Alexandre Gramfort <alexandre.gramfort at inria.fr>
Date:   Sun Feb 13 22:08:28 2011 -0500

    ENH : adding exact permutation test with all permutations
    + some tests and cosmits
---
 mne/stats/permutations.py            | 84 ++++++++++++++++++++++++++++++++----
 mne/stats/tests/test_permutations.py |  9 +++-
 2 files changed, 83 insertions(+), 10 deletions(-)

diff --git a/mne/stats/permutations.py b/mne/stats/permutations.py
index 7286304..6667962 100644
--- a/mne/stats/permutations.py
+++ b/mne/stats/permutations.py
@@ -1,7 +1,8 @@
 """T-test with permutations
 """
 
-# Author: Alexandre Gramfort <gramfort at nmr.mgh.harvard.edu>
+# Authors: Alexandre Gramfort <gramfort at nmr.mgh.harvard.edu>
+#          Fernando Perez (bin_perm_rep function)
 #
 # License: Simplified BSD
 
@@ -9,6 +10,41 @@ from math import sqrt
 import numpy as np
 
 
+def bin_perm_rep(ndim, a=0, b=1):
+    """bin_perm_rep(ndim) -> ndim permutations with repetitions of (a,b).
+
+    Returns an array with all the possible permutations with repetitions of
+    (0,1) in ndim dimensions.  The array is shaped as (2**ndim,ndim), and is
+    ordered with the last index changing fastest.  For examble, for ndim=3:
+
+    Examples:
+
+    >>> bin_perm_rep(3)
+    array([[0, 0, 0],
+           [0, 0, 1],
+           [0, 1, 0],
+           [0, 1, 1],
+           [1, 0, 0],
+           [1, 0, 1],
+           [1, 1, 0],
+           [1, 1, 1]])
+    """
+
+    # Create the leftmost column as 0,0,...,1,1,...
+    nperms = 2**ndim
+    perms = np.empty((nperms, ndim), type(a))
+    perms.fill(a)
+    half_point = nperms / 2
+    perms[half_point:, 0] = b
+    # Fill the rest of the table by sampling the pervious column every 2 items
+    for j in range(1, ndim):
+        half_col = perms[::2, j-1]
+        perms[:half_point, j] = half_col
+        perms[half_point:, j] = half_col
+
+    return perms
+
+
 def permutation_t_test(X, n_permutations=10000, tail=0):
     """One sample/paired sample permutation test based on a t-statistic.
 
@@ -27,8 +63,10 @@ def permutation_t_test(X, n_permutations=10000, tail=0):
         Data of size number of samples (aka number of observations) times
         number of tests (aka number of variables)
 
-    n_permutations : int
-        Number of permutations
+    n_permutations : int or 'exact'
+        Number of permutations. If n_permutations is 'exact' all possible
+        permutations are tested (2**n_samples).
+        If n_permutations >= 2**n_samples then the 'exact test is performed'
 
     tail : -1 or 0 or 1 (default = 0)
         If tail is 1, the alternative hypothesis is that the
@@ -51,12 +89,22 @@ def permutation_t_test(X, n_permutations=10000, tail=0):
     """
     n_samples, n_tests = X.shape
 
+    do_exact = False
+    if n_permutations is 'exact' or (n_permutations >= 2**n_samples - 1):
+        do_exact = True
+        n_permutations = 2**n_samples - 1
+
     X2 = np.mean(X**2, axis=0) # precompute moments
     mu0 = np.mean(X, axis=0)
     dof_scaling = sqrt(n_samples / (n_samples - 1.0))
     std0 = np.sqrt(X2 - mu0**2) * dof_scaling # get std with variance splitting
     T0 = np.mean(X, axis=0) / (std0 / sqrt(n_samples))
-    perms = np.sign(np.random.randn(n_permutations, n_samples))
+
+    if do_exact:
+        perms = bin_perm_rep(n_samples, a=1, b=-1)[1:,:]
+    else:
+        perms = np.sign(np.random.randn(n_permutations, n_samples))
+
     mus = np.dot(perms, X) / float(n_samples)
     stds = np.sqrt(X2[None,:] - mus**2) * dof_scaling # std with splitting
     max_abs = np.max(np.abs(mus) / (stds / sqrt(n_samples)), axis=1) # t-max
@@ -75,23 +123,43 @@ def permutation_t_test(X, n_permutations=10000, tail=0):
 
 permutation_t_test.__test__ = False # for nosetests
 
+
 if __name__ == '__main__':
     # 1 sample t-test
     n_samples, n_tests = 30, 5
+    n_permutations = 50000
+    # n_permutations = 'exact'
     X = np.random.randn(n_samples, n_tests)
-    X[:,:2] += 1
-    p_values, T0, H0 = permutation_t_test(X, n_permutations=999, tail=1)
+    X[:,:2] += 0.6
+    p_values, T0, H0 = permutation_t_test(X, n_permutations, tail=1)
     is_significant = p_values < 0.05
+    print 80*"-"
+    print "-------- 1-sample t-test :"
+    print "T stats : ", T0
+    print "p_values : ", p_values
+    print "Is significant : ", is_significant
+
+    print 80*"-"
+    print "-------- Comparison analytic vs permutation :"
+    p_values, T0, H0 = permutation_t_test(X, n_permutations, tail=1)
+    print "--- permutation_t_test :"
     print "T stats : ", T0
     print "p_values : ", p_values
     print "Is significant : ", is_significant
 
+    from scipy import stats
+    T0, p_values = stats.ttest_1samp(X[:,0], 0)
+    print "--- scipy.stats.ttest_1samp :"
+    print "T stats : ", T0
+    print "p_values : ", p_values
+
     # 2 samples t-test
-    n_samples, n_tests = 30, 5
     X1 = np.random.randn(n_samples, n_tests)
     X2 = np.random.randn(n_samples, n_tests)
     X1[:,:2] += 2
-    p_values, T0, H0 = permutation_t_test(X1 - X2, n_permutations=999)
+    p_values, T0, H0 = permutation_t_test(X1 - X2, n_permutations)
+    print 80*"-"
+    print "-------- 2-samples t-test :"
     print "T stats : ", T0
     print "p_values : ", p_values
     print "Is significant : ", is_significant
diff --git a/mne/stats/tests/test_permutations.py b/mne/stats/tests/test_permutations.py
index 8cbfeaf..640bbc6 100644
--- a/mne/stats/tests/test_permutations.py
+++ b/mne/stats/tests/test_permutations.py
@@ -5,10 +5,12 @@ from scipy import stats
 import mne
 from ..permutations import permutation_t_test
 
+
 def test_permutation_t_test():
     """Test T-test based on permutations
     """
     # 1 sample t-test
+    np.random.seed(10)
     n_samples, n_tests = 30, 5
     X = np.random.randn(n_samples, n_tests)
     X[:,:2] += 1
@@ -25,5 +27,8 @@ def test_permutation_t_test():
     is_significant = p_values < 0.05
     assert_array_equal(is_significant, [False, False, False, False, False])
 
-    assert_almost_equal(T0[0], stats.ttest_1samp(X[:,0], 0)[0], 8)
-
+    X = np.random.randn(18, 1)
+    p_values, T0, H0 = permutation_t_test(X[:, [0]], n_permutations='exact')
+    T0_scipy, p_values_scipy = stats.ttest_1samp(X[:, 0], 0)
+    assert_almost_equal(T0[0], T0_scipy, 8)
+    assert_almost_equal(p_values[0], p_values_scipy, 2)

-- 
Alioth's /usr/local/bin/git-commit-notice on /srv/git.debian.org/git/debian-med/python-mne.git