Compare evokeds confidence intervals for grads

[ktavabi]

Still not clear to me if there is a bug in how confidence intervals for gradiometers in plot_compare_evokeds are computed? Using simulation example to illustrate...

import numpy as np
import matplotlib.pyplot as plt

import mne
from mne.datasets import sample
from mne.time_frequency import fit_iir_model_raw
from mne.viz import plot_sparse_source_estimates
from mne.simulation import simulate_sparse_stc, simulate_evoked
data_path = sample.data_path()

raw = mne.io.read_raw_fif(data_path + '/MEG/sample/sample_audvis_raw.fif')
proj = mne.read_proj(data_path + '/MEG/sample/sample_audvis_ecg-proj.fif')
raw.info['projs'] += proj

fwd = mne.read_forward_solution(data_path + 
                                '/MEG/sample/sample_audvis-meg-eeg-oct-6-fwd.fif')
fwd = mne.pick_types_forward(fwd, meg=True, eeg=True, exclude=[])
cov = mne.read_cov(data_path + 
                   '/MEG/sample/sample_audvis-cov.fif')
info = mne.io.read_info(data_path + 
                        '/MEG/sample/sample_audvis-no-filter-ave.fif')

label_names = ['Aud-lh', 'Aud-rh']
labels = [mne.read_label(data_path + '/MEG/sample/labels/%s.label' % ln)
          for ln in label_names]
times = np.arange(300, dtype=np.float) / raw.info['sfreq'] - 0.1
picks = mne.pick_types(raw.info, meg=True, exclude='bads')
iir_filter = fit_iir_model_raw(raw, order=5, picks=picks, tmin=60, tmax=180)[1]
nave = 100  # simulate average of 100 epoc
rng = np.random.RandomState(42)

def data_fun(times):
    """Function to generate random source time courses"""
    return (50e-9 * np.sin(30. * times) *
            np.exp(- (times - 0.15 + 0.05 * rng.randn(1)) ** 2 / 0.01))

evokeds = list()
for i in range(10):
    stc = simulate_sparse_stc(fwd['src'], n_dipoles=2, times=times,
                              random_state=42, labels=labels, data_fun=data_fun)
    evokeds.append(simulate_evoked(fwd, stc, info, cov, nave=nave, use_cps=True,
                                   iir_filter=iir_filter))
picks = [raw.ch_names.index(j) for j in ['MEG 1811',
                                          'MEG 1812',
                                          'MEG 1813']]
mne.viz.plot_compare_evokeds(evokeds=[evokeds], picks=picks,
                             ci=.95)

[agramfort]

yes CIs is broken when you pass gradiometers and they are combined for plotting but not to estimated the confidence intervals :(

I'll try to look over the next few days unless someone beats me to it...

[jona-sassenhagen]

@dengemann

[ktavabi]

I think I got the CI estimation for grads headed in the right direction with

...
if ch_type == 'grad':
    data = [(_merge_grad_data(evoked_.data[picks, :])).
                  mean(0) for evoked_ in evokeds[condition]]
    data = np.asarray(data)
else:
    data = np.asarray([evoked_.data[picks, :].mean(0)
                                    for evoked_ in evokeds[condition]])
ci_array[condition] = _ci_fun(data) * scaling
...

But I am still unable to plot them correctly, almost like CI's are scaled differently, and it's got me stumped.

Not sure if this is an inconvenience to anyone WRT #4526 convo, but I am curious to know if I am headed in the right direction.