New exercise for point processes

pointprocesses/code/counthist.m

@@ -1,4 +1,4 @@
-function [ counts, bins ] = counthist( spikes, w )
+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
@@ -18,18 +18,17 @@ function [ counts, bins ] = counthist( spikes, w )
     end
     % histogram of spike counts:
     maxn = max( n );
-    [counts, bins ] = hist( n, 0:1:maxn+1 );
+    [counts, bins ] = hist( n, 0:1:maxn+10 );
     counts = counts / sum( counts );
     if nargout == 0
         bar( bins, counts );
         hold on;
         % Poisson distribution:
         rate = mean( r );
-        x = 0:1:20;
+        x = 0:1:maxn+10;
         l = rate*w;
         y = l.^x.*exp(-l)./factorial(x);
         plot( x, y, 'r', 'LineWidth', 3 );
-        xlim( [ 0 20 ] );
         hold off;
         xlabel( 'counts k' );
         ylabel( 'P(k)' );

pointprocesses/code/isihist.m

@@ -24,9 +24,9 @@ function isihist( isis, binwidth )
     misi = mean( isis );
     sdisi = std( isis );
     disi = sdisi^2.0/2.0/misi^3;
-    text( 0.5, 0.6, sprintf( 'mean=%.1f ms', 1000.0*misi ), 'Units', 'normalized' )
-    text( 0.5, 0.5, sprintf( 'std=%.1f ms', 1000.0*sdisi ), 'Units', 'normalized' )
-    text( 0.5, 0.4, sprintf( 'CV=%.2f', sdisi/misi ), 'Units', 'normalized' )
+    text( 0.95, 0.8, sprintf( 'mean=%.1f ms', 1000.0*misi ), 'Units', 'normalized', 'HorizontalAlignment', 'right' )
+    text( 0.95, 0.7, sprintf( 'std=%.1f ms', 1000.0*sdisi ), 'Units', 'normalized', 'HorizontalAlignment', 'right' )
+    text( 0.95, 0.6, sprintf( 'CV=%.2f', sdisi/misi ), 'Units', 'normalized', 'HorizontalAlignment', 'right' )
     %text( 0.5, 0.3, sprintf( 'D=%.1f Hz', disi ), 'Units', 'normalized' )
 end
 

pointprocesses/code/lifadaptspikes.m

