Updated pointprocess exercises and code

pointprocesses/code/counthist.m

@@ -1,24 +1,36 @@
 function [counts, bins] = counthist(spikes, w)
 % computes count histogram and compare them with Poisson distribution
-% spikes: a cell array of vectors of spike times
+%
-% w: observation window duration for computing the counts
+% [counts, bins] = counthist(spikes, w)
+%   spikes: a cell array of vectors of spike times in seconds
+%   w: observation window duration in seconds for computing the counts
+%   counts: the histogram of counts normalized to probabilities
+%   bins: the bin centers for the histogram
 
     % collect spike counts:
     tmax = spikes{1}(end);
     n = [];
     r = [];
     for k = 1:length(spikes)
-        for tk = 0:w:tmax-w
+        times = spikes{k};
-            nn = sum( ( spikes{k} >= tk ) & ( spikes{k} < tk+w ) );
+%       alternative 1: count the number of spikes in each window:
-            %nn = length( find( ( spikes{k} >= tk ) & ( spikes{k} < tk+w ) ) );
+%         for tk = 0:w:tmax-w
+%             nn = sum( ( times >= tk ) & ( times < tk+w ) );
+%             %nn = length( find( ( times >= tk ) & ( times < tk+w ) ) );
+%             n = [ n nn ];
+%         end
+%     alternative 2: use the hist function to do that!
+        tbins = 0.5*w:w:tmax-0.5*w;
+        nn = hist(times, tbins);
         n = [ n nn ];
-        end
+        % the rate of the spikes:
-        rate = (length(spikes{k})-1)/(spikes{k}(end) - spikes{k}(1));
+        rate = (length(times)-1)/(times(end) - times(1));
         r = [ r rate ];
     end
     % histogram of spike counts:
     maxn = max( n );
     [counts, bins ] = hist( n, 0:1:maxn+10 );
+    % normalize to probabilities:
     counts = counts / sum( counts );
     if nargout == 0
         bar( bins, counts );
@@ -26,12 +38,11 @@ function [counts, bins] = counthist(spikes, w)
         % Poisson distribution:
         rate = mean( r );
         x = 0:1:maxn+10;
-        l = rate*w;
+        a = rate*w;
-        y = l.^x.*exp(-l)./factorial(x);
+        y = a.^x.*exp(-a)./factorial(x);
         plot( x, y, 'r', 'LineWidth', 3 );
         hold off;
         xlabel( 'counts k' );
         ylabel( 'P(k)' );
     end
 end
-

pointprocesses/code/isihist.m

@@ -1,7 +1,9 @@
 function isihist( isis, binwidth )
-% plot histogram of isis
+% plot histogram of interspike intervals
-% isis: vector of interspike intervals
+%
-% binwidth: optional width to be used for the isi bins
+% isihist(isis, binwidth)
+%   isis: vector of interspike intervals in seconds
+%   binwidth: optional width in seconds to be used for the isi bins
 
     if nargin < 2
         nperbin = 200;  % average number of data points per bin

pointprocesses/code/isis.m

@@ -1,10 +1,15 @@
 function isivec = isis( spikes )
 % returns a single list of isis computed from all trials in spikes
-% spikes: a cell array of vectors of spike times
+%
+% isivec = isis( spikes )
+%   spikes: a cell array of vectors of spike times in seconds
+%   isivec: a column vector with all the interspike intervalls
 
     isivec = [];
     for k = 1:length(spikes)
         difftimes = diff( spikes{k} );
+        % difftimes(:) ensures a column vector
+        % regardless of the type of vector spikes{k}
         isivec = [ isivec; difftimes(:) ];
     end
 end

pointprocesses/code/isiserialcorr.m

@@ -1,13 +1,18 @@
 function isicorr = isiserialcorr( isis, maxlag )
-% serial correlation of isis
+% serial correlation of interspike intervals
-% isis: vector of interspike intervals
+%
-% maxlag: the maximum lag
+% isicorr = isiserialcorr( isis, maxlag )
+%   isis: vector of interspike intervals in seconds
+%   maxlag: the maximum lag in seconds
+%   isicorr: vector with the serial correlations for lag 0 to maxlag
 
     lags = 0:maxlag;
     isicorr = zeros( size( lags ) );
     for k = 1:length(lags)
         lag = lags(k);
