Module reference¶
Model fitting and Diagnostics¶
Model fitting classes¶

class
sails.modelfit.
VieiraMorfLinearModel
[source]¶ TODO: Description

classmethod
fit_model
(data, delay_vect)[source]¶ Estimates the multichannel autoregressive spectral estimators using a VieiraMorf algorithm. Equations are referenced to Marple, appendix 15.B.
This is the multitrial versions of the algorithm, using the AMVAR method outlined in Ding 2000.
Parameters:  data (numpy.ndarray) – array of shape [nsignals, nsamples, ntrials]
 delay_vect (numpy.ndarray) – Vector containing evenly spaced delays to be assessed in the model. Delays are represented in samples. Must start with zero.
 assess_order (bool) – Optional True/False value indicating whether additional information about each iteration of the fit should be returned. Defaults to False
Returns: A populated object containing the fitted forward and backwards coefficients and several other useful variables and methods.

classmethod
Model fitting helper functions¶
Model diagnostics¶

class
sails.diags.
DelayDiagnostics
[source]¶ Class which computes the mutual information as a function of lag from zero lag to the first zero crossing of the autocorrelation function.

classmethod
delay_search
(data, maxdelay, step, sample_rate, constant_window=True)[source]¶ Compute MI as a function of lag from zero lag to the first zero crossing of the autocorrelation function
If NafOptions().verbose is True, progress will be reported to stdout.
Parameters:  data (numpy.ndarray) – array of [signals, samples, trials]
 maxdelay (int) – maximum delay to consider
 step (int) – step to increment the delay by
 sample_rate (float) – sample rate of data
 constant_window (bool) – Flag indicating that the same number of datapoints should be included at each delay. Default is True

classmethod
Model Decomposition¶

class
sails.modal.
MvarModalDecomposition
[source]¶ Object which calculates the modal decomposition of a fitted model derived from the AbstractLinearModel class.

compute_modal_parameters
()[source]¶ Compute the order 1 or 2 autoregressive coefficients associated with each pole or conjugate pair.

compute_residues
()[source]¶ Compute the residue matrix from the eigenvectors associated with each pole.

excitation
(resid_cov)[source]¶ http://climatedynamics.org/wpcontent/uploads/2015/08/arfit.pdf

classmethod
initialise
(model, sample_rate=None, normalise=False)[source]¶ Compute the modal poleresidue decomposition of the A matrix from a fitted MVAR model.
Currently only implemented for a single realisation (A.ndim == 3)
Parameters:  model – Fitted model object. Must be derived from AbstractLinearModel.
 sample_rate – sample rate of system in Hz
 normalise – Whether to adjust the phase of the eigenvectors
Returns: Modal decomposition object

modal_timecourse
(data)[source]¶ Compute the modal timecourse from a dataset and a fitted model. This is the tranformation of a set of [nsignals x time] data observations into modal coordinates of shape [nmodes x time].
Parameters: data – TODO

modal_transfer_function
(sample_rate, modes=None)[source]¶ Compute a transfer function matrix based on the excitation period for each mode
Parameters:  sample_rate – sample rate of system in Hz
 modes – List of mode indices to evaluate over (optional)

mode_indices
¶ Returns a list of tuples of mode indices where each tuple contains the indices of either a polepair or an unpaired pole
Returns: list of tuples containing indices into the modes

per_mode_transfer_function
(sample_rate, freq_vect)[source]¶ Extracts the transfer function for each mode by summing the transfer function across polepairs where necessary
Parameters:  sample_rate – sample rate of system in Hz
 freq_vect – Vector of frequencies at which to evaluate function

period_hz
¶ This is an old and deprecated way of accessing peak_frequency

pole_plot
(ax=None, normscale=50, plottype='unitcircle')[source]¶ Plot the poles of the current system.
Parameters:  ax – Optional axis on which to plot
 normscale – Scaling factor for points on plot
 plottyle – String defining which plot to make, can be ‘unitcircle’ or ‘eigenspectrum’
Returns: Reference to axes on which plot was drawn

residue_norm
¶ TODO

transfer_function
(sample_rate, freq_vect, modes=None, sum_modes=True)[source]¶ Compute the transfer function in poleresidue form, splitting the system into modes with separate transfer functions. The full system transfer function is then a linear sum of each modal transfer function.
Parameters:  sample_rate – sample rate of system in Hz
 freq_vect – Vector of frequencies at which to evaluate function
 modes – List of mode indices to evaluate over (optional)
 sum_modes – Boolean indicating whether to sum modes or return all

Model metric estimatation¶

class
sails.mvar_metrics.
ModalMvarMetrics
[source]¶ 
classmethod
initialise
(model, sample_rate, freq_vect, nmodes=None, sum_modes=True)[source]¶ Compute some basic values from the poleresidue decomposition of the A matrix
Currently only implemented for a single realisation (A.ndim == 3)
Parameters:  model – TODO
 sample_rate – TODO
 freq_vect – TODO
 nmodes – TODO

classmethod

sails.mvar_metrics.
modal_transfer_function
(evals, evecl, evecr, nchannels, sample_rate=None, freq_vect=None)[source]¶ Comute the transfer function in poleresidue form, splitting the system into modes with separate transfer functions. The full system transfer function is then a linear sum of each modal transfer function.
Plotting¶

sails.plotting.
root_plot
(rts, ax=None, figargs={}, plotargs={})[source]¶ Plot a set of roots (complex numbers).
Parameters:  rts – Roots to plot
 ax – Optional Axes on which to place plot.
Returns: Axes object on which plot was drawn

sails.plotting.
plot_vector
(metric, x_vect, y_vect=None, x_label=None, y_label=None, title=None, labels=None, line_labels=None, F=None, cmap=<matplotlib.colors.LinearSegmentedColormap object>, triangle=None, diag=False, thresh=None, two_tail_thresh=False, font_size=10, use_latex=False)[source]¶ Function for plotting frequency domain connectivity at a single timepoint.
Params metric: matrix containing connectivity values [nsignals x signals x frequencies x participants] in which the first dimension refers to source nodes and the second dimension refers to target nodes
Params x_vect: vector of frequencies to label the x axis
Params y_vect: y_vect: vector containing the values for the yaxis
Parameters:  x_label (string [optional]) – label for the x axis
 y_label (string [optional]) – label for the y axis
 title (string [optional]) – title for the figure
 labels (list) – list of node labels for columns and vectors
 line_labels (list) – list of labels for each separate line (participant dimension in metric)
 F (figurehandle [optional]) – handle of existing figure to plot within
 triangle (string [optional]) – string to indicate whether only the ‘upper’ or ‘lower’ triangle of the matrix should be plotted
 diag (bool [optional]) – flag to indicate whether the diagonal elements should be plotted
 thresh (ndarray [optional]) – matrix containing thesholds to be plotted alongside connectivity values [nsignals x nsignals x frequencies]
 two_tailed_thresh (bool [optional]) – flag to indicate whether both signs (+/) of the threshold should be plotted
 font_size (int [optional]) – override the default font size

sails.plotting.
plot_matrix
(metric, x_vect, y_vect, x_label=None, y_label=None, z_vect=None, title=None, labels=None, F=None, cmap=<matplotlib.colors.LinearSegmentedColormap object>, font_size=8, use_latex=False, diag=True)[source]¶ Function for plotting frequency domain connectivity over many time points
Parameters:  metric (ndarray) – matrix containing connectivity values [nsignals x signals x frequencies x participants] in which the first dimension refers to source nodes and the second dimension refers to target nodes
 x_vect (1d array) – vector of frequencies to label the x axis
 y_vect (1d array [optional]) – vector containing the values for the yaxis
 z_vect (1d array [optional]) – vector containing values for the colour scale
 x_label (string [optional]) – label for the x axis
 y_label (string [optional]) – label for the y axis
 title (string [optional]) – title for the figure
 labels (list) – list of node labels for columns and vectors
 F (figurehandle [optional]) – handle of existing figure to plot within
 cmap (matplotlib colormap [optional]) – matplotlib.cm.<colormapname> to use for colourscale (redundant for plot vector??)
 font_size (int [optional]) – override the default font size
Simulation¶

