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.