-        if length( isis ) > lag+10
+        if length( isis ) > lag+10        % ensure "enough" data
+	    % DANGER: the arguments to corr must be column vectors!
+            % We insure this in the isis() function that generats the isis.
             isicorr(k) = corr( isis(1:end-lag), isis(lag+1:end) );
         end
     end

pointprocesses/code/plotspikestats.m

@@ -1,5 +1,4 @@
 %% load data:
-
 clear all
 % alternative 1:
 % pro: no structs. contra: global unknown variables
@@ -22,35 +21,37 @@ x = load( 'lifadapt.mat' );
 lifadaptspikes = x.spikes;
 
 %% spike raster plots:
+tmax = 1.0;
 subplot(1, 3, 1);
-spikeraster(poissonspikes, 1.0);
+spikeraster(poissonspikes, tmax);
 title('Poisson');
 
 subplot(1, 3, 2);
-spikeraster(pifouspikes, 1.0);
+spikeraster(pifouspikes, tmax);
 title('PIF OU');
 
 subplot(1, 3, 3);
-spikeraster(lifadaptspikes, 1.0);
+spikeraster(lifadaptspikes, tmax);
 title('LIF adapt');
 
 %% isi histograms:
 maxisi = 300.0;
+binwidth = 0.002;
 subplot(1, 3, 1);
 poissonisis = isis(poissonspikes);
-isihist(poissonisis, 0.001);
+isihist(poissonisis, binwidth);
 xlim([0, maxisi])
 title('Poisson');
 
 subplot(1, 3, 2);
 pifouisis = isis(pifouspikes);
-isihist(pifouisis, 0.001);
+isihist(pifouisis, binwidth);
 xlim([0, maxisi])
 title('PIF OU');
 
 subplot(1, 3, 3);
 lifadaptisis = isis(lifadaptspikes);
-isihist(lifadaptisis, 0.001);
+isihist(lifadaptisis, binwidth);
 xlim([0, maxisi])
 title('LIF adapt');
 
@@ -71,3 +72,29 @@ subplot(1, 3, 3);
 isiserialcorr(lifadaptisis, maxlag);
 ylim(rrange)
 title('LIF adapt');
+
+%% spike counts:
+w = 0.1;
+cmax = 8;
+pmax = 0.5;
+subplot(1, 3, 1);
+counthist(poissonspikes, w);
+xlim([0 cmax])
+set(gca, 'XTick', 0:2:cmax)
+ylim([0 pmax])
+title('Poisson');
+
+subplot(1, 3, 2);
+counthist(pifouspikes, w);
+xlim([0 cmax])
+set(gca, 'XTick', 0:2:cmax)
+ylim([0 pmax])
+title('PIF OU');
+
+subplot(1, 3, 3);
+counthist(lifadaptspikes, w);
+xlim([0 cmax])
+set(gca, 'XTick', 0:2:cmax)
+ylim([0 pmax])
+title('LIF adapt');
+savefigpdf(gcf, 'counthist.pdf', 20, 7);

pointprocesses/code/poissonspikes.m

@@ -1,9 +1,11 @@
 function spikes = poissonspikes( trials, rate, tmax )
 % Generate spike times of a homogeneous poisson process
+%
+% spikes = poissonspikes( trials, rate, tmax )
 %   trials: number of trials that should be generated
 %   rate: the rate of the Poisson process in Hertz
 %   tmax: the duration of each trial in seconds
-% returns a cell array of vectors of spike times
+%   spikes: a cell array of vectors of spike times in seconds
 
     dt = 3.33e-5;
     p = rate*dt; % probability of event per bin of width dt

pointprocesses/code/spikeraster.m

@@ -1,6 +1,8 @@
 function spikeraster(spikes, tmax)
 % Display a spike raster of the spike times given in spikes.
-% spikes: a cell array of vectors of spike times
+%
+% spikeraster(spikes, tmax)
+%   spikes: a cell array of vectors of spike times in seconds
 %   tmax: plot spike raster upto tmax seconds
 
 ntrials = length(spikes);
@@ -16,8 +18,10 @@ for k = 1:ntrials
 end
 if tmax < 1.5
     xlabel( 'Time [ms]' );
+    xlim([0.0 1000.0*tmax]);
 else
     xlabel( 'Time [s]' );
+    xlim([0.0 tmax]);
 end
 ylabel( 'Trials');
 ylim( [ 0.3 ntrials+0.7 ] )