class
sails.simulate.
AbstractSigGen
(extraargs=None)[source]¶ Abstract base class for signal generation

generate_signal
(meg_reader)[source]¶ Generate a custom signal for embedding in a dataset
Parameters: meg_reader (naf.meg.abstractdata.AbstractMEGReader) – MEG data for which to generate signal Return type: numpy.ndarray Returns: numpy array of shape (samples). Note that samples must be equal to the value of meg_reader.num_samples


class
sails.simulate.
AbstractLagSigGen
(extraargs=None)[source]¶ Generic class for making lagged signals

generate_basesignal
(nsamples)[source]¶ Function defining the starting signal which will be lagged throughout the network

generate_signal
(sample_rate=None, num_samples=None, noise_dB=None)[source]¶ Method for generating a set of time courses based on either the structure of a meg_reader object or passed parameters and the self.generate_parameters method defined in this class
Returns a timeseries of dimension [samples x signals]

get_lags
()[source]¶ Method returning the lagging of the signal through the system. Should return a [nsignals x nsignals x nlags] matrix in which the first dimension refers to the source node and the second dimension the target node. Each nonzero value indicates that the signal in the first dimension will be added into the signal in the second dimension at the specified weight.


class
sails.simulate.
AbstractMVARSigGen
(extraargs=None)[source]¶ Generic class for making MVAR signals

check_stability
()[source]¶ Checks the stablility of the parameter matrix by the magnitude of the largest eigenvalue of the first lag in the parameter matrix. Lutkepohl 2005 ch 1 for more detail

generate_parameters
()[source]¶ Method to generate the parameters of the MVAR system. Should return a [nsignals x nsignals x nlags] parameter matrix in which the first dimension refers to the source node and the second dimension the target node

generate_signal
(sample_rate=None, num_samples=None, noise_dB=None, num_realisations=1)[source]¶ Method for generating a set of time courses based on either the structure of a meg_reader object or passed parameters and the self.generate_parameters method defined in this class
Returns a timeseries of dimension [samples x signals]


class
sails.simulate.
Baccala2001_fig1
(extraargs=None)[source]¶

class
sails.simulate.
Baccala2001_fig2
(extraargs=None)[source]¶

class
sails.simulate.
Baccala2001_fig3
(extraargs=None)[source]¶

class
sails.simulate.
Baccala2001_fig4
(extraargs=None)[source]¶

class
sails.simulate.
Baccala2001_fig5
(extraargs=None)[source]¶

class
sails.simulate.
Fasoula2013_eqn26
(extraargs=None)[source]¶

class
sails.simulate.
Schelter2006_fig1
(extraargs=None)[source]¶

class
sails.simulate.
PascualMarqui2014_fig3
(extraargs=None)[source]¶

class
sails.simulate.
Korzeniewska2003Lag
(extraargs=None)[source]¶ Generate realisations of the timelagged network defined in Korzeniewska 2003 figure 1.
https://doi.org/10.1016/S01650270(03)000529

generate_basesignal
(nsamples)[source]¶ Generate oscillatory signal to be propogated through the netowrk. This creates a single resonance at onehalf of the Nyqvist frequency via direct pole placement.
Note: the original basesignal in the paper is a realisation from an AR model. I couldn’t find the parameters for this mode, so use pole placement here. The spectral profile will be different to the paper, but the connectivity patterns are the same.

get_lags
()[source]¶ Method returning the lagging of the signal through the system. Should return a [nsignals x nsignals x nlags] matrix in which the first dimension refers to the source node and the second dimension the target node. Each nonzero value indicates that the signal in the first dimension will be added into the signal in the second dimension at the specified weight.


class
sails.simulate.
SailsTutorialExample
(extraargs=None)[source]¶ 
generate_basesignal
(f, r, sample_rate, num_samples)[source]¶ Generate a simple signal by pole placement

generate_signal
(f1, f2, sample_rate, num_samples, noise_dB=None, num_realisations=1, magnitude_noise=0.02)[source]¶ Method for generating a set of time courses based on either the structure of a meg_reader object or passed parameters and the self.generate_parameters method defined in this class
Returns a timeseries of dimension [samples x signals x realisations]

