import sys from numpy import * import matplotlib.pyplot as plt #Initialise T = 100 #length of trial r = 10 #rate of firing in Hz if (len(sys.argv) == 3): T = int(sys.argv[1]) r = int(sys.argv[2]) print "Using given values. Length = %is, Rate = %iHz" % (T, r) else: print "Using default values. Length = %is, Rate = %iHz" % (T, r) #Generate Data intervals = [] while (sum(intervals) < T): #add random intervals using Poisson distribution until over length of trial intervals.append(-log(random.rand())/r) intervals.pop() #remove last one as it's over the length of the trial #Output Results print "Number of spikes: %i (expected approx. %i)" % (len(intervals), r*T) #Coefficient of variation = deviation of interval / average interval CV = std(intervals) / mean(intervals) print "Coefficient of Variation (should be approx 1) = %f" % CV sp_times = cumsum(intervals) #ordered array of each spike time. #Plotting the generated spike chain. y_axis = [1] * len(sp_times) plt.bar(sp_times, y_axis, width=(1/(float(T*r)))) plt.show()
The results show the problems with not modelling the refractory period with many spikes very close when rate is set high.
Python script for checking generation leads to Gaussian distribution for number of spikes over the time period:
import sys from numpy import * import matplotlib.pyplot as plt #Initialise T = 10 #length of trial r = 10 #rate of firing in Hz trials = 10000 if (len(sys.argv) >= 3): T = int(sys.argv[1]) r = int(sys.argv[2]) if (len(sys.argv) >= 4): trials = int(sys.argv[3]) print "Using given values. Length = %is, Rate = %iHz for %i trials" % (T, r, trials) else: print "Using default values. Length = %is, Rate = %iHz for %i trials" % (T, r, trials) #Generate Data no_spikes = [] while (len(no_spikes) < trials): if(len(no_spikes) % 10 == 0): print '{0}{1}{2}\r'.format(len(no_spikes)/float(trials)*100, "%", " completed..."), intervals = [] while (sum(intervals) < T): #add random intervals using Poisson distribution until over length of trial intervals.append(-log(random.rand())/r) no_spikes.append(len(intervals)-1) #Output Results print "Average number of spikes: %f (expected approx. %i)" % (mean(no_spikes), r*T) density = [] for i in range(0, max(no_spikes)+1): density.append(no_spikes.count(i)) print density plt.plot(density, 'o') plt.show()
Sample output: