Last active
May 13, 2017 16:25
-
-
Save leon-nn/66f046c4ba829f39850c70f0c2484834 to your computer and use it in GitHub Desktop.
Nonhomogenous Poisson spiking train in a spiking neural network framework
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
%% Non-homogenous Poisson spike train in a spiking neural network framework | |
% We generate inhomogenous Poisson spike trains for a set of output neurons | |
% in a spiking neural network. We use the "thinning" approach introduced in | |
% "Simulation of nonhomogeneous poisson processes by thinning" by Lewis and | |
% Shedler (1979). | |
% Number of output neurons | |
numOut = 10; | |
% Define a supremum for lambda(t). I.e., this should be the smallest upper | |
% bound on the values that lambda(t) can take. | |
lambdaMax = 1; | |
% Runtime of simulation | |
tSim = 2500; | |
% Initialize (but note that we don't need to store these because we only | |
% depend on their values at the current time t) | |
lambda = zeros(numOut, tSim); | |
u = zeros(numOut, tSim); | |
% Inhibitory function | |
I = 2000 * exp(-5 * (0: tSim + 4)') + 550; | |
% Initialize a time counter for the inhibitory function | |
tInh = 5 * ones(numOut, 1); | |
% Initialize the time of next spike for the output neurons | |
tNext = tSim * ones(numOut, 1); | |
for t = 1: tSim | |
% Generate EPSP u(t) | |
% YOUR CODE HERE | |
% If there is a spike, reset the time of the inhibitory function for | |
% the neurons that did not spike | |
if any(tNext == t) | |
% Find the index of the first neuron that is spiking this time t | |
ind = find(tNext == t, 1); | |
% Reset the time of the inhibitory function on the other neurons | |
tInh([1:ind-1, ind+1:10]) = 1; | |
end | |
% Poisson rate at time t for each of the output neurons | |
lambda(:, t) = exp(u(:, t) - I(tInh)); | |
% Find the indices to update: if the next spiking time is the default | |
% initialization or the current time, then we have to calculate a new | |
% next spiking time. | |
update = tNext == tSim | tNext == t; | |
% Find the next time of spike assuming a homogenous Poisson rate | |
% lambdaMax. This uses inverse transform sampling of the exponential | |
% distribution to find the time interval to the next Poisson event. | |
% Note that the time between Poisson events is exponentially | |
% distributed. | |
tNext(update) = t - log(rand(sum(update), 1)) / lambdaMax; | |
% Reject the next time of spike if a Unif ~ [0, 1] random number is | |
% greater than the ratio lambda(t)/lambdaMax | |
tNext(update & (rand(numOut, 1) > lambda(:, t) / lambdaMax)) = tSim; | |
% Perform postsynaptic potential weight update | |
if any(tNext == t) | |
% YOUR CODE HERE | |
end | |
% Update time on inhibition function | |
tInh = tInh + 1; | |
end |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment