1

Forgive me if this seems obvious, I am VERY new to R. So I am trying to get random Spike Train that looks like a vector [ 1 0 1 1 0 0 0 1] and for one instance, I have been able to do it using the following code:

fr = 100 #Firing rate of 100Hz
dt = 1/1000 #Short duration of time dt
nBins = 10 #10msSpikeTrain 
x = runif(nBins) #Creating a "nBins" length string of random variable 
#falling uniformly between 0 and 1
x
y<- numeric(nBins) # creating an empty vector of size nBins
MyPoissonSpikeTrain = function(x){
for (i in 1:nBins){
if (x[i] < fr*dt)
y[i]=1
else
y[i]=0
}
return(y)
}  
#Creating a function that that returns its values in vector y as 1 if the 
#spike was fired 
#and 0 if the spike wasn't fired. 
#Spike was fired, if randomly generated x value is smaller than fr*dt. 
MyPoissonSpikeTrain(x)

This works all well, but I want to create a matrix where 1st row would be one instance of the above procedure, 2nd row would be 2nd instance and so on.

I want this process happening (say) 20 times, and then put each of the individual y I find as a row of a Bigger Matrix, say Y. I understand that I need to modify my y to be a matrix and not a vector and then modify my MyPoissonSpikeTrain function accordingly to have the output create a matrix but I can't figure out how to do it. Any help will be appreciated.

14thTimeLord
  • 363
  • 1
  • 14

1 Answers1

0

This function can be simplified a lot, check this out:

MyPoissonSpikeTrain = function(fr = 100, dt = 1/1000, nBins = 10){
  x = runif(nBins)
  y = ifelse(x < dt * fr, 0, 1)
  return(y)
}  

now it is vectorized, and you can change the arguments if you wish. Now to generate the desired matrix:

set.seed(1)
replicate(10, MyPoissonSpikeTrain())
#output
      [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
 [1,]    1    1    1    1    1    1    1    1    1     1
 [2,]    1    1    1    1    1    1    1    1    1     0
 [3,]    1    1    1    1    1    1    1    1    1     1
 [4,]    1    1    1    1    1    1    1    1    1     1
 [5,]    1    1    1    1    1    0    1    1    1     1
 [6,]    1    1    1    1    1    0    1    1    1     1
 [7,]    1    1    0    1    0    1    1    1    1     1
 [8,]    1    1    1    1    1    1    1    1    1     1
 [9,]    1    1    1    1    1    1    0    1    1     1
[10,]    0    1    1    1    1    1    1    1    1     1

Here each column is one repetition, if you want it the other way round use: t

t(replicate(10, MyPoissonSpikeTrain()))

If you want to change some argument:

set.seed(1)
replicate(10, MyPoissonSpikeTrain(fr = 500))
#output
      [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
 [1,]    1    1    1    1    1    1    0    1    0     1
 [2,]    0    1    0    0    1    1    1    1    0     0
 [3,]    0    0    0    0    0    0    0    1    1     0
 [4,]    1    0    0    1    0    0    1    0    1     1
 [5,]    1    0    1    1    1    1    1    0    1     0
 [6,]    0    0    0    1    0    0    0    1    0     1
 [7,]    0    1    1    1    0    1    0    1    1     0
 [8,]    0    0    0    1    1    0    0    1    1     1
 [9,]    1    0    0    1    0    0    1    1    1     0
[10,]    1    1    1    1    1    0    0    1    1     1

or even better:

MyPoissonSpikeTrain = function(fr = 100, dt = 1/1000, nBins = 10, ...){
  x = runif(nBins, ...)
  y = ifelse(x < dt * fr, 0, 1)
  return(y)
}  

so you can pass additional argument to runif

set.seed(1)
replicate(10, MyPoissonSpikeTrain(max = 0.5))
missuse
  • 19,056
  • 3
  • 25
  • 47
  • Hi @missuse , thank you so much for your help! I checked out your initial code (the simplification), it didn't give me the desired result (it returned a string of 1's no matter what the input x was). However, I figured out the "mistake" if you can call it that. The simplified function that we are using actually doesn't have a variable. So I combined your code and mine and mentioned the constants (fr, dt, nBins) and then made MyPoissonSpikeTrain a function of x but used the ifelse operator for the output. It's shorter but is also working well. I will try the rest of the code now. Thanks again. – 14thTimeLord Mar 28 '18 at 09:15
  • @R Tomar I just copied the function from my post to a clean R session and the output is not all 1, all the time. But it is mostly one since it is testing how many random samples from the uniform distribution 0 - 1 is less than 0.1. This is 10%, so usually there are 0-2 zeros. – missuse Mar 28 '18 at 09:43
  • I tested it with different dt, which made the limit around 0.5 because I had the same doubt. But maybe I did something wrong in copying :) – 14thTimeLord Mar 28 '18 at 11:00
  • Rest of it is working wonderfully, thank you so much for your time and efforts. I really appreciate it @missuse – 14thTimeLord Mar 28 '18 at 11:22