## ranksurprise.R --- Rank Surprise method for burst detection.
## Author: Zhengzheng Zhang
## Copyright: GPL.
## Following method in Gourevitch and Eggermont (2007)
## "Rank Surprise" method for burst detection. The method and algorithm are
## described in the paper, a nonparametric approach for detection of bursts
## in spike trains. See reference.
val2rk <- function(values) {
## Convert values to ranks, with mean of ranks for tied values. Alternative
## and faster version of "tiedrank" in statistical toolbox of Matlab 6-7.
lp = length(values)
rk = rk2 = rep(0,lp)
S = sort(values, index.return = TRUE)
y = S$x
cl = S$ix
rk[cl] = c(1:lp)
cl2 = sort(-values, index.return = TRUE)$ix
rk2[cl2] = c(1:lp)
ranks = (lp+1-rk2+rk)/2
ranks
}
rs.find.bursts <- function(s, limit = NULL, RSalpha = -log(0.01)) {
nspikes = length(s)
if (nspikes <= 10) {
archive_burst = NA
}else{
## General parameters
## limit for using the real distribution
q_lim = 30
## minimum length of a burst in spikes
l_min = 3
## General vectors
## vector (-1)^k
alternate = rep(1, 400)
alternate[c(1:200)*2] = -1
## log factorials
log_fac = cumsum(log(1:q_lim))
## Ranks computation
## compute the ISI
ISI = diff(s)
N = length(ISI)
## ISI value not to include in a burst
if (length(limit)==0){
## percentile 75% (default)
limit=summary(ISI)["3rd Qu."][[1]]
}
## compute ranks
R = val2rk(ISI)
## Find sequences of ISI under 'limit'
D = rep(0, N)
D[ISI<limit]=1
ISI_limit = diff(D)
## first time stamp of these intervals
begin_int = which(ISI_limit==1)+1
## manage the first ISI
if (ISI[1]<limit){
begin_int = c(1, begin_int) # the first IS is under limit
}
## last time stamp of these intervals
end_int = which(ISI_limit==-1)
## manage the last ISI
if (length(end_int)<length(begin_int)){
end_int=c(end_int, N)
}
## length of intervals of interest
length_int=end_int-begin_int+1
## Initializations
archive_burst_RS=c()
archive_burst_length=c()
archive_burst_start=c()
archive_burst_end=c()
archive_burst_IBI=c()
archive_burst_durn=c()
archive_burst_mean.isis=c()
## Going through the intervals of interest
indic=0
for (n_j in begin_int){
indic=indic+1
p_j=length_int[indic]
subseq_RS = NULL
subseq_RS_RS=c()
subseq_RS_i =c()
subseq_RS_q =c()
## test each set of spikes
for (i in 0:(p_j-(l_min-1))){
## length of burst tested
q=l_min-2
while (q<(p_j-i)){
q=q+1
## statistic
u=sum(R[(n_j+i):(n_j+i+q-1)])
u=floor(u)
if (q<q_lim){
## exact discrete distribution
k=c(0:((u-q)/N))
length_k=length(k)
SUM = 0
for (j in 1:q){
SUM = SUM + log(u-matrix(rep(k,q),q,length_k, byrow = TRUE)*N
- matrix(rep(c(0:(q-1)),length_k),length(c(0:(q-1))),length_k))[j,]
}
if (length_k<2){
prob = exp((SUM- log_fac[1]-log_fac[q-k])-q*log(N))%*%alternate[1:length_k]
}else{
prob = exp((SUM- log_fac[c(1, k[2:length_k])]-log_fac[q-k])-q*log(N))%*%alternate[1:length_k]
}
}else{
## approximate Gaussian distribution
prob=pnorm((u-q*(N+1)/2)/sqrt(q*(N^2-1)/12))
}
RS=-log(prob)
## archive results for each subsequence [RSstatistic beginning length]
if (RS>RSalpha){
subseq_RS_RS[length(subseq_RS_RS)+1]=RS
subseq_RS_i[length(subseq_RS_i)+1]=i
subseq_RS_q[length(subseq_RS_q)+1]=q
subseq_RS = rbind(subseq_RS_RS, subseq_RS_i, subseq_RS_q)
subseq_RS = t(as.matrix(subseq_RS))
}
}
}
## vet results archive to extract most significant bursts
if (length(subseq_RS)!=0){
## sort RS for all subsequences
if (!is.vector(subseq_RS)){
ind = sort(subseq_RS[,1], decreasing = TRUE, index = TRUE)$ix
subseq_RS = subseq_RS[ind,]
}
while (length(subseq_RS)!=0){
## extract most surprising burst
if (is.vector(subseq_RS)){
current_burst=subseq_RS
archive_burst_RS[length(archive_burst_RS)+1]=current_burst[1]
archive_burst_length[length(archive_burst_length)+1]=current_burst[3]+1 #number of ISI involved + 1
archive_burst_start[length(archive_burst_start)+1]=n_j+current_burst[2]
archive_burst_end[length(archive_burst_end)+1]=n_j+current_burst[2]+current_burst[3]
subseq_RS=NULL
}else{
current_burst=subseq_RS[1,]
archive_burst_RS[length(archive_burst_RS)+1]=current_burst[1]
archive_burst_length[length(archive_burst_length)+1]=current_burst[3]+1 #number of ISI involved + 1
archive_burst_start[length(archive_burst_start)+1]=n_j+current_burst[2]
archive_burst_end[length(archive_burst_end)+1]=n_j+current_burst[2]+current_burst[3]
## remove most surprising burst from the set
## subseq_RS=subseq_RS(2:end,:);
## keep only other bursts non-overlapping with this burst
D = (subseq_RS[,2]+subseq_RS[,3]-1)<current_burst[2]
E = subseq_RS[,2]>(current_burst[2]+current_burst[3]-1)
subseq_RS=subseq_RS[D|E,]
}
}
}
}
## sort bursts by ascending time
if (length(archive_burst_RS)!=0){
S = sort(archive_burst_start, index.return = TRUE)
archive_burst_start = S$x
ind_sort = S$ix
archive_burst_RS=archive_burst_RS[ind_sort]
archive_burst_length=archive_burst_length[ind_sort]
archive_burst_end=archive_burst_end[ind_sort]
last.end = archive_burst_end[-length(archive_burst_end)]
last.end = c(NA, s[last.end])
archive_burst_IBI = s[archive_burst_start]-last.end
archive_burst_durn = s[archive_burst_end]-s[archive_burst_start]
archive_burst_mean.isis = archive_burst_durn/(archive_burst_length-1)
archive_burst=cbind(beg=archive_burst_start, end=archive_burst_end, IBI=archive_burst_IBI,
len=archive_burst_length, durn=archive_burst_durn,
mean.isis=archive_burst_mean.isis, SI=archive_burst_RS)
}else{
archive_burst = NA
}
}
archive_burst
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.