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 73713a854ef7adb50f34637cd15709a67b0a8fb6 Author: Alexandre Gramfort <[email protected]> Date: Tue Mar 15 12:31:39 2011 -0400 bug fix in permutation_cluster_t_test with tail = -1 --- mne/stats/cluster_level.py | 68 +++++++++++++++++++++++------------ mne/stats/tests/test_cluster_level.py | 57 +++++++++++++++++------------ 2 files changed, 80 insertions(+), 45 deletions(-) diff --git a/mne/stats/cluster_level.py b/mne/stats/cluster_level.py index fb87998..098cde7 100644 --- a/mne/stats/cluster_level.py +++ b/mne/stats/cluster_level.py @@ -63,6 +63,29 @@ def _find_clusters(x, threshold, tail=0): return clusters, sums +def _pval_from_histogram(T, H0, tail): + """Get p-values from stats values given an H0 distribution + + For each stat compute a p-value as percentile of its statistics + within all statistics in surrogate data + """ + if not tail in [-1, 0, 1]: + raise ValueError('invalid tail parameter') + + pval = np.array([percentileofscore(H0, t) for t in T]) + + # from pct to fraction + if tail == -1: # up tail + pval = pval / 100.0 + elif tail == 1: # low tail + pval = (100.0 - pval) / 100.0 + elif tail == 0: # both tails + pval = 100.0 - pval + pval += np.array([percentileofscore(H0, -t) for t in T]) + + return pval + + def permutation_cluster_test(X, stat_fun=f_oneway, threshold=1.67, n_permutations=1000, tail=0): """Cluster-level statistical permutation test @@ -127,7 +150,6 @@ def permutation_cluster_test(X, stat_fun=f_oneway, threshold=1.67, # Step 2: If we have some clusters, repeat process on permuted data # ------------------------------------------------------------------- if len(clusters) > 0: - cluster_pv = np.ones(len(clusters), dtype=np.float) H0 = np.zeros(n_permutations) # histogram for i_s in range(n_permutations): np.random.shuffle(X_full) @@ -140,12 +162,7 @@ def permutation_cluster_test(X, stat_fun=f_oneway, threshold=1.67, else: H0[i_s] = 0 - # for each cluster in original data, calculate p-value as percentile - # of its cluster statistics within all cluster statistics in surrogate - # data - cluster_pv[:] = [percentileofscore(H0, cluster_stats[i_cl]) - for i_cl in range(len(clusters))] - cluster_pv[:] = (100.0 - cluster_pv[:]) / 100.0 # from pct to fraction + cluster_pv = _pval_from_histogram(cluster_stats, H0, tail) return T_obs, clusters, cluster_pv, H0 else: return T_obs, np.array([]), np.array([]), np.array([]) @@ -211,8 +228,7 @@ def permutation_cluster_t_test(X, threshold=1.67, n_permutations=1000, tail=0): # Step 2: If we have some clusters, repeat process on permuted data # ------------------------------------------------------------------- if len(clusters) > 0: - cluster_pv = np.ones(len(clusters), dtype=np.float) - H0 = np.zeros(n_permutations) # histogram + H0 = np.empty(n_permutations) # histogram for i_s in range(n_permutations): # new surrogate data with random sign flip signs = np.sign(0.5 - np.random.rand(n_samples, *shape_ones)) @@ -223,16 +239,13 @@ def permutation_cluster_t_test(X, threshold=1.67, n_permutations=1000, tail=0): _, perm_clusters_sums = _find_clusters(T_obs_surr, threshold, tail) if len(perm_clusters_sums) > 0: - H0[i_s] = np.max(perm_clusters_sums) + idx_max = np.argmax(np.abs(perm_clusters_sums)) + H0[i_s] = perm_clusters_sums[idx_max] # get max with sign info else: H0[i_s] = 0 - # for each cluster in original data, calculate p-value as percentile - # of its cluster statistics within all cluster statistics in surrogate - # data - cluster_pv[:] = [percentileofscore(H0, cluster_stats[i_cl]) - for i_cl in range(len(clusters))] - cluster_pv[:] = (100.0 - cluster_pv[:]) / 100.0 # from pct to fraction + cluster_pv = _pval_from_histogram(cluster_stats, H0, tail) + return T_obs, clusters, cluster_pv, H0 else: return T_obs, np.array([]), np.array([]), np.array([]) @@ -268,16 +281,25 @@ permutation_cluster_t_test.__test__ = False # condition1[..., -3:] = np.random.randn(*shape1) * noiselevel # condition2[..., -3:] = np.random.randn(*shape2) * noiselevel # -# fs, clusters, cluster_p_values, histogram = permutation_cluster_t_test( -# condition1, n_permutations=100) +# # X, threshold, tail = condition1, 1.67, 1 +# # X, threshold, tail = -condition1, -1.67, -1 +# # # X, threshold, tail = condition1, 1.67, 0 +# # fs, clusters, cluster_p_values, histogram = permutation_cluster_t_test( +# # condition1, n_permutations=500, tail=tail, +# # threshold=threshold) # -# # fs, clusters, cluster_p_values, histogram = permutation_cluster_test( -# # [condition1, condition2], n_permutations=1000) +# # import pylab as pl +# # pl.close('all') +# # pl.hist(histogram) +# # pl.show() # +# fs, clusters, cluster_p_values, histogram = permutation_cluster_test( +# [condition1, condition2], n_permutations=1000) +# # # Plotting for a better understanding # import pylab as pl # pl.close('all') -# +# # if condition1.ndim == 2: # pl.subplot(211) # pl.plot(condition1.mean(axis=0), label="Condition 1") @@ -299,7 +321,7 @@ permutation_cluster_t_test.__test__ = False # for c, p_val in zip(clusters, cluster_p_values): # if p_val <= 0.05: # fs_plot[c] = fs[c] -# +# # pl.imshow(fs.T, cmap=pl.cm.gray) # pl.imshow(fs_plot.T, cmap=pl.cm.jet) # # pl.imshow(fs.T, cmap=pl.cm.gray, alpha=0.6) @@ -307,5 +329,5 @@ permutation_cluster_t_test.__test__ = False # pl.xlabel('time') # pl.ylabel('Freq') # pl.colorbar() -# +# # pl.show() diff --git a/mne/stats/tests/test_cluster_level.py b/mne/stats/tests/test_cluster_level.py index 8e2f8e2..5c7d78a 100644 --- a/mne/stats/tests/test_cluster_level.py +++ b/mne/stats/tests/test_cluster_level.py @@ -1,27 +1,29 @@ import numpy as np -from numpy.testing import assert_equal +from numpy.testing import assert_equal, assert_array_equal -from ..cluster_level import permutation_cluster_test, permutation_cluster_t_test +from ..cluster_level import permutation_cluster_test, \ + permutation_cluster_t_test +noiselevel = 20 -def test_cluster_permutation_test(): - """Test cluster level permutations tests. - """ - noiselevel = 20 +normfactor = np.hanning(20).sum() + +condition1 = np.random.randn(40, 350) * noiselevel +for c in condition1: + c[:] = np.convolve(c, np.hanning(20), mode="same") / normfactor - normfactor = np.hanning(20).sum() +condition2 = np.random.randn(33, 350) * noiselevel +for c in condition2: + c[:] = np.convolve(c, np.hanning(20), mode="same") / normfactor - condition1 = np.random.randn(40, 350) * noiselevel - for c in condition1: - c[:] = np.convolve(c, np.hanning(20), mode="same") / normfactor +pseudoekp = 5 * np.hanning(150)[None,:] +condition1[:, 100:250] += pseudoekp +condition2[:, 100:250] -= pseudoekp - condition2 = np.random.randn(33, 350) * noiselevel - for c in condition2: - c[:] = np.convolve(c, np.hanning(20), mode="same") / normfactor - pseudoekp = 5 * np.hanning(150)[None,:] - condition1[:, 100:250] += pseudoekp - condition2[:, 100:250] -= pseudoekp +def test_cluster_permutation_test(): + """Test cluster level permutations tests. + """ T_obs, clusters, cluster_p_values, hist = permutation_cluster_test( [condition1, condition2], n_permutations=500, @@ -33,12 +35,23 @@ def test_cluster_permutation_test(): tail=0) assert_equal(np.sum(cluster_p_values < 0.05), 1) - T_obs, clusters, cluster_p_values, hist = permutation_cluster_test( - [condition1, condition2], n_permutations=500, - tail=-1) - assert_equal(np.sum(cluster_p_values < 0.05), 0) - condition1 = condition1[:,:,None] # to test 2D also +def test_cluster_permutation_t_test(): + """Test cluster level permutations T-test. + """ + + my_condition1 = condition1[:,:,None] # to test 2D also T_obs, clusters, cluster_p_values, hist = permutation_cluster_t_test( - condition1, n_permutations=500, tail=0) + my_condition1, n_permutations=500, tail=0) assert_equal(np.sum(cluster_p_values < 0.05), 1) + + T_obs_pos, _, cluster_p_values_pos, _ = permutation_cluster_t_test( + my_condition1, n_permutations=500, tail=1, + threshold=1.67) + + T_obs_neg, _, cluster_p_values_neg, _ = permutation_cluster_t_test( + -my_condition1, n_permutations=500, tail=-1, + threshold=-1.67) + assert_array_equal(T_obs_pos, -T_obs_neg) + assert_array_equal(cluster_p_values_pos < 0.05, + cluster_p_values_neg < 0.05) -- Alioth's /usr/local/bin/git-commit-notice on /srv/git.debian.org/git/debian-med/python-mne.git _______________________________________________ debian-med-commit mailing list [email protected] http://lists.alioth.debian.org/cgi-bin/mailman/listinfo/debian-med-commit