Support routines¶
Abstract Classes¶

class
sails.modelfit.
AbstractLinearModel
[source]¶ TODO: Description

companion
¶ Parameter matrix in companion form

compute_diagnostics
(data, compute_pc=False)[source]¶ Compute several diagnostic metrics from a fitted MVAR model and a dataset.
Return type: diags.ModelDiagnostics Returns: ModelDiagnostics object containing LL, AIC, BIC, stability ratio, durbinwatson, R squared and percent consistency measures

get_residuals
(data)[source]¶ Returns the prediction error from a fitted model. This is a wrapper function for get_residuals()
Parameters: data – TODO

nsignals
¶ Number of signals in fitted model

order
¶ Order of fitted model


class
sails.mvar_metrics.
AbstractMVARMetrics
[source]¶ TODO: Description

PSD
¶ Power Spectral Density matrix

S
¶ Spectral Matrix of shape [nsignals, nsignals, order, epochs]

coherency
¶ TODO: Document

d_directed_transfer_function
¶ TODO: Document

directed_transfer_function
¶ TODO: Document

ff_directed_transfer_function
¶ TODO: Document

geweke_granger_causality
¶ TODO: Document

imaginary_coherence
¶ TODO: Document

inv_S
¶ TODO: Document

isolated_effective_coherence
¶ TODO: Document

magnitude_squared_coherence
¶ TODO: Document

partial_coherence
¶ TODO: Document

partial_directed_coherence
¶ TODO: Document

phase_coherence
¶ TODO: Document

Statistical Functions¶

sails.mvar_metrics.
sdf_spectrum
(A, sigma, sample_rate, freq_vect)[source]¶ Estimate of the Spectral Density as found on wikipedia and Quirk & Liu 1983
This assumes that the spectral representation of A is invertable

sails.mvar_metrics.
psd_spectrum
(A, sigma, sample_rate, freq_vect)[source]¶ Estimate the PSD representation of a set of MVAR coefficients as stated in Penny 2000 Signal processing course array.ps 7.41
This assumes that the spectral representation of A is invertable
TODO: This doesn’t behave as expected for univariate?, sdf_spectrum recommended
Parameters:  A – TODO
 sigma – TODO
 sample_rate – TODO
 freq_vect – TODO
Returns: TODO

sails.mvar_metrics.
ar_spectrum
(A, sigma, sample_rate, freq_vect)[source]¶ Estimate the spectral representation of a set of MVAR coefficients as suggested by Baccala & Sameshima 2001, the equation without a number just below equation 13.
Parameters:  A – TODO
 sigma – TODO
 sample_rate – TODO
 freq_vect – TODO
Returns: TODO

sails.mvar_metrics.
transfer_function
(Af)[source]¶ Function for computing the transfer function of a system
Parameters: Af – Frequency domain version of parameters (can be calculated using ar_spectrum function) [nchannels x nchannels x nfreqs x nepochs] Returns: numpy array of dimensions [nsignals x nsignals x frequency x epochs]

sails.mvar_metrics.
spectral_matrix
(H, noise_cov)[source]¶ Function for computing the spectral matrix, the matrix of spectra and crossspectra.
Parameters:  H – The transfer matrix [nchannels x nchannels x nfreqs x nepochs]
 noise_cov – The noise covariance matrix of the system [nchannels x nchannels x nepochs]
Returns: numpy array of dimensions [nsignals x nsignals x frequency x epochs]

sails.mvar_metrics.
coherency
(S)[source]¶ Method for computing the Coherency. This is the complex form of coherence, from which metrics such as magnitude squared coherence can be derived
Params S: The spectral matrix: 3D or 4D [nchannels x nchannels x nfreqs x nepochs] Returns: numpy array of dimensions [nsignals x nsignals x frequency x epochs]

sails.mvar_metrics.
partial_coherence
(inv_S)[source]¶ Method for computing the Partial Coherence.
Params inv_S: The inverse spectral matrix [nchannels x nchannels x nfreqs x nepochs] Returns: numpy array of dimensions [nsignals x nsignals x frequency x epochs]

