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

elephant.statistics.
complexity_pdf
(spiketrains, binsize)[source]¶ 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
References
[1]Gruen, S., Abeles, M., & Diesmann, M. (2008). Impact of higherorder correlations on coincidence distributions of massively parallel data. In Dynamic Brainfrom Neural Spikes to Behaviors (pp. 96114). Springer Berlin Heidelberg.

elephant.statistics.
cost_function
(x, N, w, dt)[source]¶ 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)[source]¶ 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}^{N1}
 frac{2isi_{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 interspike 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), 18061814.

elephant.statistics.
fanofactor
(spiketrains)[source]¶ 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 crosstrial nonstationarity. In theory for a timestationary 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)[source]¶ 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)[source]¶ 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 ofKernel
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)[source]¶ Return an array containing the interspike 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)[source]¶ 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}^{N1}
 frac{3(isi_iisi_{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 interspike 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)[source]¶ 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 trialaveraged spike trains.
Exponential and alpha kernels may also be used to represent postynaptic currents / potentials in a linear (currentbased) 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 cutoff 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.
References
[1] Meier R, Egert U, Aertsen A, Nawrot MP, “FIND  a unified framework for neural data analysis”; Neural Netw. 2008 Oct; 21(8):108593. [2] Nawrot M, Aertsen A, Rotter S, “Singletrial estimation of neuronal firing rates  from single neuron spike trains to population activity”; J. Neurosci Meth 94: 8192; 1999. Examples
To obtain single trial rate function of trial one should use:
r = norm * scipy.signal.fftconvolve(sua, kernel)
To obtain trialaveraged 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)[source]¶ 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.
oldfct_instantaneous_rate
(spiketrain, sampling_period, form, sigma='auto', t_start=None, t_stop=None, acausal=True, trim=False)[source]¶ 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 cutoff 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
References
..[1] H. Shimazaki, S. Shinomoto, J Comput Neurosci (2010) 29:171–182.

elephant.statistics.
sskernel
(spiketimes, tin=None, w=None, bootstrap=False)[source]¶ 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 (12): 17182. doi:10.1007/s1082700901804.

elephant.statistics.
time_histogram
(spiketrains, binsize, t_start=None, t_stop=None, output='counts', binary=False)[source]¶ 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