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

Question

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

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:

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? Any help is appreciated, even if it's just about how to start/where to look.

CPak · Accepted Answer · 2018-04-16T12:41:15.323

1

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)

edited Apr 16 '18 at 12:41

answered Apr 10 '18 at 13:53

CPak

13,260
3
30
48

Hi @CPak , thank you so much for your time and efforts. I am very new at R and all of the functions you have used, I haven't used before. I am trying to make sense of it. Could you please elaborate on the lapply function and how is it working? – 14thTimeLord Apr 11 '18 at 12:53
`lapply` is like `for` loop. `for (i in 1:N) { do stuff }` is equal to `lapply(1:N, function(i) { do stuff }`. `seq_len(N)` is a function that returns `1:N`. I use `setNames({values}, {names})` to generate a named vector. Hope that helps a little. Obviously, you're expected to have to a little Google-fu on your own but let me know if something is super confusing. – CPak Apr 11 '18 at 12:55
Yes of course @CPak , I went through the help section of R, and `seq_along` was the only part that didn't make complete sense to me, as to how is it contributing. Also the `nearest_points` you have mentioned in `nearest_columns`; I am sorry if I sound naive but we didn't define that during our code and I can't find anything related in help because it's not a function. Again, I have started really recently and so many things that might seem obvious to you, are some times very confusing to me. Hope I will be better in this, soon. – 14thTimeLord Apr 11 '18 at 13:03
`seq_along(list)` or `seq_along(vector)` returns `1:length(list)` or `1:length(vector)`. `nearest_points` should have been `distances`. I changed the names in the middle of writing the code. Let me know if you're getting an error. – CPak Apr 11 '18 at 14:53
thank you. The code ran perfectly. `sapply` gives me: `6 35 ... 9 15 10` i.e. no of bins between successive 1's in the neighborhood of `RT` . For confirmation I took values of `nearest_columns`. `RT` has the value **531**, so in theory it should fall between column no **(531-x)** and **(531+x)**. Eg for the last trial it falls between 524 and 534 and the value after `sapply` is 10 but for the trial before that it falls between 526 and 521 and `sapply` gives 15. And throughout, the `sapply` result only matches `nearest_columns` where the value of falls between 531(-) and 531(+). Please help. – 14thTimeLord Apr 16 '18 at 12:02
My mistake. I assumed that the nearest 2 values would always be bookending, but they might both be to the left or right of `RT`. Updated my answer. Let me know if it doesn't look right. – CPak Apr 16 '18 at 12:40
Ah okay, I get what was going wrong. Thank you for the code. It works perfectly and confirms my results too :) – 14thTimeLord Apr 17 '18 at 09:09

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

1 Answers1

Linked