sails.mvar_metrics.
partial_directed_coherence
(Af)[source]¶ Function to estimate the partial directed coherence from a set of multivariate parameters.
Parameters: Af – Frequency domain version of parameters (can be calculated using ar_spectrum function) Returns: PDC : ndarray containing the PDC estimates from the parameters of dimensions [nchannels x nchannels x order]

sails.mvar_metrics.
isolated_effective_coherence
(Af, noise_cov)[source]¶ Function for estimating the Isolated Effective Coherence as defined in PascualMarqui et al 2014
Parameters:  Af – Frequency domain version of parameters (can be calculated using ar_spectrum function) [nchannels x nchannels x nfreqs x nepochs]
 noise_cov – The noise covariance matrix of the system [nchannels x nchannels x nepochs]
Returns: iec: numpy array of dimensions [nsignals x nsignals x frequency x epochs]

sails.mvar_metrics.
directed_transfer_function
(H)[source]¶ Method for computing the Directed Transfer Function as defined in Kaminski & Blinowska 1991
Parameters: H – The transfer matrix [nchannels x nchannels x nfreqs x nepochs] Returns: numpy array of dimensions [nsignals x nsignals x frequency x epochs]

sails.modal.
adjust_phase
(x)[source]¶ Method for normalising eigenvector matrices to meet the conditions in eqn 21 of Neumaier & Schneider 2001.
Parameters: x – Vector of complex numbers to adjust the phase of Returns: Vector of numbers with adjusted phase

sails.stats.
find_cov_multitrial
(a, b, ddof=1)[source]¶ Estimate the average covariance between two datasets across many trials
Parameters:  a (ndarray) – a multivariate time series of shape nsignals x nsamples x ntrials
 b (ndarray) – a multivariate time series of shape nsignals x nsamples x ntrials
 ddof (None orint) – if None, no normalisation will be performed. If an integer, covariance will be normalised by (nsamples  ddof)
Return type: ndarray
Returns: A covariance matrix between a and b of shape nsignals x nsignals

sails.stats.
find_cov
(a, b, ddof=1)[source]¶ Estimate the covariance between two datasets.
Parameters:  a (ndarray) – a multivariate time series of shape nsignals x nsamples
 b (ndarray) – a multivariate time series of shape nsignals x nsamples
 ddof (int) – if None, no normalisation will be performed. If an integer, covariance will be normalised by (nsamples  ddof)
Return type: ndarray
Returns: A covariance matrix between a and b of shape nsignals x nsignals

sails.stats.
mutual_information
(A, B, r=None, nbins=21, base=<ufunc 'log2'>)[source]¶ Calculate the mutual information between two signals.
Parameters:  A (ndarray) – First signal, of shape [channels x samples]
 B (ndarray) – Second signal, of shape [channels x samples]
 r (tuple (float, float)) – Range of the random variables used for bin calculations. If not set, takes the minimum and maximum values from the combination of A and B
 nbins (int) – Number of bins to compute the frequencies within
 base (function) – Base to use for MI calculation  defaults to np.log2
Return type: ndarray
Returns: the mutual information within each signal passed

sails.stats.
durbin_watson
(residuals, step=1)[source]¶ Test for autocorrelation in the residuals of a model.
The result is a value between 0 and 4. A value of 0 indicates a positive autocorrelaiton, 2 indicates no correlation and a value of 4 indicates a negative autocorrelation.
If many trials are found the statistic is estimated for each trial separately.
Parameters:  residuals (ndarray) – Residuals of shape [channels, samples, ntrials] If only two dimensions, will be treated as [channels, samples] If only one dimension, will be treated as [samples]
 step (int) – stepsize to use through residuals when calculating DW. Defaults to 1.
Return type: ndarray
Returns: durbinwatson index

sails.stats.
percent_consistency
(data, residuals)[source]¶ Estimate how well a model retains the correlational structure of a dataset.
This uses the percent consistency measure outlined in Ding et al 2000. All auto and cross correlations within a multivariate data set are computed and compared to the auto and cross correlations of the data as predicted by a model (in this case the data  model residuals, as these are easily obtained).
Parameters:  data (numpy.ndarray) – An array containing the observed data, of shape [nsignals x nsamples x ntrials]
 residuals – An array containing the remaining variance of data after model fitting, of shape [nsignals x nsamples x ntrials]