@@ -3,8 +3,8 @@ function spikes = lifadaptspikes( trials, input, tmaxdt, D, tauadapt, adaptincr
 % with an adaptation current
 % trials: the number of trials to be generated
 % input: the stimulus either as a single value or as a vector
-% tmaxdt: in case of a single value stimulus the duration of a trial
-%         in case of a vector as a stimulus the time step
+% tmaxdt: in case of a single value stimulus: the duration of a trial
+%         in case of a vector as a stimulus: the time step
 % D: the strength of additive white noise
 % tauadapt: adaptation time constant
 % adaptincr: adaptation strength

pointprocesses/code/lifspikes.m

@@ -2,8 +2,8 @@ function spikes = lifspikes( trials, input, tmaxdt, D )
 % Generate spike times of a leaky integrate-and-fire neuron
 % trials: the number of trials to be generated
 % input: the stimulus either as a single value or as a vector
-% tmaxdt: in case of a single value stimulus the duration of a trial
-%         in case of a vector as a stimulus the time step
+% tmaxdt: in case of a single value stimulus: the duration of a trial
+%         in case of a vector as a stimulus: the time step
 % D: the strength of additive white noise
     
     tau = 0.01;

pointprocesses/code/lifspikesoustim.m

@@ -0,0 +1,20 @@
+function spikes = lifspikesoustim(trials, tmax, D, Iou, Dou, tauou )
+% Generate spike times of a leaky integrate-and-fire neuron with frozen
+% Ohrnstein-Uhlenbeck stimulus
+% trials: the number of trials to be generated
+% tmax: the duration of a trial
+% D: the strength of additive white noise
+% Iou: the mean input
+% Dou: noise strength of the frozen OU noise
+% tauou: time constant of the OU noise
+    
+    dt = 1e-4;
+    input = zeros(round(tmax/dt), 1);
+    n = 0.0;
+    noise = sqrt(2.0*Dou)*randn(length(input), 1)/sqrt(dt);
+    for i=1:length(noise)
+        n = n + ( - n + noise(i))*dt/tauou;
+        input(i) = Iou + n;
+    end
+    spikes = lifspikes(trials, input, dt, D );
+end

pointprocesses/code/plotcounthist.m

@@ -0,0 +1,24 @@
+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/plotisih.m

@@ -0,0 +1,19 @@
+maxisi = 300.0;
+subplot(1, 3, 1);
+poissonisis = isis(poissonspikes);
+isihist(poissonisis, 0.001);
+xlim([0, maxisi])
+title('Poisson');
+
+subplot(1, 3, 2);
+pifouisis = isis(pifouspikes);
+isihist(pifouisis, 0.001);
+xlim([0, maxisi])
+title('PIF OU');
+
+subplot(1, 3, 3);
+lifadaptisis = isis(lifadaptspikes);
+isihist(lifadaptisis, 0.001);
+xlim([0, maxisi])
+title('LIF adapt');
+savefigpdf(gcf, 'isihist.pdf', 20, 7);

pointprocesses/code/plotserialcorr.m

@@ -0,0 +1,17 @@
+maxlag = 10;
+rrange = [-0.5, 1.05];
+subplot(1, 3, 1);
+isiserialcorr(poissonisis, maxlag);
+ylim(rrange)
+title('Poisson');
+
+subplot(1, 3, 2);
+isiserialcorr(pifouisis, maxlag);
+ylim(rrange)
+title('PIF OU');
+
+subplot(1, 3, 3);
+isiserialcorr(lifadaptisis, maxlag);
+ylim(rrange)
+title('LIF adapt');
+savefigpdf(gcf, 'serialcorr.pdf', 20, 7);

pointprocesses/code/plotspikeraster.m

@@ -0,0 +1,13 @@
+subplot(1, 3, 1);
+spikeraster(poissonspikes, 1.0);
+title('Poisson');
+
+subplot(1, 3, 2);
+spikeraster(pifouspikes, 1.0);
+title('PIF OU');
+
+subplot(1, 3, 3);
+spikeraster(lifadaptspikes, 1.0);
+title('LIF adapt');
+
+savefigpdf(gcf, 'spikeraster.pdf', 15, 5);

pointprocesses/code/poissonspikes.m

@@ -12,9 +12,9 @@ function spikes = poissonspikes( trials, rate, tmax )
         p = 0.1
         dt = p/rate;
     end
-    spikes = cell( trials, 1 );
+    spikes = cell(trials, 1);
     for k=1:trials
-        x = rand( 1, round(tmax/dt) ); % uniform random numbers for each bin
-        spikes{k} = find( x < p ) * dt;
+        x = rand(round(tmax/dt), 1); % uniform random numbers for each bin
+        spikes{k} = find(x < p) * dt;
     end
 end

pointprocesses/code/psth.m

@@ -0,0 +1,14 @@
+function p = psth(spikes, dt, tmax)
+% plots a PSTH of the spikes with binwidth dt
+    t = 0.0:dt:tmax+dt;
+    p = zeros(1, length(t));
+    for k=1:length(spikes)
+        times = spikes{k};
+        [h, b] = hist(times, t);
+        p = p + h;
+    end
+    p = p/length(spikes)/dt;
+    t(end) = [];
+    p(end) = [];
+    plot(t, p);
+end

pointprocesses/code/spikeraster.m

@@ -1,15 +1,24 @@
-function spikeraster( spikes )
+function spikeraster(spikes, tmax)
 % Display a spike raster of the spike times given in spikes.
 % spikes: a cell array of vectors of spike times
+% tmax: plot spike raster upto tmax seconds
 
 ntrials = length(spikes);
 for k = 1:ntrials
-    times = 1000.0*spikes{k};  % conversion to ms
+    times = spikes{k};
+    times = times(times<tmax);
+    if tmax < 1.5
+        times = 1000.0*times; % conversion to ms
+    end
     for i = 1:length( times )
       line([times(i) times(i)],[k-0.4 k+0.4], 'Color', 'k' ); 
     end
 end
-xlabel( 'Time [ms]' );
+if tmax < 1.5
+    xlabel( 'Time [ms]' );
+else
+    xlabel( 'Time [s]' );
+end
 ylabel( 'Trials');
 ylim( [ 0.3 ntrials+0.7 ] )
 

pointprocesses/code/spikerate.m

@@ -0,0 +1,12 @@
+function r = spikerate(spikes, duration)
+% returns the average spike rate of the spikes
+% for the first duration seconds
+% spikes: a cell array of vectors of spike times
+
+    rates = zeros(length(spikes),1);
+    for k = 1:length(spikes)
+        times = spikes{k};
+        rates(k) = sum(times<duration)/duration;
+    end
+    r = mean(rates);
+end