Skip to content

Commit

Permalink
[Feature] Add revised local variability LvR (#346)
Browse files Browse the repository at this point in the history 
* Add LvR statistic, move the order of statistics to have cv, cv2, lv, lvr,
change CV -> Cv, and LV -> Lv, in agreement with the notation by 
Shinomoto et al.
  • Loading branch information
@morales-gregorio
morales-gregorio authored Sep 28, 2020
1 parent 5d43e0a commit f171f28
Show file tree
Hide file tree
Showing 2 changed files with 174 additions and 35 deletions.
156 changes: 121 additions & 35 deletions elephant/statistics.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -40,8 +40,9 @@
isi
cv
lv
cv2
lv
lvr
Statistics across spike trains
Expand Down Expand Up @@ -81,9 +82,10 @@
"isi",
"mean_firing_rate",
"fanofactor",
"lv",
"cv",
"cv2",
"lv",
"lvr",
"instantaneous_rate",
"time_histogram",
"complexity_pdf",
Expand Down Expand Up @@ -322,22 +324,91 @@ def __variation_check(v, with_nan):
return None


@deprecated_alias(v='time_intervals')
def cv2(time_intervals, with_nan=False):
r"""
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
----------
time_intervals : pq.Quantity or np.ndarray or list
Vector of consecutive time intervals.
with_nan : bool, optional
If True, `cv2` of a spike train with less than two spikes results in a
np.NaN value and a warning is raised.
If False, `ValueError` exception is raised with a spike train with
less than two spikes.
Default: True.
Returns
-------
float
The Cv2 of the inter-spike interval of the input sequence.
Raises
------
ValueError
If an empty list is specified, or if the sequence has less than two
entries and `with_nan` is False.
If a matrix is passed to the function. Only vector inputs are
supported.
Warns
-----
UserWarning
If `with_nan` is True and `cv2` is calculated for a sequence with less
than two entries, generating a np.NaN.
References
----------
.. [1] G. R. Holt, W. R. Softky, C. Koch, & R. J. Douglas, "Comparison of
discharge variability in vitro and in vivo in cat visual cortex
neurons," Journal of Neurophysiology, vol. 75, no. 5, pp. 1806-1814,
1996.
"""
# convert to array, cast to float
time_intervals = np.asarray(time_intervals)
np_nan = __variation_check(time_intervals, with_nan)
if np_nan is not None:
return np_nan

# calculate Cv2 and return result
cv_i = np.diff(time_intervals) / (time_intervals[:-1] + time_intervals[1:])
return 2. * np.mean(np.abs(cv_i))


@deprecated_alias(v='time_intervals')
def lv(time_intervals, with_nan=False):
r"""
Calculate the measure of local variation LV for a sequence of time
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:
Given a vector I 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}
Lv := \frac{1}{N} \sum_{i=1}^{N-1}
\frac{3(I_i-I_{i+1})^2}
{(I_i+I_{i+1})^2}
The LV is typically computed as a substitute for the classical coefficient
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
rate fluctuation. As with the Cv, Lv=1 for a sequence of intervals
generated by a Poisson process.
Parameters
Expand All @@ -354,7 +425,7 @@ def lv(time_intervals, with_nan=False):
Returns
-------
float
The LV of the inter-spike interval of the input sequence.
The Lv of the inter-spike interval of the input sequence.
Raises
------
Expand Down Expand Up @@ -388,40 +459,44 @@ def lv(time_intervals, with_nan=False):
return 3. * np.mean(np.power(cv_i, 2))


@deprecated_alias(v='time_intervals')
def cv2(time_intervals, with_nan=False):
def lvr(time_intervals, R=5*pq.ms, with_nan=False):
r"""
Calculate the measure of CV2 for a sequence of time intervals between
events.
Calculate the measure of revised local variation LvR for a sequence of time
intervals between events.
Given a vector v containing a sequence of intervals, the CV2 is defined
as:
Given a vector :math:`I` containing a sequence of intervals, the LvR 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|}
LvR := \frac{3}{N-1} \sum_{i=1}^{N-1}
\left(1-\frac{4 I_i I_{i+1}}
{(I_i+I_{i+1})^2}\right)
\left(1+\frac{4 R}{I_i+I_{i+1}}\right)
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.
The LvR is a revised version of the Lv, with enhanced invariance to firing
rate fluctuations by introducing a refractoriness constant R. The LvR with
`R=5ms` was shown to outperform other ISI variability measures in spike
trains with firing rate fluctuatins and sensory stimuli [1]_.
Parameters
----------
time_intervals : pq.Quantity or np.ndarray or list
Vector of consecutive time intervals.
R : pq.Quantity or int or float
Refractoriness constant (R >= 0). If no quantity is passed `ms` are
assumed.
Default: 5 ms.
with_nan : bool, optional
If True, `cv2` of a spike train with less than two spikes results in a
If True, LvR of a spike train with less than two spikes results in a
np.NaN value and a warning is raised.
If False, `ValueError` exception is raised with a spike train with
If False, a `ValueError` exception is raised with a spike train with
less than two spikes.
Default: True.
Returns
-------
float
The CV2 of the inter-spike interval of the input sequence.
The LvR of the inter-spike interval of the input sequence.
Raises
------
Expand All @@ -435,26 +510,37 @@ def cv2(time_intervals, with_nan=False):
Warns
-----
UserWarning
If `with_nan` is True and `cv2` is calculated for a sequence with less
than two entries, generating a np.NaN.
If `with_nan` is True and the `lvr` is calculated for a spike train
with less than two spikes, generating a np.NaN.
If R is passed without any units attached milliseconds are assumed.
References
----------
.. [1] G. R. Holt, W. R. Softky, C. Koch, & R. J. Douglas, "Comparison of
discharge variability in vitro and in vivo in cat visual cortex
neurons," Journal of Neurophysiology, vol. 75, no. 5, pp. 1806-1814,
1996.
.. [1] S. Shinomoto, H. Kim, T. Shimokawa et al. "Relating Neuronal Firing
Patterns to Functional Differentiation of Cerebral Cortex" PLOS
Computational Biology 5(7): e1000433, 2009.
"""
if isinstance(R, pq.Quantity):
R = R.rescale('ms').magnitude
else:
warnings.warn('No units specified for R, assuming milliseconds (ms)')

if R < 0:
raise ValueError('R must be >= 0')

# convert to array, cast to float
time_intervals = np.asarray(time_intervals)
np_nan = __variation_check(time_intervals, with_nan)
if np_nan is not None:
return np_nan

# calculate CV2 and return result
cv_i = np.diff(time_intervals) / (time_intervals[:-1] + time_intervals[1:])
return 2. * np.mean(np.abs(cv_i))
N = len(time_intervals)
t = time_intervals[:-1] + time_intervals[1:]
frac1 = 4 * time_intervals[:-1] * time_intervals[1:] / t**2
frac2 = 4 * R / t
lvr = (3 / (N-1)) * np.sum((1-frac1) * (1+frac2))
return lvr


def instantaneous_rate(spiketrain, sampling_period, kernel='auto',
Expand Down
53 changes: 53 additions & 0 deletions elephant/test/test_statistics.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -391,6 +391,59 @@ def test_2short_spike_train(self):
self.assertTrue(math.isnan(statistics.lv(seq, with_nan=True)))


class LVRTestCase(unittest.TestCase):
def setUp(self):
self.test_seq = [1, 28, 4, 47, 5, 16, 2, 5, 21, 12,
4, 12, 59, 2, 4, 18, 33, 25, 2, 34,
4, 1, 1, 14, 8, 1, 10, 1, 8, 20,
5, 1, 6, 5, 12, 2, 8, 8, 2, 8,
2, 10, 2, 1, 1, 2, 15, 3, 20, 6,
11, 6, 18, 2, 5, 17, 4, 3, 13, 6,
1, 18, 1, 16, 12, 2, 52, 2, 5, 7,
6, 25, 6, 5, 3, 15, 4, 3, 16, 3,
6, 5, 24, 21, 3, 3, 4, 8, 4, 11,
5, 7, 5, 6, 8, 11, 33, 10, 7, 4]

self.target = 2.1845363464753134

def test_lvr_with_quantities(self):
seq = pq.Quantity(self.test_seq, units='ms')
assert_array_almost_equal(statistics.lvr(seq), self.target, decimal=9)

def test_lvr_with_plain_array(self):
seq = np.array(self.test_seq)
assert_array_almost_equal(statistics.lvr(seq), self.target, decimal=9)

def test_lvr_with_list(self):
seq = self.test_seq
assert_array_almost_equal(statistics.lvr(seq), self.target, decimal=9)

def test_lvr_raise_error(self):
seq = self.test_seq
self.assertRaises(ValueError, statistics.lvr, [])
self.assertRaises(ValueError, statistics.lvr, 1)
self.assertRaises(ValueError, statistics.lvr, np.array([seq, seq]))
self.assertRaises(ValueError, statistics.lvr, seq, -1)

@unittest.skipUnless(python_version_major == 3, "assertWarns requires 3.2")
def test_lvr_refractoriness_kwarg(self):
seq = np.array(self.test_seq)
with self.assertWarns(UserWarning):
assert_array_almost_equal(statistics.lvr(seq, R=5),
self.target, decimal=9)

@unittest.skipUnless(python_version_major == 3, "assertWarns requires 3.2")
def test_2short_spike_train(self):
seq = [1]
with self.assertWarns(UserWarning):
"""
Catches UserWarning: Input size is too small. Please provide
an input with more than 1 entry.
"""
self.assertTrue(math.isnan(statistics.lvr(seq, with_nan=True)))



class CV2TestCase(unittest.TestCase):
def setUp(self):
self.test_seq = [1, 28, 4, 47, 5, 16, 2, 5, 21, 12,
Expand Down

0 comments on commit f171f28

Please sign in to comment.