elephant.statistics.instantaneous_rate

elephant.statistics.instantaneous_rate(spiketrains, sampling_period, kernel='auto', cutoff=5.0, t_start=None, t_stop=None, trim=False, center_kernel=True, border_correction=False, pool_trials=False, pool_spike_trains=False)[source]

Estimates instantaneous firing rate by kernel convolution.

Visualization of this function is covered in Viziphant: viziphant.statistics.plot_instantaneous_rates_colormesh().

Parameters:
spiketrainsneo.SpikeTrain, list of neo.SpikeTrain or elephant.trials.Trials # noqa

Input spike train(s) for which the instantaneous firing rate is calculated. If a list of spike trains is supplied, the parameter pool_spike_trains determines the behavior of the function. If a Trials object is supplied, the behavior is determined by the parameters pool_spike_trains (within a trial) and pool_trials (across trials).

sampling_periodpq.Quantity

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

kernel‘auto’ or Kernel, optional

The string ‘auto’ or callable object of class kernels.Kernel. The kernel is used for convolution with the spike train and its standard deviation determines the time resolution of the instantaneous rate estimation. Currently, implemented kernel forms are rectangular, triangular, epanechnikovlike, gaussian, laplacian, exponential, and alpha function. If ‘auto’, the optimized kernel width for the rate estimation is calculated according to (Shimazaki and Shinomoto, 2010) and a Gaussian kernel is constructed with this width. Automatized calculation of the kernel width is not available for other than Gaussian kernel shapes.

Note: The kernel width is not adaptive, i.e., it is calculated as global optimum across the data.

Default: ‘auto’

cutofffloat, optional

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_startpq.Quantity, optional

Start time of the interval used to compute the firing rate. If None, t_start is assumed equal to t_start attribute of spiketrain.

Default: None

t_stoppq.Quantity, optional

End time of the interval used to compute the firing rate. If None, t_stop is assumed equal to t_stop attribute of spiketrain.

Default: None

trimbool, optional

Accounts for the asymmetry of a kernel. If False, the output of the Fast Fourier Transformation being a longer vector than the input vector (ouput = input + kernel - 1) is reduced back to the original size of the considered time interval of the spiketrain using the median of the kernel. False (no trimming) is equivalent to ‘same’ convolution mode for symmetrical kernels. 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. True (trimming) is equivalent to ‘valid’ convolution mode for symmetrical kernels.

Default: False

center_kernelbool, optional

If set to True, the kernel will be translated such that its median is centered on the spike, thus putting equal weight before and after the spike. If False, no adjustment is performed such that the spike sits at the origin of the kernel.

Default: True

border_correctionbool, optional

Apply a border correction to prevent underestimating the firing rates at the borders of the spike trains, i.e., close to t_start and t_stop. The correction is done by estimating the mass of the kernel outside these spike train borders under the assumption that the rate does not change strongly. Only possible in the case of a Gaussian kernel.

Default: False

pool_trials: bool, optional

If true, calculate firing rates averaged over trials if spiketrains is of type elephant.trials.Trials Has no effect for single spike train or lists of spike trains.

Default: False

pool_spike_trains: bool, optional

If true, calculate firing rates averaged over spike trains. If the input is a Trials object, spike trains are pooled across spike trains within each trial, and pool_trials determines whether spike trains are additionally pooled across trials. Has no effect for a single spike train.

Default: False

Returns:
rateneo.AnalogSignal

2D matrix that contains the rate estimation in unit hertz (Hz) of shape (time, len(spiketrains)) or (time, 1) in case of a single input spiketrain. rate.times 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 neo.SpikeTrain.

  • If sampling_period is not a pq.Quantity.

  • If sampling_period is not larger than zero.

  • If kernel is neither instance of kernels.Kernel nor string ‘auto’.

  • If cutoff is neither float nor int.

  • If t_start and t_stop are neither None nor a pq.Quantity.

  • If trim is not bool.

ValueError

  • If sampling_period is smaller than zero.

  • If kernel is ‘auto’ and the function was unable to calculate optimal kernel width for instantaneous rate from input data.

Warns:
UserWarning

  • If cutoff is less than min_cutoff attribute of kernel, the width of the kernel is adjusted to a minimally allowed width.

Notes

  • The resulting instantaneous firing rate values smaller than 0, which may happen due to machine precision errors, are clipped to zero.

  • The instantaneous firing rate estimate is calculated based on half-open intervals [), except the last one e.g. if t_start = 0s, t_stop = 4s and sampling_period = 1s, the intervals are:

    [0, 1) [1, 2) [2, 3) [3, 4].

    This induces a sampling bias, which can lead to a time shift of the estimated rate, if the sampling_period is chosen large relative to the duration (t_stop - t_start). One possibility to counteract this is to choose a smaller sampling_period.

  • The last interval of the given duration (t_stop - t_start) is dropped if it is shorter than sampling_period, e.g. if t_start = 0s, t_stop = 4.5s and sampling_period = 1s, the intervals considered are:

    [0, 1) [1, 2) [2, 3) [3, 4],

    the last interval [4, 4.5] is excluded from all calculations.

Examples

Example 1. Automatic kernel estimation.

>>> import neo
>>> import quantities as pq
>>> from elephant import statistics
>>> spiketrain = neo.SpikeTrain([0.3, 4.5, 6.7, 9.3], t_stop=10, units='s')
>>> rate = statistics.instantaneous_rate(spiketrain,
...                                      sampling_period=10 * pq.ms,
...                                      kernel='auto')
>>> rate.annotations['kernel']
{'type': 'GaussianKernel', 'sigma': '7.273225922958104 s', 'invert': False}
>>> print(rate.sampling_rate)
0.1 1/ms

Example 2. Manually set kernel.

>>> from elephant import kernels
>>> spiketrain = neo.SpikeTrain([0], t_stop=1, units='s')
>>> kernel = kernels.GaussianKernel(sigma=300 * pq.ms)
>>> rate = statistics.instantaneous_rate(spiketrain,
...        sampling_period=200 * pq.ms, kernel=kernel, t_start=-1 * pq.s)
>>> rate
<AnalogSignal(array([[0.01007419],
       [0.05842767],
       [0.22928759],
       [0.60883028],
       [1.0938699 ],
       [1.3298076 ],
       [1.0938699 ],
       [0.60883028],
       [0.22928759],
       [0.05842767]]) * Hz, [-1.0 s, 1.0 s], sampling rate: 0.005 1/ms)>

>>> rate.magnitude
array([[0.01007419],
       [0.05842767],
       [0.22928759],
       [0.60883028],
       [1.0938699 ],
       [1.3298076 ],
       [1.0938699 ],
       [0.60883028],
       [0.22928759],
       [0.05842767]])