...
 
Commits (3)
......@@ -8,6 +8,7 @@ which has its own page:
:maxdepth: 2
api/simulation
api/preprocessing
api/fitting
api/diagnostics
api/decomposition
......
......@@ -19,6 +19,8 @@ Model fitting helper functions
In addition to basic model fitting routines, some helper
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.pca_reduced_fit
Preprocessing Functions
=======================
.. autoclass:: sails.utils.PCA
:members:
......@@ -574,14 +574,43 @@ __all__.append('sliding_window_fit')
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
back to original data dimensions and the pca object used for reduction
"""
Helper for computing an MVAR on dimensionality reduced data (using PCA).
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]
from .utils import PCA
pc = PCA(X.T, ndim)
Xred = pc.scores.T
Xred = Xred[:, :, None]
......@@ -591,13 +620,17 @@ def pca_reduced_fit(X, delay_vect, ndim, linear_model=VieiraMorfLinearModel):
# Project back to data dimensions
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)):
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]
# Store projected model
proj_model = AbstractLinearModel()
proj_model = linear_model()
proj_model.parameters = Afull
proj_model.resid_cov = resid_cov_full
proj_model.delay_vect = delay_vect
......
......@@ -6,6 +6,8 @@ import math
import numpy as np
from scipy import signal, stats
from .anam import AbstractAnam, register_class
__all__ = []
......@@ -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',
segment_len=100, gesd_args=None, ret_mode='bad_inds'):
"""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
find artefacts in.
reject_dim : {'channels','samples','trials'}
......@@ -334,20 +338,19 @@ def detect_artefacts(X, reject_dim='channels', reject_metric='std',
return out
class PCA:
"""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
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.
class PCA(AbstractAnam):
"""
Class for handling PCA-based reduction before fitting a model
"""
Can also be the basis for some of the methods
in: https://doi.org/10.1016/j.jneumeth.2016.12.016 if interesting...
hdf5_outputs = ['npcs', 'data_mean', '_U', '_s', '_VT',
'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
if data is not None:
self.data_mean = data.mean(axis=0)
self._pca_svd(data, npcs)
......@@ -414,3 +417,7 @@ class PCA:
Compute a projection of a set of scores back to original data space
"""
return self.data_mean + np.dot(scores, self.components)
__all__.append('PCA')
register_class(PCA)