{
  "nbformat_minor": 0, 
  "nbformat": 4, 
  "cells": [
    {
      "execution_count": null, 
      "cell_type": "code", 
      "source": [
        "%matplotlib inline"
      ], 
      "outputs": [], 
      "metadata": {
        "collapsed": false
      }
    }, 
    {
      "source": [
        "\n\n# Non-parametric 1 sample cluster statistic on single trial power\n\n\nThis script shows how to estimate significant clusters\nin time-frequency power estimates. It uses a non-parametric\nstatistical procedure based on permutations and cluster\nlevel statistics.\n\nThe procedure consists in:\n\n  - extracting epochs\n  - compute single trial power estimates\n  - baseline line correct the power estimates (power ratios)\n  - compute stats to see if ratio deviates from 1.\n\n\n"
      ], 
      "cell_type": "markdown", 
      "metadata": {}
    }, 
    {
      "execution_count": null, 
      "cell_type": "code", 
      "source": [
        "# Authors: Alexandre Gramfort <alexandre.gramfort@telecom-paristech.fr>\n#\n# License: BSD (3-clause)\n\nimport numpy as np\nimport matplotlib.pyplot as plt\n\nimport mne\nfrom mne.time_frequency import tfr_morlet\nfrom mne.stats import permutation_cluster_1samp_test\nfrom mne.datasets import sample\n\nprint(__doc__)"
      ], 
      "outputs": [], 
      "metadata": {
        "collapsed": false
      }
    }, 
    {
      "source": [
        "Set parameters\n--------------\n\n"
      ], 
      "cell_type": "markdown", 
      "metadata": {}
    }, 
    {
      "execution_count": null, 
      "cell_type": "code", 
      "source": [
        "data_path = sample.data_path()\nraw_fname = data_path + '/MEG/sample/sample_audvis_raw.fif'\ntmin, tmax, event_id = -0.3, 0.6, 1\n\n# Setup for reading the raw data\nraw = mne.io.read_raw_fif(raw_fname)\nevents = mne.find_events(raw, stim_channel='STI 014')\n\ninclude = []\nraw.info['bads'] += ['MEG 2443', 'EEG 053']  # bads + 2 more\n\n# picks MEG gradiometers\npicks = mne.pick_types(raw.info, meg='grad', eeg=False, eog=True,\n                       stim=False, include=include, exclude='bads')\n\n# Load condition 1\nevent_id = 1\nepochs = mne.Epochs(raw, events, event_id, tmin, tmax, picks=picks,\n                    baseline=(None, 0), preload=True,\n                    reject=dict(grad=4000e-13, eog=150e-6))\n\n# Take only one channel\nch_name = 'MEG 1332'\nepochs.pick_channels([ch_name])\n\nevoked = epochs.average()\n\n# Factor to down-sample the temporal dimension of the TFR computed by\n# tfr_morlet. Decimation occurs after frequency decomposition and can\n# be used to reduce memory usage (and possibly computational time of downstream\n# operations such as nonparametric statistics) if you don't need high\n# spectrotemporal resolution.\ndecim = 5\nfrequencies = np.arange(8, 40, 2)  # define frequencies of interest\nsfreq = raw.info['sfreq']  # sampling in Hz\ntfr_epochs = tfr_morlet(epochs, frequencies, n_cycles=4., decim=decim,\n                        average=False, return_itc=False, n_jobs=1)\n\n# Baseline power\ntfr_epochs.apply_baseline(mode='logratio', baseline=(-.100, 0))\n\n# Crop in time to keep only what is between 0 and 400 ms\nevoked.crop(0., 0.4)\ntfr_epochs.crop(0., 0.4)\n\nepochs_power = tfr_epochs.data[:, 0, :, :]  # take the 1 channel"
      ], 
      "outputs": [], 
      "metadata": {
        "collapsed": false
      }
    }, 
    {
      "source": [
        "Compute statistic\n-----------------\n\n"
      ], 
      "cell_type": "markdown", 
      "metadata": {}
    }, 
    {
      "execution_count": null, 
      "cell_type": "code", 
      "source": [
        "threshold = 2.5\nT_obs, clusters, cluster_p_values, H0 = \\\n    permutation_cluster_1samp_test(epochs_power, n_permutations=100,\n                                   threshold=threshold, tail=0)"
      ], 
      "outputs": [], 
      "metadata": {
        "collapsed": false
      }
    }, 
    {
      "source": [
        "View time-frequency plots\n-------------------------\n\n"
      ], 
      "cell_type": "markdown", 
      "metadata": {}
    }, 
    {
      "execution_count": null, 
      "cell_type": "code", 
      "source": [
        "evoked_data = evoked.data\ntimes = 1e3 * evoked.times\n\nplt.figure()\nplt.subplots_adjust(0.12, 0.08, 0.96, 0.94, 0.2, 0.43)\n\n# Create new stats image with only significant clusters\nT_obs_plot = np.nan * np.ones_like(T_obs)\nfor c, p_val in zip(clusters, cluster_p_values):\n    if p_val <= 0.05:\n        T_obs_plot[c] = T_obs[c]\n\nvmax = np.max(np.abs(T_obs))\nvmin = -vmax\nplt.subplot(2, 1, 1)\nplt.imshow(T_obs, cmap=plt.cm.gray,\n           extent=[times[0], times[-1], frequencies[0], frequencies[-1]],\n           aspect='auto', origin='lower', vmin=vmin, vmax=vmax)\nplt.imshow(T_obs_plot, cmap=plt.cm.RdBu_r,\n           extent=[times[0], times[-1], frequencies[0], frequencies[-1]],\n           aspect='auto', origin='lower', vmin=vmin, vmax=vmax)\nplt.colorbar()\nplt.xlabel('Time (ms)')\nplt.ylabel('Frequency (Hz)')\nplt.title('Induced power (%s)' % ch_name)\n\nax2 = plt.subplot(2, 1, 2)\nevoked.plot(axes=[ax2])\nplt.show()"
      ], 
      "outputs": [], 
      "metadata": {
        "collapsed": false
      }
    }
  ], 
  "metadata": {
    "kernelspec": {
      "display_name": "Python 2", 
      "name": "python2", 
      "language": "python"
    }, 
    "language_info": {
      "mimetype": "text/x-python", 
      "nbconvert_exporter": "python", 
      "name": "python", 
      "file_extension": ".py", 
      "version": "2.7.13", 
      "pygments_lexer": "ipython2", 
      "codemirror_mode": {
        "version": 2, 
        "name": "ipython"
      }
    }
  }
}