Spike train statistics

Statistical measures of spike trains (e.g., Fano factor) and functions to estimate firing rates.

elephant.statistics.complexity_pdf(spiketrains, binsize)

Complexity Distribution [1] of a list of neo.SpikeTrain objects.

Probability density computed from the complexity histogram which is the histogram of the entries of the population histogram of clipped (binary) spike trains computed with a bin width of binsize. It provides for each complexity (== number of active neurons per bin) the number of occurrences. The normalization of that histogram to 1 is the probability density.

Parameters:

spiketrains : List of neo.SpikeTrain objects

Spiketrains with a common time axis (same `t_start` and `t_stop`)

binsize : quantities.Quantity

Width of the histogram’s time bins.

Returns:

time_hist : neo.AnalogSignal

A neo.AnalogSignal object containing the histogram values.

`AnalogSignal[j]` is the histogram computed between .

See also

elephant.conversion.BinnedSpikeTrain

References

[1]Gruen, S., Abeles, M., & Diesmann, M. (2008). Impact of higher-order correlations on coincidence distributions of massively parallel data. In Dynamic Brain-from Neural Spikes to Behaviors (pp. 96-114). Springer Berlin Heidelberg.

elephant.statistics.cost_function(x, N, w, dt)

The cost function Cn(w) = sum_{i,j} int k(x - x_i) k(x - x_j) dx - 2 sum_{i~=j} k(x_i - x_j)

elephant.statistics.cv2(v)

Calculate the measure of CV2 for a sequence of time intervals between events.

Given a vector v containing a sequence of intervals, the CV2 is defined as:

.math $$ CV2 := frac{1}{N}sum_{i=1}^{N-1}

frac{2|isi_{i+1}-isi_i|}

{|isi_{i+1}+isi_i|} $$

The CV2 is typically computed as a substitute for the classical coefficient of variation (CV) for sequences of events which include some (relatively slow) rate fluctuation. As with the CV, CV2=1 for a sequence of intervals generated by a Poisson process.

Parameters:

v : quantity array, numpy array or list

Vector of consecutive time intervals

Returns:

cv2 : float

The CV2 of the inter-spike interval of the input sequence.

Raises:

AttributeError :

If an empty list is specified, or if the sequence has less than two entries, an AttributeError will be raised.

AttributeError :

Only vector inputs are supported. If a matrix is passed to the function an AttributeError will be raised.

References

..[1] Holt, G. R., Softky, W. R., Koch, C., & Douglas, R. J. (1996). Comparison of discharge variability in vitro and in vivo in cat visual cortex neurons. Journal of neurophysiology, 75(5), 1806-1814.

elephant.statistics.fanofactor(spiketrains)

Evaluates the empirical Fano factor F of the spike counts of a list of neo.core.SpikeTrain objects.

Given the vector v containing the observed spike counts (one per spike train) in the time window [t0, t1], F is defined as:

F := var(v)/mean(v).

The Fano factor is typically computed for spike trains representing the activity of the same neuron over different trials. The higher F, the larger the cross-trial non-stationarity. In theory for a time-stationary Poisson process, F=1.

Parameters:

spiketrains : list of neo.SpikeTrain objects, quantity arrays, numpy arrays or lists

Spike trains for which to compute the Fano factor of spike counts.

Returns:

fano : float or nan

The Fano factor of the spike counts of the input spike trains. If an empty list is specified, or if all spike trains are empty, F:=nan.

elephant.statistics.fftkernel(x, w)

