Source code for elephant.cubic
# -*- coding: utf-8 -*-
'''
CuBIC is a statistical method for the detection of higher order of
correlations in parallel spike trains based on the analysis of the
cumulants of the population count.
Given a list sts of SpikeTrains, the analysis comprises the following
steps:
1) compute the population histogram (PSTH) with the desired bin size
>>> binsize = 5 * pq.ms
>>> pop_count = elephant.statistics.time_histogram(sts, binsize)
2) apply CuBIC to the population count
>>> alpha = 0.05 # significance level of the tests used
>>> xi, p_val, k = cubic(data, ximax=100, alpha=0.05, errorval=4.):
:copyright: Copyright 2016 by the Elephant team, see AUTHORS.txt.
:license: BSD, see LICENSE.txt for details.
'''
# -*- coding: utf-8 -*-
from __future__ import division
import scipy.stats
import scipy.special
import math
import warnings
# Based on matlab code by Benjamin Staude
# Adaptation to python by Pietro Quaglio and Emiliano Torre
[docs]def cubic(data, ximax=100, alpha=0.05):
'''
Performs the CuBIC analysis [1] on a population histogram, calculated from
a population of spiking neurons.
The null hypothesis :math:`H_0: k_3(data)<=k^*_{3,\\xi}` is iteratively
tested with increasing correlation order :math:`\\xi` (correspondent to
variable xi) until it is possible to accept, with a significance level alpha,
that :math:`\\hat{\\xi}` (corresponding to variable xi_hat) is the minimum
order of correlation necessary to explain the third cumulant
:math:`k_3(data)`.
:math:`k^*_{3,\\xi}` is the maximized third cumulant, supposing a Compund
Poisson Process (CPP) model for correlated spike trains (see [1])
with maximum order of correlation equal to :math:`\\xi`.
Parameters
----------
data : neo.AnalogSignal
The population histogram (count of spikes per time bin) of the entire
population of neurons.
ximax : int
The maximum number of iteration of the hypothesis test:
if it is not possible to compute the :math:`\\hat{\\xi}` before ximax
iteration the CuBIC procedure is aborted.
Default: 100
alpha : float
The significance level of the hypothesis tests perfomed.
Default: 0.05
Returns
-------
xi_hat : int
The minimum correlation order estimated by CuBIC, necessary to
explain the value of the third cumulant calculated from the population.
p : list
The ordred list of all the p-values of the hypothesis tests that have
been performed. If the maximum number of iteration ximax is reached the
last p-value is set to -4
kappa : list
The list of the first three cumulants of the data.
test_aborted : bool
Wheter the test was aborted because reached the maximum number of
iteration ximax
References
----------
[1]Staude, Rotter, Gruen, (2009) J. Comp. Neurosci
'''
# alpha in in the interval [0,1]
if alpha < 0 or alpha > 1:
raise ValueError(
'the significance level alpha (= %s) has to be in [0,1]' % alpha)
if not isinstance(ximax, int) or ximax < 0:
raise ValueError(
'The maximum number of iterations ximax(= %i) has to be a positive'
% alpha + ' integer')
# dict of all possible rate functions
try:
data = data.magnitude
except AttributeError:
pass
L = len(data)
# compute first three cumulants
kappa = _kstat(data)
xi_hat = 1
xi = 1
pval = 0.
p = []
test_aborted = False
# compute xi_hat iteratively
while pval < alpha:
xi_hat = xi
if xi > ximax:
warnings.warn('Test aborted, xihat= %i > ximax= %i' % (xi, ximax))
test_aborted = True
break
# compute p-value
pval = _H03xi(kappa, xi, L)
p.append(pval)
xi = xi + 1
return xi_hat, p, kappa, test_aborted
def _H03xi(kappa, xi, L):
'''
Computes the p_value for testing the :math:`H_0: k_3(data)<=k^*_{3,\\xi}`
hypothesis of CuBIC in the stationary rate version
Parameters
-----
kappa : list
The first three cumulants of the populaton of spike trains
xi : int
The the maximum order of correlation :math:`\\xi` supposed in the
hypothesis for which is computed the p value of :math:`H_0`
L : float
The length of the orginal population histogram on which is performed
the CuBIC analysis
Returns
-----
p : float
The p-value of the hypothesis tests
'''
# Check the order condition of the cumulants necessary to perform CuBIC
if kappa[1] < kappa[0]:
# p = errorval
kstar = [0]
raise ValueError(
'H_0 can not be tested:'
'kappa(2)= %f<%f=kappa(1)!!!' % (kappa[1], kappa[0]))
else:
# computation of the maximized cumulants
kstar = [_kappamstar(kappa[:2], i, xi) for i in range(2, 7)]
k3star = kstar[1]
# variance of third cumulant (from Stuart & Ord)
sigmak3star = math.sqrt(
kstar[4] / L + 9 * (kstar[2] * kstar[0] + kstar[1] ** 2) /
(L - 1) + 6 * L * kstar[0] ** 3 / ((L - 1) * (L - 2)))
# computation of the p-value (the third cumulant is supposed to
# be gaussian istribuited)
p = 1 - scipy.stats.norm(k3star, sigmak3star).cdf(kappa[2])
return p
def _kappamstar(kappa, m, xi):
'''
Computes maximized cumulant of order m
Parameters
-----
kappa : list
The first two cumulants of the data
xi : int
The :math:`\\xi` for which is computed the p value of :math:`H_0`
m : float
The order of the cumulant
Returns
-----
k_out : list
The maximized cumulant of order m
'''
if xi == 1:
kappa_out = kappa[1]
else:
kappa_out = \
(kappa[1] * (xi ** (m - 1) - 1) -
kappa[0] * (xi ** (m - 1) - xi)) / (xi - 1)
return kappa_out
def _kstat(data):
'''
Compute first three cumulants of a population count of a population of
spiking
See http://mathworld.wolfram.com/k-Statistic.html
Parameters
-----
data : numpy.aray
The population histogram of the population on which are computed
the cumulants
Returns
-----
kappa : list
The first three cumulants of the population count
'''
L = len(data)
if L == 0:
raise ValueError('The input data must be a non-empty array')
S = [(data ** r).sum() for r in range(1, 4)]
kappa = []
kappa.append(S[0] / float(L))
kappa.append((L * S[1] - S[0] ** 2) / (L * (L - 1)))
kappa.append(
(2 * S[0] ** 3 - 3 * L * S[0] * S[1] + L ** 2 * S[2]) / (
L * (L - 1) * (L - 2)))
return kappa