Detection of non-stationary processes

This algorithm determines if a spike train spk can be considered as stationary process (constant firing rate) or not as stationary process (i.e. presence of one or more points at which the rate increases or decreases). In case of non-stationarity, the output is a list of detected Change Points (CPs). Essentially, a det of two-sided window of width h (_filter(t, h, spk)) slides over the spike train within the time [h, t_final-h]. This generates a _filter_process(time_step, h, spk) that assigns at each time t the difference between a spike lying in the right and left window. If at any time t this difference is large ‘enough’ is assumed the presence of a rate Change Point in a neighborhood of t. A threshold test_quantile for the maximum of the filter_process (max difference of spike count between the left and right window) is derived based on asymptotic considerations. The procedure is repeated for an arbitrary set of windows, with different size h.

Examples

The following applies multiple_filter_test to a spike trains.

>>> import quantities as pq
>>> import neo
>>> from elephant.change_point_detection import multiple_filter_test
...
>>> test_array = [1.1,1.2,1.4,   1.6,1.7,1.75,1.8,1.85,1.9,1.95]
>>> st = neo.SpikeTrain(test_array, units='s', t_stop = 2.1)
>>> window_size = [0.5]*pq.s
>>> t_fin = 2.1*pq.s
>>> alpha = 5.0
>>> num_surrogates = 10000
>>> change_points = multiple_filter_test(window_size, st, t_fin, alpha,
...                 num_surrogates, time_step = 0.5*pq.s)

References

Messer, M., Kirchner, M., Schiemann, J., Roeper, J., Neininger, R., & Schneider, G. (2014). A multiple filter test for the detection of rate changes in renewal processes with varying variance. The Annals of Applied Statistics, 8(4),2027-2067.

Original code

Adapted from the published R implementation: DOI: 10.1214/14-AOAS782SUPP;.r

elephant.change_point_detection.empirical_parameters(window_sizes, t_final, alpha, n_surrogates, time_step=None)[source]

This function generates the threshold and the null parameters. The`_filter_process_h` has been proved to converge (for t_fin, h–>infinity) to a continuous functional of a Brownaian motion (‘limit_process’). Using a MonteCarlo technique, maxima of these limit_processes are collected.

The threshold is defined as the alpha quantile of this set of maxima. Namely: test_quantile := alpha quantile of {max_(h in window_size)[

max_(t in [h, t_final-h])_limit_process_h(t)]}
Parameters:
window_sizeslist of quantity objects

set of windows’ size

t_finalquantity object

final time of the spike

alphafloat

alpha-quantile in range [0, 100]

n_surrogatesinteger

numbers of simulated limit processes

time_stepquantity object

resolution, time step at which the windows are slided

Returns:
test_quantilefloat

threshold for the maxima of the filter derivative processes, if any of these maxima is larger than this value, it is assumed the presence of a cp at the time corresponding to that maximum

test_paramnp.array 3 * num of window,

first row: list of h, second and third rows: empirical means and variances of the limit process correspodning to h. This will be used to normalize the filter_process in order to give to the every maximum the same impact on the global statistic.

elephant.change_point_detection.multiple_filter_test(window_sizes, spiketrain, t_final, alpha, n_surrogates, test_quantile=None, test_param=None, time_step=None)[source]

Detects change points.

This function returns the detected change points, that correspond to the maxima of the _filter_processes. These are the processes generated by sliding the windows of step time_step; at each step the difference between spike on the right and left window is calculated.

Parameters:
window_sizeslist of quantity objects

list that contains windows sizes

spiketrainneo.SpikeTrain, numpy array or list

spiketrain objects to analyze

t_finalquantity

final time of the spike train which is to be analysed

alphafloat

alpha-quantile in range [0, 100] for the set of maxima of the limit processes

n_surrogatesinteger

numbers of simulated limit processes

test_quantilefloat

threshold for the maxima of the filter derivative processes, if any of these maxima is larger than this value, it is assumed the presence of a cp at the time corresponding to that maximum

time_stepquantity

resolution, time step at which the windows are slided

test_paramnp.array of shape (3, num of window),

first row: list of h, second and third rows: empirical means and variances of the limit process correspodning to h. This will be used to normalize the filter_process in order to give to the every maximum the same impact on the global statistic.

Returns:
cpslist of list

one list for each window size h, containing the points detected with the corresponding filter_process. N.B.: only cps whose h-neighborhood does not include previously detected cps (with smaller window h) are added to the list.