Function `fftkernel’ applies the Gauss kernel smoother to an input signal using FFT algorithm.

Input argument x: Sample signal vector. w: Kernel bandwidth (the standard deviation) in unit of the sampling resolution of x.

Output argument y: Smoothed signal.

MAY 5/23, 2012 Author Hideaki Shimazaki RIKEN Brain Science Insitute http://2000.jukuin.keio.ac.jp/shimazaki

Ported to Python: Subhasis Ray, NCBS. Tue Jun 10 10:42:38 IST 2014

elephant.statistics.instantaneous_rate(spiketrain, sampling_period, kernel='auto', cutoff=5.0, t_start=None, t_stop=None, trim=False)

Estimates instantaneous firing rate by kernel convolution.

Parameters:

spiketrain : ‘neo.SpikeTrain’

Neo object that contains spike times, the unit of the time stamps and t_start and t_stop of the spike train.

sampling_period : Time Quantity

Time stamp resolution of the spike times. The same resolution will be assumed for the kernel

kernel : string ‘auto’ or callable object of Kernel from module

‘kernels.py’. Currently implemented kernel forms are rectangular, triangular, epanechnikovlike, gaussian, laplacian, exponential, and alpha function. Example: kernel = kernels.RectangularKernel(sigma=10*ms, invert=False) The kernel is used for convolution with the spike train and its standard deviation determines the time resolution of the instantaneous rate estimation. Default: ‘auto’. In this case, the optimized kernel width for the rate estimation is calculated according to [1] and with this width a gaussian kernel is constructed. Automatized calculation of the kernel width is not available for other than gaussian kernel shapes.

cutoff : float

This factor determines the cutoff of the probability distribution of the kernel, i.e., the considered width of the kernel in terms of multiples of the standard deviation sigma. Default: 5.0

t_start : Time Quantity (optional)

Start time of the interval used to compute the firing rate. If None assumed equal to spiketrain.t_start Default: None

t_stop : Time Quantity (optional)

End time of the interval used to compute the firing rate (included). If None assumed equal to spiketrain.t_stop Default: None

trim : bool

if False, the output of the Fast Fourier Transformation being a longer vector than the input vector by the size of the kernel is reduced back to the original size of the considered time interval of the spiketrain using the median of the kernel. if True, only the region of the convolved signal is returned, where there is complete overlap between kernel and spike train. This is achieved by reducing the length of the output of the Fast Fourier Transformation by a total of two times the size of the kernel, and t_start and t_stop are adjusted. Default: False

Returns:

rate : neo.AnalogSignal

Contains the rate estimation in unit hertz (Hz). Has a property ‘rate.times’ which contains the time axis of the rate estimate. The unit of this property is the same as the resolution that is given via the argument ‘sampling_period’ to the function.

Raises:

TypeError:

If spiketrain is not an instance of SpikeTrain of Neo. If sampling_period is not a time quantity. If kernel is neither instance of Kernel or string ‘auto’. If cutoff is neither float nor int. If t_start and t_stop are neither None nor a time quantity. If trim is not bool.

ValueError:

If sampling_period is smaller than zero.

References

..[1] H. Shimazaki, S. Shinomoto, J Comput Neurosci (2010) 29:171–182.

elephant.statistics.isi(spiketrain, axis=-1)

Return an array containing the inter-spike intervals of the SpikeTrain.

Accepts a Neo SpikeTrain, a Quantity array, or a plain NumPy array. If either a SpikeTrain or Quantity array is provided, the return value will be a quantities array, otherwise a plain NumPy array. The units of the quantities array will be the same as spiketrain.

Parameters:

spiketrain : Neo SpikeTrain or Quantity array or NumPy ndarray

The spike times.

axis : int, optional

The axis along which the difference is taken. Default is the last axis.

Returns:

NumPy array or quantities array.

elephant.statistics.lv(v)

Calculate the measure of local variation LV for a sequence of time intervals between events.

Given a vector v containing a sequence of intervals, the LV is defined as:

.math $$ LV := frac{1}{N}sum_{i=1}^{N-1}

frac{3(isi_i-isi_{i+1})^2}

{(isi_i+isi_{i+1})^2} $$

The LV is typically computed as a substitute for the classical coefficient of variation for sequences of events which include some (relatively slow) rate fluctuation. As with the CV, LV=1 for a sequence of intervals generated by a Poisson process.

Parameters:

v : quantity array, numpy array or list

Vector of consecutive time intervals

Returns:

lvar : float

The LV of the inter-spike interval of the input sequence.

Raises:

AttributeError :

If an empty list is specified, or if the sequence has less than two entries, an AttributeError will be raised.

ValueError :

Only vector inputs are supported. If a matrix is passed to the function a ValueError will be raised.

References

..[1] Shinomoto, S., Shima, K., & Tanji, J. (2003). Differences in spiking patterns among cortical neurons. Neural Computation, 15, 2823–2842.

elephant.statistics.make_kernel(form, sigma, sampling_period, direction=1)

Creates kernel functions for convolution.

Constructs a numeric linear convolution kernel of basic shape to be used for data smoothing (linear low pass filtering) and firing rate estimation from single trial or trial-averaged spike trains.

Exponential and alpha kernels may also be used to represent postynaptic currents / potentials in a linear (current-based) model.

Parameters:

form : {‘BOX’, ‘TRI’, ‘GAU’, ‘EPA’, ‘EXP’, ‘ALP’}

Kernel form. Currently implemented forms are BOX (boxcar), TRI (triangle), GAU (gaussian), EPA (epanechnikov), EXP (exponential), ALP (alpha function). EXP and ALP are asymmetric kernel forms and assume optional parameter direction.

sigma : Quantity

Standard deviation of the distribution associated with kernel shape. This parameter defines the time resolution of the kernel estimate and makes different kernels comparable (cf. [1] for symmetric kernels). This is used here as an alternative definition to the cut-off frequency of the associated linear filter.

sampling_period : float

Temporal resolution of input and output.

direction : {-1, 1}

Asymmetric kernels have two possible directions. The values are -1 or 1, default is 1. The definition here is that for direction = 1 the kernel represents the impulse response function of the linear filter. Default value is 1.

Returns:

kernel : numpy.ndarray

Array of kernel. The length of this array is always an odd number to represent symmetric kernels such that the center bin coincides with the median of the numeric array, i.e for a triangle, the maximum will be at the center bin with equal number of bins to the right and to the left.

norm : float

For rate estimates. The kernel vector is normalized such that the sum of all entries equals unity sum(kernel)=1. When estimating rate functions from discrete spike data (0/1) the additional parameter norm allows for the normalization to rate in spikes per second.

For example: rate = norm * scipy.signal.lfilter(kernel, 1, spike_data)

m_idx : int

Index of the numerically determined median (center of gravity) of the kernel function.

See also

elephant.statistics.instantaneous_rate

References

[1]	Meier R, Egert U, Aertsen A, Nawrot MP, “FIND - a unified framework for neural data analysis”; Neural Netw. 2008 Oct; 21(8):1085-93.

[2]	Nawrot M, Aertsen A, Rotter S, “Single-trial estimation of neuronal firing rates - from single neuron spike trains to population activity”; J. Neurosci Meth 94: 81-92; 1999.

Examples

To obtain single trial rate function of trial one should use:

r = norm * scipy.signal.fftconvolve(sua, kernel)

To obtain trial-averaged spike train one should use:

r_avg = norm * scipy.signal.fftconvolve(sua, np.mean(X,1))

where X is an array of shape (l,n), n is the number of trials and l is the length of each trial.

elephant.statistics.mean_firing_rate(spiketrain, t_start=None, t_stop=None, axis=None)

Return the firing rate of the SpikeTrain.

Accepts a Neo SpikeTrain, a Quantity array, or a plain NumPy array. If either a SpikeTrain or Quantity array is provided, the return value will be a quantities array, otherwise a plain NumPy array. The units of the quantities array will be the inverse of the spiketrain.

The interval over which the firing rate is calculated can be optionally controlled with t_start and t_stop

Parameters:

spiketrain : Neo SpikeTrain or Quantity array or NumPy ndarray

The spike times.

t_start : float or Quantity scalar, optional

The start time to use for the interval. If not specified, retrieved from the``t_start` attribute of spiketrain. If that is not present, default to 0. Any value from spiketrain below this value is ignored.

