% Poisson_neuron

% APMA 0410
% Fall 2012

% Last updated October 5, 2012

clear
close all

duration_of_recording = 3;              % [sec]
delta_t = 10^(-6);                      % time resolution [sec}
r = 10;                                 % firing rate [spikes/sec]

N_bins = duration_of_recording/delta_t; % number of time intervals
time_axis = 1:N_bins;
prob_spike_per_bin = r*delta_t;

random_numbers = rand(1,N_bins);
spikes = random_numbers < prob_spike_per_bin;
N_spikes = sum(spikes);                 % number of spikes
figure(1)
plot(time_axis,spikes,'-')
set(gca, 'fontsize',16)
axis([0 N_bins -1 2])
lambda_hat = N_spikes/duration_of_recording; % estimated firing rate
xlabel(['time in ' num2str(delta_t) 'sec'])
ylabel('spike')
title(['estimated firing rate = ',num2str(lambda_hat), ' spikes/sec'])


%% Multiple iterations of above experiment
%  to study the number of spikes in one experiment


clear

duration_of_recording = 1;        % [sec]
delta_t = 10^(-3);                % time resolution [sec}
r = 5;                            % firing rate [spikes/sec]

N_bins = duration_of_recording/delta_t;      % number of time intervals
prob_spike_per_bin = r*delta_t;

N_iterations = 5000;        % number of experiments

random_numbers = rand(N_iterations,N_bins);
spikes = random_numbers < prob_spike_per_bin;
N_spikes = sum(spikes,2);   % number of spikes in each experiment

figure(2)
range_spike_number = 0:2*max(N_spikes);    % range of spike number
empirical_distrib = hist(N_spikes,range_spike_number)/N_iterations;
bar(range_spike_number,empirical_distrib)
set(gca, 'fontsize',16)
xlabel('number of spikes')
ylabel('probability')
title('distribution of spike numbers')
lambda = r*duration_of_recording;       % parameter of Poisson distribution
poisson_distribution = poisspdf(range_spike_number,lambda);
binomial_distribution = binopdf(range_spike_number,N_bins,prob_spike_per_bin);
hold on
plot(range_spike_number,binomial_distribution,'gd','MarkerSize',10,'MarkerFaceColor','g')
plot(range_spike_number,poisson_distribution, 'r*','MarkerSize',20)
hold off
legend('empirical','binomial','Poisson')


%% Probability distribution of inter-spike intervals

clear

duration_of_recording = 10000;         % [sec]
delta_t = 10^(-3);                      % time resolution [sec}
r = 10;                                 % firing rate [spikes/sec]

N_bins = duration_of_recording/delta_t;      % number of time intervals
time_axis = 1:N_bins;
prob_spike_per_bin = r*delta_t;

random_numbers = rand(1,N_bins);
spikes = random_numbers < prob_spike_per_bin;
spike_times = find(spikes);         % spike times in delta_t units
spike_times = spike_times*delta_t;  % spike times in seconds
ISI = diff(spike_times);            % Inter-Spike Intervals
N_spikes = sum(spikes);             % number of spikes
N_ISI = N_spikes - 1;               % number of ISIs

figure(3)
range_ISI = 0:delta_t:max(ISI);         % range of ISI
empirical_distrib = hist(ISI,range_ISI)/(N_ISI*delta_t);
bar(range_ISI,empirical_distrib)
set(gca, 'fontsize',20)
xlim([0 max(ISI)])
hold on
beta = 1/r;                         % parameter of exponential distribution
exponential_distribution = exppdf(range_ISI,beta);
plot(range_ISI,exponential_distribution, 'r', 'linewidth',3)
hold off
legend('empirical','exponential')
xlabel(['ISI in sec'])
title('probability distribution of ISI')
