[med-svn] [python-mne] 115/376: new example of TF cluster
Yaroslav Halchenko
debian at onerussian.com
Fri Nov 27 17:22:18 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 740f72ba1131fdc9c30a73034cb3b58d98df8c3e
Author: Alexandre Gramfort <alexandre.gramfort at inria.fr>
Date: Mon Mar 7 11:38:06 2011 -0500
new example of TF cluster
---
examples/stats/plot_cluster_stats_evoked.py | 21 ++--
.../stats/plot_cluster_stats_time_frequency.py | 117 +++++++++++++++++++++
mne/tfr.py | 5 +-
3 files changed, 131 insertions(+), 12 deletions(-)
diff --git a/examples/stats/plot_cluster_stats_evoked.py b/examples/stats/plot_cluster_stats_evoked.py
index 53dbe64..ff8e390 100644
--- a/examples/stats/plot_cluster_stats_evoked.py
+++ b/examples/stats/plot_cluster_stats_evoked.py
@@ -15,8 +15,6 @@ with cluster level permutation test.
print __doc__
-import numpy as np
-
import mne
from mne import fiff
from mne.stats import permutation_cluster_test
@@ -42,14 +40,17 @@ include = [channel]
# Read epochs for the channel of interest
picks = fiff.pick_types(raw['info'], meg=False, include=include)
event_id = 1
-data1, times, channel_names = mne.read_epochs(raw, events, event_id,
+epochs1 = mne.Epochs(raw, events, event_id,
tmin, tmax, picks=picks, baseline=(None, 0))
-condition1 = np.squeeze(np.array([d['epoch'] for d in data1])) # as 3D matrix
+condition1 = epochs1.get_data() # as 3D matrix
event_id = 2
-data2, times, channel_names = mne.read_epochs(raw, events, event_id,
+epochs2 = mne.Epochs(raw, events, event_id,
tmin, tmax, picks=picks, baseline=(None, 0))
-condition2 = np.squeeze(np.array([d['epoch'] for d in data2])) # as 3D matrix
+condition2 = epochs2.get_data() # as 3D matrix
+
+condition1 = condition1[:,0,:] # take only one channel to get a 2D array
+condition2 = condition2[:,0,:] # take only one channel to get a 2D array
###############################################################################
# Compute statistic
@@ -60,6 +61,7 @@ T_obs, clusters, cluster_p_values, H0 = \
###############################################################################
# Plot
+times = epochs1.times
import pylab as pl
pl.close('all')
pl.subplot(211)
@@ -69,11 +71,12 @@ pl.plot(times, condition1.mean(axis=0) - condition2.mean(axis=0),
pl.ylabel("MEG (T / m)")
pl.legend()
pl.subplot(212)
-for i_c, (start, stop) in enumerate(clusters):
+for i_c, c in enumerate(clusters):
+ c = c[0]
if cluster_p_values[i_c] <= 0.05:
- h = pl.axvspan(times[start], times[stop-1], color='r', alpha=0.3)
+ h = pl.axvspan(times[c.start], times[c.stop-1], color='r', alpha=0.3)
else:
- pl.axvspan(times[start], times[stop-1], color=(0.3, 0.3, 0.3),
+ pl.axvspan(times[c.start], times[c.stop-1], color=(0.3, 0.3, 0.3),
alpha=0.3)
hf = pl.plot(times, T_obs, 'g')
pl.legend((h, ), ('cluster p-value < 0.05', ))
diff --git a/examples/stats/plot_cluster_stats_time_frequency.py b/examples/stats/plot_cluster_stats_time_frequency.py
new file mode 100644
index 0000000..a5793b0
--- /dev/null
+++ b/examples/stats/plot_cluster_stats_time_frequency.py
@@ -0,0 +1,117 @@
+"""
+======================================================
+Non-parametric cluster statistic on single trial power
+======================================================
+
+This script shows how to estimate significant clusters
+in time-frequency power estimates. It uses a non-parametric
+statistical procedure based on permutations and cluster
+level statistics.
+
+"""
+# Authors: Alexandre Gramfort <gramfort at nmr.mgh.harvard.edu>
+#
+# License: BSD (3-clause)
+
+print __doc__
+
+import numpy as np
+
+import mne
+from mne import fiff
+from mne.tfr import single_trial_power
+from mne.stats import permutation_cluster_test
+from mne.datasets import sample
+
+###############################################################################
+# Set parameters
+data_path = sample.data_path('..')
+raw_fname = data_path + '/MEG/sample/sample_audvis_raw.fif'
+event_fname = data_path + '/MEG/sample/sample_audvis_raw-eve.fif'
+event_id = 1
+tmin = -0.2
+tmax = 0.5
+
+# Setup for reading the raw data
+raw = fiff.setup_read_raw(raw_fname)
+events = mne.read_events(event_fname)
+
+include = []
+exclude = raw['info']['bads'] + ['MEG 2443', 'EEG 053'] # bads + 2 more
+
+# picks MEG gradiometers
+picks = fiff.pick_types(raw['info'], meg='grad', eeg=False,
+ stim=False, include=include, exclude=exclude)
+
+picks = [picks[97]]
+ch_name = raw['info']['ch_names'][picks[0]]
+
+# Load condition 1
+event_id = 1
+epochs_condition_1 = mne.Epochs(raw, events, event_id,
+ tmin, tmax, picks=picks, baseline=(None, 0))
+data_condition_1 = epochs_condition_1.get_data() # as 3D matrix
+data_condition_1 *= 1e13 # change unit to fT / cm
+
+# Load condition 1
+event_id = 2
+epochs_condition_2 = mne.Epochs(raw, events, event_id,
+ tmin, tmax, picks=picks, baseline=(None, 0))
+data_condition_2 = epochs_condition_2.get_data() # as 3D matrix
+data_condition_2 *= 1e13 # change unit to fT / cm
+
+# Time vector
+times = 1e3 * epochs_condition_1.times # change unit to ms
+
+frequencies = np.arange(7, 30, 3) # define frequencies of interest
+Fs = raw['info']['sfreq'] # sampling in Hz
+epochs_coefs_1 = single_trial_power(data_condition_1, Fs=Fs,
+ frequencies=frequencies,
+ n_cycles=2, use_fft=False)
+
+epochs_coefs_2 = single_trial_power(data_condition_2, Fs=Fs,
+ frequencies=frequencies,
+ n_cycles=2, use_fft=False)
+
+epochs_coefs_1 = epochs_coefs_1[:,0,:,:] # only 1 channel to get a 3D matrix
+epochs_coefs_2 = epochs_coefs_2[:,0,:,:] # only 1 channel to get a 3D matrix
+ ###############################################################################
+# Compute statistic
+threshold = 6.0
+T_obs, clusters, cluster_p_values, H0 = \
+ permutation_cluster_test([epochs_coefs_1, epochs_coefs_2],
+ n_permutations=100, threshold=threshold, tail=0)
+
+###############################################################################
+# View time-frequency plots
+import pylab as pl
+pl.clf()
+pl.subplots_adjust(0.12, 0.08, 0.96, 0.94, 0.2, 0.43)
+pl.subplot(2, 1, 1)
+evoked_contrast = np.mean(data_condition_1, 0) - np.mean(data_condition_2, 0)
+pl.plot(times, evoked_contrast.T)
+pl.title('Contrast of evoked response (%s)' % ch_name)
+pl.xlabel('time (ms)')
+pl.ylabel('Magnetic Field (fT/cm)')
+pl.xlim(times[0], times[-1])
+pl.ylim(-100, 200)
+
+pl.subplot(2, 1, 2)
+
+# Create new stats image with only significant clusters
+T_obs_plot = np.nan * np.ones_like(T_obs)
+for c, p_val in zip(clusters, cluster_p_values):
+ if p_val <= 0.05:
+ T_obs_plot[c] = T_obs[c]
+
+pl.imshow(T_obs, cmap=pl.cm.gray, extent=[times[0], times[-1],
+ frequencies[0], frequencies[-1]],
+ aspect='auto', origin='lower')
+pl.imshow(T_obs_plot, cmap=pl.cm.jet, extent=[times[0], times[-1],
+ frequencies[0], frequencies[-1]],
+ aspect='auto', origin='lower')
+
+pl.xlabel('time (ms)')
+pl.ylabel('Frequency (Hz)')
+pl.title('Induced power (%s)' % ch_name)
+pl.show()
diff --git a/mne/tfr.py b/mne/tfr.py
index 2ac0b55..986f45b 100644
--- a/mne/tfr.py
+++ b/mne/tfr.py
@@ -189,8 +189,7 @@ def _time_frequency(X, Ws, use_fft):
return psd, plf
-def single_trial_power(epochs, Fs, frequencies, use_fft=True, n_cycles=25,
- n_jobs=1):
+def single_trial_power(epochs, Fs, frequencies, use_fft=True, n_cycles=7):
"""Compute time-frequency power on single epochs
"""
n_frequencies = len(frequencies)
@@ -210,7 +209,7 @@ def single_trial_power(epochs, Fs, frequencies, use_fft=True, n_cycles=25,
for k, e in enumerate(epochs):
mode = 'same'
- power[k] = np.abs(_cwt(e, Ws, mode))**2
+ power[k] = np.abs(list(_cwt(e, Ws, mode)))**2
return power
--
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