t_stop : float or Quantity scalar, optional

The stop time to use for the time points. If not specified, retrieved from the t_stop attribute of spiketrain. If that is not present, default to the maximum value of spiketrain. Any value from spiketrain above this value is ignored.

axis : int, optional

The axis over which to do the calculation. Default is None, do the calculation over the flattened array.

Returns:

float, quantities scalar, NumPy array or quantities array.

Raises:

TypeError

If spiketrain is a NumPy array and t_start or t_stop is a quantity scalar.

Notes

If spiketrain is a Quantity or Neo SpikeTrain and t_start or t_stop are not, t_start and t_stop are assumed to have the same units as spiketrain.

elephant.statistics.nextpow2(x)

Return the smallest integral power of 2 that >= x

elephant.statistics.oldfct_instantaneous_rate(spiketrain, sampling_period, form, sigma='auto', t_start=None, t_stop=None, acausal=True, trim=False)

Estimate instantaneous firing rate by kernel convolution.

Parameters:

spiketrain: ‘neo.SpikeTrain’

Neo object that contains spike times, the unit of the time stamps and t_start and t_stop of the spike train.

sampling_period : Quantity

time stamp resolution of the spike times. the same resolution will be assumed for the kernel

form : {‘BOX’, ‘TRI’, ‘GAU’, ‘EPA’, ‘EXP’, ‘ALP’}

Kernel form. Currently implemented forms are BOX (boxcar), TRI (triangle), GAU (gaussian), EPA (epanechnikov), EXP (exponential), ALP (alpha function). EXP and ALP are asymmetric kernel forms and assume optional parameter direction.

