...
 
Commits (3)
...@@ -8,6 +8,7 @@ which has its own page: ...@@ -8,6 +8,7 @@ which has its own page:
:maxdepth: 2 :maxdepth: 2
api/simulation api/simulation
api/preprocessing
api/fitting api/fitting
api/diagnostics api/diagnostics
api/decomposition api/decomposition
......
...@@ -19,6 +19,8 @@ Model fitting helper functions ...@@ -19,6 +19,8 @@ Model fitting helper functions
In addition to basic model fitting routines, some helper In addition to basic model fitting routines, some helper
functions exist; for instance to help with fitting a number of functions exist; for instance to help with fitting a number of
models to a data series following a "sliding window" pattern. models to a data series following a "sliding window" pattern
or by applying PCA before the fit.
.. autofunction:: sails.modelfit.sliding_window_fit .. autofunction:: sails.modelfit.sliding_window_fit
.. autofunction:: sails.modelfit.pca_reduced_fit
Preprocessing Functions
=======================
.. autoclass:: sails.utils.PCA
:members:
...@@ -574,14 +574,43 @@ __all__.append('sliding_window_fit') ...@@ -574,14 +574,43 @@ __all__.append('sliding_window_fit')
def pca_reduced_fit(X, delay_vect, ndim, linear_model=VieiraMorfLinearModel): def pca_reduced_fit(X, delay_vect, ndim, linear_model=VieiraMorfLinearModel):
""" First pass at a helper for computing an MVAR on dimensionality reduced """
data. Returns a model fitted to reduced data, the reduced model projected Helper for computing an MVAR on dimensionality reduced data (using PCA).
back to original data dimensions and the pca object used for reduction
Returns a
model fitted to reduced data, the reduced model projected back to original
data dimensions and the pca object used for reduction
Parameters
----------
X : ndarray
The data to compute the reduced model fit on
delay_vect : ndarray
A vector of lags specifying which lags to fit
ndim : int
Number of components to reduce data to.
linear_model : class
Subclass of AbstractLinearModel to use for fit.
Defaults to VieiraMorfLinearModel.
(Default value = True)
Returns
-------
red_model : AbstractLinearModel
An instance of the linear model passed into the function containing MVAR
parameters for all windows for the reduced model.
proj_model : AbstractLinearModel
An instance of the linear model passed into the function containing MVAR
parameters for all windows for the projected model.
pc : sails.utils.PCA
PCA class with information about the PCA projection included
""" """
# PCA requires [samples x channels] whereas data is [channels x samples] # PCA requires [samples x channels] whereas data is [channels x samples]
from .utils import PCA from .utils import PCA
pc = PCA(X.T, ndim) pc = PCA(X.T, ndim)
Xred = pc.scores.T Xred = pc.scores.T
Xred = Xred[:, :, None] Xred = Xred[:, :, None]
...@@ -591,13 +620,17 @@ def pca_reduced_fit(X, delay_vect, ndim, linear_model=VieiraMorfLinearModel): ...@@ -591,13 +620,17 @@ def pca_reduced_fit(X, delay_vect, ndim, linear_model=VieiraMorfLinearModel):
# Project back to data dimensions # Project back to data dimensions
Afull = np.zeros((X.shape[0], X.shape[0], len(delay_vect))) Afull = np.zeros((X.shape[0], X.shape[0], len(delay_vect)))
Afull[:, :, 0] = -np.eye(X.shape[0]) # Don't project identity at start
# Don't project identity at start
Afull[:, :, 0] = -np.eye(X.shape[0])
for ii in range(1, len(delay_vect)): for ii in range(1, len(delay_vect)):
Afull[:, :, ii] = np.dot(pc.components.T, red_model.parameters[:, :, ii]).dot(pc.components) Afull[:, :, ii] = np.dot(pc.components.T, red_model.parameters[:, :, ii]).dot(pc.components)
resid_cov_full = np.dot(pc.components.T, red_model.resid_cov[:, :]).dot(pc.components)[:, :, None] resid_cov_full = np.dot(pc.components.T, red_model.resid_cov[:, :]).dot(pc.components)[:, :, None]
# Store projected model # Store projected model
proj_model = AbstractLinearModel() proj_model = linear_model()
proj_model.parameters = Afull proj_model.parameters = Afull
proj_model.resid_cov = resid_cov_full proj_model.resid_cov = resid_cov_full
proj_model.delay_vect = delay_vect proj_model.delay_vect = delay_vect
......
...@@ -6,6 +6,8 @@ import math ...@@ -6,6 +6,8 @@ import math
import numpy as np import numpy as np
from scipy import signal, stats from scipy import signal, stats
from .anam import AbstractAnam, register_class
__all__ = [] __all__ = []
...@@ -212,8 +214,10 @@ def gesd(x, alpha=0.05, p_out=.1, outlier_side=0): ...@@ -212,8 +214,10 @@ def gesd(x, alpha=0.05, p_out=.1, outlier_side=0):
def detect_artefacts(X, reject_dim='channels', reject_metric='std', def detect_artefacts(X, reject_dim='channels', reject_metric='std',
segment_len=100, gesd_args=None, ret_mode='bad_inds'): segment_len=100, gesd_args=None, ret_mode='bad_inds'):
"""Detect bad channels or segments in M/EEG data """Detect bad channels or segments in M/EEG data
Compute a Principal Components Analysis on a given dataset using the SVD
method. Parameters
----------
X : ndarray
Array of size [channels x samples] or [channels x samples x trials] to Array of size [channels x samples] or [channels x samples x trials] to
find artefacts in. find artefacts in.
reject_dim : {'channels','samples','trials'} reject_dim : {'channels','samples','trials'}
...@@ -334,22 +338,21 @@ def detect_artefacts(X, reject_dim='channels', reject_metric='std', ...@@ -334,22 +338,21 @@ def detect_artefacts(X, reject_dim='channels', reject_metric='std',
return out return out
class PCA: class PCA(AbstractAnam):
"""The class structure here could be seriously cleaned up and merged into """
main sails. Longer term it would be good to use this for dimensionality Class for handling PCA-based reduction before fitting a model
reduction of sensor data, and projecting A and associated metrics back into """
full space. Currently lots of redundancy and comments/links for my own
education whilst figuring this out.
Can also be the basis for some of the methods hdf5_outputs = ['npcs', 'data_mean', '_U', '_s', '_VT',
in: https://doi.org/10.1016/j.jneumeth.2016.12.016 if interesting... 'explained_variance', 'components', 'loadings', 'scores']
I'll update and merge into main sails def __init__(self, data=None, npcs=None):
""" AbstractAnam.__init__(self)
def __init__(self, data, npcs): # Need to cope with no-argument initialisation for anamnesis
self.data_mean = data.mean(axis=0) if data is not None:
self._pca_svd(data, npcs) self.data_mean = data.mean(axis=0)
self._pca_svd(data, npcs)
def _pca_svd(self, data, npcs=10): def _pca_svd(self, data, npcs=10):
""" """
...@@ -414,3 +417,7 @@ class PCA: ...@@ -414,3 +417,7 @@ class PCA:
Compute a projection of a set of scores back to original data space Compute a projection of a set of scores back to original data space
""" """
return self.data_mean + np.dot(scores, self.components) return self.data_mean + np.dot(scores, self.components)
__all__.append('PCA')
register_class(PCA)