Rtype float: Returns: A value between 0 and 100 indicating the percentage of the autocorrelation in the data that is retained in the model.

sails.stats.
rsquared_from_residuals
(data, residuals, per_signal=False)[source]¶ Estimate the variance explained by a model from the original data and the model residuals
If data and residuals are 2D, they will be treated as [channels, nsamples]
If data and residuals are 1D, they will be treated as [nsamples]
Parameters:  data (ndarray) – Array containing the original data, of shape [channels, nsamples, ntrials]
 residuals (ndarray) – Array containing the model residuals, of shape [channels, nsamples, ntrials]
 per_signal (bool) – Boolean indicating whether the variance explained should be computed separately per channel
Return type: float
Returns: Value between 0 and 1 indicating the proportion of variance explained by the model.

sails.stats.
mahalanobis
(X, Y=None, sigma=None)[source]¶ Estimate the mahalanobis distance of each point in an array across the samples dimension.
Parameters:  X (ndarray) – Array containing data of points: shape [channels x samples]
 Y (ndarray or None) – Optional Array containing data of second points: shape [channels x samples]. If not provided, the origin will be used as the point for comparison.
 sigma (ndarray or None) – A precomputed covariance matrix to use in the calculation, of shape [channels x channels]. If not provided, sigma will be calculated from X
Return type: ndarray
Returns: The mahalanobis distance for each sample, of shape [samples]

sails.stats.
profile_likelihood
(X)[source]¶ Find the ‘elbow’ point in a curve.
This function uses the automatic dimensionality selection method from Zhu et al 2006 to estimate the elbow or inflection point in a given curve.
Method can be found in Mu Zhu, Ali Ghodsi. Computational Statistics & Data Analysis 51 (2006) 918  930
Parameters: X (ndarray) – 1d array containing the curve to be analysed, of shape [samples] Return type: ndarray Returns: The log likelihood that each point on the curve separates the distributions to the left and right, of shape [samples]

sails.stats.
stability_index
(A)[source]¶ Compute the stability index as defined in Lutkepohl 2006 pg 11.
This is a proxy for the stationarity assumption of MVAR models, if the magnitude of the largest principal component of the model parameters is less than 1, the model is stable and thus, stationary.
Parameters: a (numpy.ndarray) – The MVAR parameter matrix Return type: float Returns: The stability index

sails.stats.
stability_ratio
(A)[source]¶ Compute the stability ratio as a measure of stationarity
A stronger test for stationarity than the stability index. Computes the ratio between the largest eigenvalue of the parameters and all the others. If the largest parameter is larger than all other parmeters (indicated by a value < 1) we have a superstable system.
Parameters: a (numpy.ndarray) – The MVAR parameter matrix Rtype float: Returns: The stability ratio

sails.modelvalidation.
approx_model_selection_criteria
(A, EF, method='ss')[source]¶ Estimate Akaike’s Information Criterion and Bayes Information Criterion from a model fit using an approximation of the loglikelihood.
 Two approximations are available:
 ‘ss’: The sum square of the residuals. This approximation is outlined in Burnham & Anderson 2004 and Chatfield 2003. ‘det’: The determinant of the residual covariance matrix.
The sum squares method is used as default.
Parameters:  A – Array containing the parameters of a fitted model, of shape [channels x channels x maxorder]
 EF – Array containing the residuals of the fitted model, of shape [channels x samples x trials]
 method – string indicating which method should be used to approximate the loglikelihood (‘ss’ or ‘det’)
Returns: tuple of (aLL, AIC, BIC) aLL: Approximation of the loglikelihood (determined by the method defined above) AIC: Estimate of Akaike’s Information Criterion BIC: Estimate of the Bayesian Information Criterion
Tutorial helper functions¶

sails.tutorial_utils.
generate_pink_roots
(power, order=24)[source]¶ Generates the poles of a system with a pink spectrum.
Parameters:  power – Power of pole
 order – Order of system to generate
Returns: Roots of the system (length of order)