[med-svn] [python-mne] 68/376: ENH : adding exact permutation test with all permutations + some tests and cosmits
Yaroslav Halchenko
debian at onerussian.com
Fri Nov 27 17:22:09 UTC 2015
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
More information about the debian-med-commit
mailing list