[med-svn] [python-mne] 228/353: ENH : add basic support for temporal whitening with AR model
Yaroslav Halchenko
debian at onerussian.com
Fri Nov 27 17:25:03 UTC 2015
This is an automated email from the git hooks/post-receive script.
yoh pushed a commit to tag 0.4
in repository python-mne.
commit e915ad04a348c03ec1a7379454e142398d60e336
Author: Alexandre Gramfort <alexandre.gramfort at inria.fr>
Date: Mon Jun 25 15:10:03 2012 +0200
ENH : add basic support for temporal whitening with AR model
---
doc/source/whats_new.rst | 2 +
examples/time_frequency/plot_temporal_whitening.py | 63 ++++++++++++
mne/time_frequency/__init__.py | 1 +
mne/time_frequency/ar.py | 111 +++++++++++++++++++++
mne/time_frequency/tests/test_ar.py | 42 ++++++++
5 files changed, 219 insertions(+)
diff --git a/doc/source/whats_new.rst b/doc/source/whats_new.rst
index 6dc2680..1a6d0aa 100644
--- a/doc/source/whats_new.rst
+++ b/doc/source/whats_new.rst
@@ -7,6 +7,8 @@ Current
Changelog
~~~~~~~~~
+ - Fit AR models to raw data for temporal whitening by `Alex Gramfort`_.
+
- speedup + reduce memory of mne.morph_data by `Alex Gramfort`_.
- Backporting scipy.signa.firwin2 so filtering works with old scipy by `Alex Gramfort`_.
diff --git a/examples/time_frequency/plot_temporal_whitening.py b/examples/time_frequency/plot_temporal_whitening.py
new file mode 100644
index 0000000..eeebbfd
--- /dev/null
+++ b/examples/time_frequency/plot_temporal_whitening.py
@@ -0,0 +1,63 @@
+"""
+================================
+Temporal whitening with AR model
+================================
+
+This script shows how to fit an AR model to data and use it
+to temporally whiten the signals.
+
+"""
+# Authors: Alexandre Gramfort <gramfort at nmr.mgh.harvard.edu>
+#
+# License: BSD (3-clause)
+
+print __doc__
+
+import numpy as np
+from scipy import signal
+import pylab as pl
+
+import mne
+from mne.time_frequency import ar_raw
+from mne.datasets import sample
+data_path = sample.data_path('..')
+
+raw_fname = data_path + '/MEG/sample/sample_audvis_raw.fif'
+proj_fname = data_path + '/MEG/sample/sample_audvis_ecg_proj.fif'
+
+raw = mne.fiff.Raw(raw_fname)
+proj = mne.read_proj(proj_fname)
+raw.info['projs'] += proj
+raw.info['bads'] = ['MEG 2443', 'EEG 053'] # mark bad channels
+
+# Set up pick list: Gradiometers - bad channels
+picks = mne.fiff.pick_types(raw.info, meg='grad', exclude=raw.info['bads'])
+
+order = 5 # define model order
+picks = picks[:5]
+
+# Estimate AR models on raw data
+coefs = ar_raw(raw, order=order, picks=picks, tmin=60, tmax=180)
+mean_coefs = np.mean(coefs, axis=0) # mean model accross channels
+
+filt = np.r_[1, -mean_coefs] # filter coefficient
+d, times = raw[0, 1e4:2e4] # look at one channel from now on
+d = d.ravel() # make flat vector
+innovation = signal.convolve(d, filt, 'valid')
+d_ = signal.lfilter([1], filt, innovation) # regenerate the signal
+d_ = np.r_[d_[0] * np.ones(order), d_] # dummy samples to keep signal length
+
+###############################################################################
+# Plot the different time series and PSDs
+pl.close('all')
+pl.figure()
+pl.plot(d[:100], label='signal')
+pl.plot(d_[:100], label='regenerated signal')
+pl.legend()
+
+pl.figure()
+pl.psd(d, Fs=raw.info['sfreq'], NFFT=2048)
+pl.psd(innovation, Fs=raw.info['sfreq'], NFFT=2048)
+pl.psd(d_, Fs=raw.info['sfreq'], NFFT=2048, linestyle='--')
+pl.legend(('Signal', 'Innovation', 'Regenerated signal'))
+pl.show()
diff --git a/mne/time_frequency/__init__.py b/mne/time_frequency/__init__.py
index 34395c5..4826123 100644
--- a/mne/time_frequency/__init__.py
+++ b/mne/time_frequency/__init__.py
@@ -3,3 +3,4 @@
from .tfr import induced_power, single_trial_power
from .psd import compute_raw_psd
+from .ar import yule_walker, ar_raw
diff --git a/mne/time_frequency/ar.py b/mne/time_frequency/ar.py
new file mode 100644
index 0000000..3daafc1
--- /dev/null
+++ b/mne/time_frequency/ar.py
@@ -0,0 +1,111 @@
+# Authors: Alexandre Gramfort <gramfort at nmr.mgh.harvard.edu>
+# The statsmodels folks for AR yule_walker
+#
+# License: BSD (3-clause)
+
+import numpy as np
+from scipy.linalg import toeplitz
+
+
+# XXX : Back ported from statsmodels
+
+def yule_walker(X, order=1, method="unbiased", df=None, inv=False, demean=True):
+ """
+ Estimate AR(p) parameters from a sequence X using Yule-Walker equation.
+
+ Unbiased or maximum-likelihood estimator (mle)
+
+ See, for example:
+
+ http://en.wikipedia.org/wiki/Autoregressive_moving_average_model
+
+ Parameters
+ ----------
+ X : array-like
+ 1d array
+ order : integer, optional
+ The order of the autoregressive process. Default is 1.
+ method : string, optional
+ Method can be "unbiased" or "mle" and this determines denominator in
+ estimate of autocorrelation function (ACF) at lag k. If "mle", the
+ denominator is n=X.shape[0], if "unbiased" the denominator is n-k.
+ The default is unbiased.
+ df : integer, optional
+ Specifies the degrees of freedom. If `df` is supplied, then it is assumed
+ the X has `df` degrees of freedom rather than `n`. Default is None.
+ inv : bool
+ If inv is True the inverse of R is also returned. Default is False.
+ demean : bool
+ True, the mean is subtracted from `X` before estimation.
+
+ Returns
+ -------
+ rho
+ The autoregressive coefficients
+ sigma
+ TODO
+
+ """
+#TODO: define R better, look back at notes and technical notes on YW.
+#First link here is useful
+#http://www-stat.wharton.upenn.edu/~steele/Courses/956/ResourceDetails/YuleWalkerAndMore.htm
+ method = str(method).lower()
+ if method not in ["unbiased", "mle"]:
+ raise ValueError("ACF estimation method must be 'unbiased' or 'MLE'")
+ X = np.array(X)
+ if demean:
+ X -= X.mean() # automatically demean's X
+ n = df or X.shape[0]
+
+ if method == "unbiased": # this is df_resid ie., n - p
+ denom = lambda k: n - k
+ else:
+ denom = lambda k: n
+ if X.ndim > 1 and X.shape[1] != 1:
+ raise ValueError("expecting a vector to estimate AR parameters")
+ r = np.zeros(order+1, np.float64)
+ r[0] = (X**2).sum() / denom(0)
+ for k in range(1,order+1):
+ r[k] = (X[0:-k]*X[k:]).sum() / denom(k)
+ R = toeplitz(r[:-1])
+
+ rho = np.linalg.solve(R, r[1:])
+ sigmasq = r[0] - (r[1:]*rho).sum()
+ if inv == True:
+ return rho, np.sqrt(sigmasq), np.linalg.inv(R)
+ else:
+ return rho, np.sqrt(sigmasq)
+
+
+def ar_raw(raw, order, picks, tmin=None, tmax=None):
+ """Fit AR model on raw data
+
+ Fit AR models for each channels and returns the models
+ coefficients for each of them.
+
+ Parameters
+ ----------
+ raw : Raw instance
+ The raw data
+ order : int
+ The AR model order
+ picks : array of int
+ The channels indices to include
+ tmin : float
+ The beginning of time interval in seconds.
+ tmax : float
+ The end of time interval in seconds.
+
+ Returns
+ -------
+ coefs : array
+ Sets of coefficients for each channel
+ """
+ start, stop = raw.time_to_index(tmin, tmax)
+ data, times = raw[picks, start:(stop + 1)]
+
+ coefs = np.empty((len(data), order))
+ for k, d in enumerate(data):
+ this_coefs, _ = yule_walker(d, order=order)
+ coefs[k, :] = this_coefs
+ return coefs
diff --git a/mne/time_frequency/tests/test_ar.py b/mne/time_frequency/tests/test_ar.py
new file mode 100644
index 0000000..dfd6819
--- /dev/null
+++ b/mne/time_frequency/tests/test_ar.py
@@ -0,0 +1,42 @@
+import os.path as op
+import numpy as np
+from numpy.testing import assert_array_almost_equal
+from nose.tools import assert_true
+import nose
+
+from ... import fiff
+from .. import yule_walker, ar_raw
+
+raw_fname = op.join(op.dirname(__file__), '..', '..', 'fiff', 'tests', 'data',
+ 'test_raw.fif')
+
+
+def test_yule_walker():
+ """Test Yule-Walker against statsmodels
+ """
+ try:
+ from statsmodels.regression.linear_model import yule_walker as sm_yw
+ d = np.random.randn(100)
+ sm_rho, sm_sigma = sm_yw(d, order=2)
+ rho, sigma = yule_walker(d, order=2)
+ assert_array_almost_equal(sm_sigma, sigma)
+ assert_array_almost_equal(sm_rho, rho)
+ except ImportError:
+ raise nose.SkipTest("XFailed Test")
+
+
+def test_ar_raw():
+ raw = fiff.Raw(raw_fname)
+
+ # picks MEG gradiometers
+ picks = fiff.pick_types(raw.info, meg='grad')
+
+ picks = picks[:2]
+
+ tmin, tmax = 0, 10 # use the first s of data
+ order = 2
+ coefs = ar_raw(raw, picks=picks, order=order, tmin=tmin, tmax=tmax)
+ mean_coefs = np.mean(coefs, axis=0)
+
+ assert_true(coefs.shape == (len(picks), order))
+ assert_true(0.9 < mean_coefs[0] < 1.1)
--
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