Search code examples
rplottime-seriesdistribution

How to get a random observation point at a specific time over multiple trials in R?

I am working on Spike Trains and my code to get a spike train like this: SpikeTrains

for 20 trials is written below. The image is representational for 5 trials.

fr = 100
dt = 1/1000 #dt in milisecond
duration = 2 #no of duration in s
nBins = 2000 #10msSpikeTrain
nTrials = 20 #NumberOfSimulations
MyPoissonSpikeTrain = function(p, fr= 100) {
 p = runif(nBins)
 q = ifelse(p < fr*dt, 1, 0)
 return(q)
  }

set.seed(1)
SpikeMat <- t(replicate(nTrials, MyPoissonSpikeTrain()))

plot(x=-1,y=-1, xlab="time (s)", ylab="Trial",
    main="Spike trains",
    ylim=c(0.5, nTrials+1), xlim=c(0, duration))
for (i in 1: nTrials)
 {
   clip(x1 = 0, x2= duration, y1= (i-0.2), y2= (i+0.4))
   abline(h=i, lwd= 1/4)
   abline(v= dt*which( SpikeMat[i,]== 1))
  }

Each trial has spikes occuring at random time points. Now what I am trying to work towards, is getting a random sample time point that works for all 20 trials and I want to get the vector consisting of length of the intervals this point falls into, for each trial. The code to get the time vector for the points where the spikes occur is,

A <- numeric()
for (i in 1: nTrials)
 {
  ISI <- function(i){
   spike_times <- c(dt*which( SpikeMat[i, ]==1))
   ISI1vec <- c(diff(spike_times))
   A <- c(A, ISI1vec)
   return(A)}
    }

Then you call ISI(i) for whichever trial you wish to see the Interspike interval vector for. A visual representation of what I want is:

NewSpikeTrain

I want to get a vector that has the lengths of the interval where this points fall into, for each trial. I want to figure out it's distribution as well, but that's for later. Can anybody help me figure out how to code my way to this?

Solution

  • Your data

    set.seed(1)
SpikeMat <- t(replicate(nTrials, MyPoissonSpikeTrain()))

    I suggest transforming your sparse matrix data into a list of indices where spikes occur

    L <- lapply(seq_len(nrow(SpikeMat)), function(i) setNames(which(SpikeMat[i, ] == 1), seq_along(which(SpikeMat[i, ] == 1))))

    Grab random timepoint

    set.seed(1)
RT <- round(runif(1) * ncol(SpikeMat)) 
# 531

    Result

    distances contains the distances to the 2 nearest spikes - each element of the list is a named vector where the values are the distances (to RT) and their names are their positions in the vector. nearest_columns shows the original timepoint (column number) of each spike in SpikeMat. 

    bookend_values <- function(vec) {
    lower_val <- head(sort(vec[sign(vec) == 1]), 1)
    upper_val <- head(sort(abs(vec[sign(vec) == -1])), 1)
    return(c(lower_val, upper_val))
}

distances <- lapply(L, function(i) bookend_values(RT-i))
nearest_columns <- lapply(seq_along(distances), function(i) L[[i]][names(distances[[i]])])

    Note that the inter-spike interval of the two nearest spikes that bookend RT can be obtained with

    sapply(distances, sum)