sigma : string or Quantity

Standard deviation of the distribution associated with kernel shape. This parameter defines the time resolution of the kernel estimate and makes different kernels comparable (cf. [1] for symmetric kernels). This is used here as an alternative definition to the cut-off frequency of the associated linear filter. Default value is ‘auto’. In this case, the optimized kernel width for the rate estimation is calculated according to [1]. Note that the automatized calculation of the kernel width ONLY works for gaussian kernel shapes!

t_start : Quantity (Optional)

start time of the interval used to compute the firing rate, if None assumed equal to spiketrain.t_start Default:None

t_stop : Qunatity

End time of the interval used to compute the firing rate (included). If none assumed equal to spiketrain.t_stop Default:None

acausal : bool

if True, acausal filtering is used, i.e., the gravity center of the filter function is aligned with the spike to convolve Default:None

m_idx : int

index of the value in the kernel function vector that corresponds to its gravity center. this parameter is not mandatory for symmetrical kernels but it is required when asymmetrical kernels are to be aligned at their gravity center with the event times if None is assumed to be the median value of the kernel support Default : None

trim : bool

if True, only the ‘valid’ region of the convolved signal are returned, i.e., the points where there isn’t complete overlap between kernel and spike train are discarded NOTE: if True and an asymmetrical kernel is provided the output will not be aligned with [t_start, t_stop]

Returns:

rate : neo.AnalogSignal

Contains the rate estimation in unit hertz (Hz). Has a property ‘rate.times’ which contains the time axis of the rate estimate. The unit of this property is the same as the resolution that is given as an argument to the function.

Raises:

TypeError:

If argument value for the parameter sigma is not a quantity object or string ‘auto’.

See also

elephant.statistics.make_kernel

References

..[1] H. Shimazaki, S. Shinomoto, J Comput Neurosci (2010) 29:171–182.

elephant.statistics.sskernel(spiketimes, tin=None, w=None, bootstrap=False)

Calculates optimal fixed kernel bandwidth.

spiketimes: sequence of spike times (sorted to be ascending).

tin: (optional) time points at which the kernel bandwidth is to be estimated.

w: (optional) vector of kernel bandwidths. If specified, optimal bandwidth is selected from this.

bootstrap (optional): whether to calculate the 95% confidence interval. (default False)

Returns

A dictionary containing the following key value pairs:

‘y’: estimated density, ‘t’: points at which estimation was computed, ‘optw’: optimal kernel bandwidth, ‘w’: kernel bandwidths examined, ‘C’: cost functions of w, ‘confb95’: (lower bootstrap confidence level, upper bootstrap confidence level), ‘yb’: bootstrap samples.

Ref: Shimazaki, Hideaki, and Shigeru Shinomoto. 2010. Kernel Bandwidth Optimization in Spike Rate Estimation. Journal of Computational Neuroscience 29 (1-2): 171-82. doi:10.1007/s10827-009-0180-4.

elephant.statistics.time_histogram(spiketrains, binsize, t_start=None, t_stop=None, output='counts', binary=False)

Time Histogram of a list of neo.SpikeTrain objects.

Parameters:

spiketrains : List of neo.SpikeTrain objects

Spiketrains with a common time axis (same t_start and t_stop)

binsize : quantities.Quantity

Width of the histogram’s time bins.

t_start, t_stop : Quantity (optional)

Start and stop time of the histogram. Only events in the input spiketrains falling between t_start and t_stop (both included) are considered in the histogram. If t_start and/or t_stop are not specified, the maximum t_start of all :attr:spiketrains is used as t_start, and the minimum t_stop is used as t_stop. Default: t_start = t_stop = None

output : str (optional)

Normalization of the histogram. Can be one of: * counts’: spike counts at each bin (as integer numbers) * `mean: mean spike counts per spike train * rate: mean spike rate per spike train. Like ‘mean’, but the

counts are additionally normalized by the bin width.

binary : bool (optional)

If True, indicates whether all spiketrain objects should first binned to a binary representation (using the BinnedSpikeTrain class in the conversion module) and the calculation of the histogram is based on this representation. Note that the output is not binary, but a histogram of the converted, binary representation. Default: False

Returns:

time_hist : neo.AnalogSignal

A neo.AnalogSignal object containing the histogram values. AnalogSignal[j] is the histogram computed between t_start + j * binsize and t_start + (j + 1) * binsize.

See also

elephant.conversion.BinnedSpikeTrain