[pointprocesses] moved exercise codes to exercise folder

pointprocesses/code/counts.m

@@ -0,0 +1,25 @@
+function n = counts(spikes, w)
+% Count spikes in time windows.
+%
+% Arguments:
+%   spikes: a cell array of vectors of spike times in seconds
+%   w: duration of window in seconds for computing the counts
+%
+% Returns:
+%   n: vector with spike counts
+    tmax = spikes{1}(end);
+    n = [];
+    for k = 1:length(spikes)
+        times = spikes{k};
+%       alternative 1: count the number of spikes in each window:
+%         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
+end

pointprocesses/code/plotcounthist.m

@@ -1,4 +1,4 @@
-function counthist(spikes, w)
+function plotcounthist(spikes, w)
 % Plot histogram of spike counts.
 %
 % Arguments:

pointprocesses/code/spikecounts.m

@@ -1,30 +0,0 @@
-function counts = spikecounts(spikes, w)
-% Compute vector of spike counts.
-%
-% counts = spikecounts(spikes, w)
-%
-% Arguments:
-%   spikes: a cell array of vectors of spike times in seconds
-%   w: observation window duration in seconds for computing the counts
-%
-% Returns:
-%   counts: vector of spike counts
-
-    % collect spike counts:
-    tmax = spikes{1}(end);
-    counts = [];
-    for k = 1:length(spikes)
-        times = spikes{k};
-%       method 1: count the number of spikes in each window:
-%         for tk = 0:w:tmax-w
-%             nn = sum((times >= tk) & (times < tk+w));
-%             %nn = length(times((times >= tk) & (times < tk+w)));
-%             %nn = length(find((times >= tk) & (times < tk+w)));
-%             counts = [counts nn];
-%         end
-%     method 2: use the hist() function to do that!
-        tbins = 0.5*w:w:tmax-0.5*w;
-        nn = hist(times, tbins);
-        counts = [counts nn];
-    end
-end

pointprocesses/code/spikecountshists.m

@@ -1,17 +0,0 @@
-w = 0.1;
-bins = 0.0:1.0:10.0;
-
-counts{1} = spikecounts(poissonspikes, w);
-counts{2} = spikecounts(pifouspikes, w);
-counts{3} = spikecounts(lifadaptspikes, w);
-titles = {'Poisson', 'PIF OU', 'LIF adapt'};
-for k = 1:3
-    subplot(1, 3, k);
-    [h, b] = hist(counts{k}, bins);
-    bar(b, h/sum(h));
-    title(titles{k})
-    xlabel('Spike count');
-    xlim([-0.5, 8.5]);
-    ylim([0.0 0.8]);
-end
-savefigpdf(gcf, 'spikecountshists.pdf', 20, 7);

pointprocesses/exercises/counts.m

@@ -0,0 +1,25 @@
+function n = counts(spikes, w)
+% Count spikes in time windows.
+%
+% Arguments:
+%   spikes: a cell array of vectors of spike times in seconds
+%   w: duration of window in seconds for computing the counts
+%
+% Returns:
+%   n: vector with spike counts
+    tmax = spikes{1}(end);
+    n = [];
+    for k = 1:length(spikes)
+        times = spikes{k};
+%       alternative 1: count the number of spikes in each window:
+%         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
+end

pointprocesses/exercises/fanoplot.m

@@ -0,0 +1,32 @@
+function fanoplot(spikes, titles)
+% Plots fano factor as a function of window size.
+%
+% Arguments:
+%   spikes: a cell array of vectors of spike times
+%   titles: string that is used as a title for the plots
+    windows = logspace(-3.0, -0.5, 100);
+    mc = zeros(1, length(windows));
+    vc = zeros(1, length(windows));
+    for j = 1:length(windows)
+        w = windows(j);
+        spikecounts = counts(spikes, w);
+        mc(j) = mean(spikecounts);
+        vc(j) = var(spikecounts);
+    end
+    
+    subplot(1, 2, 1);
+    scatter(mc, vc, 'filled');
+    title(titles);
+    xlabel('Mean count');
+    ylabel('Count variance');
+
+    subplot(1, 2, 2);
+    scatter(1000.0*windows, vc ./ mc, 'filled');
+    title(titles);
+    xlabel('Window [ms]');
+    ylabel('Fano factor');
+    xlim(1000.0*[windows(1) windows(end)])
+    ylim([0.0 1.1]);
+    set(gca, 'XScale', 'log');
+end
+

pointprocesses/exercises/fanoplots.m

@@ -0,0 +1,9 @@
+spikes{1} = poissonspikes;
+spikes{2} = pifouspikes;
+spikes{3} = lifadaptspikes;
+idents = {'poisson', 'pifou', 'lifadapt'};
+for k = 1:3
+    figure(k)
+    fanoplot(spikes{k}, titles{k});
+    savefigpdf(gcf, sprintf('fanoplots%s.pdf', idents{k}), 20, 7);
+end

pointprocesses/exercises/isihist.m

@@ -0,0 +1,29 @@
+function [pdf, centers] = isihist(isis, binwidth)
+% Compute normalized histogram of interspike intervals.
+%
+% [pdf, centers] = isihist(isis, binwidth)
+%
+% Arguments:
+%   isis: vector of interspike intervals in seconds
+%   binwidth: optional width in seconds to be used for the isi bins
+%
+% Returns:
+%   pdf: vector with pdf of interspike intervals in Hertz
+%   centers: vector with centers of interspikeintervalls in seconds
+
+    if nargin < 2
+        % compute good binwidth:
+        nperbin = 200;     % average number of data points per bin
+        bins = length(isis)/nperbin;   % number of bins
+        binwidth = max(isis)/bins;
+        if binwidth < 5e-4             % half a millisecond
+            binwidth = 5e-4;
+        end
+    end
+    bins = 0.5*binwidth:binwidth:max(isis);
+    % histogram data:
+    [nelements, centers] = hist(isis, bins);
+    % normalization (integral = 1):
+    pdf = nelements / sum(nelements) / binwidth;
+end
+

pointprocesses/exercises/isis.m

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

pointprocesses/exercises/isiserialcorr.m

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

pointprocesses/exercises/plotisihist.m

@@ -0,0 +1,31 @@
+function plotisihist(isis, binwidth)
+% Plot and annotate histogram of interspike intervals.
+%
+% plotisihist(isis, binwidth)
+%
+% Arguments:
+%   isis: vector of interspike intervals in seconds
+%   binwidth: optional width in seconds to be used for the isi bins
+
+    % compute normalized histogram:
+    if nargin < 2
+        [pdf, centers] = isihist(isis);
+    else
+        [pdf, centers] = isihist(isis, binwidth);
+    end
+
+    % plot:
+    bar(1000.0*centers, pdf);     % milliseconds on x-axis
+    xlabel('ISI [ms]')
+    ylabel('p(ISI) [1/s]')
+    % annotation:
+    misi = mean(isis);
+    sdisi = std(isis);
+    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')
+end
+

pointprocesses/exercises/plotisiserialcorr.m

@@ -0,0 +1,14 @@
+function isicorr = plotisiserialcorr(isivec, maxlag)
+% plot serial correlation of interspike intervals
+%
+% Arguments:
+%   isivec: vector of interspike intervals in seconds
+%   maxlag: the maximum lag
+    [isicorr, lags] = isiserialcorr(isivec, maxlag);
+    plot(lags, isicorr, '-b');
+    hold on;
+    scatter(lags, isicorr, 20.0, 'b', 'filled');
+    hold off;
+    xlabel('Lag k')
+    ylabel('\rho_k')
+end

pointprocesses/exercises/plotpoissonisih.m

@@ -0,0 +1,12 @@
+maxisi = 300.0;
+poissonisis = isis(poissonspikes);
+rate = 1.0/mean(poissonisis);
+plotisihist(poissonisis, 0.001);
+tt = linspace(0.0, 0.001*maxisi, 200);
+pexp = rate*exp(-tt*rate);
+hold on;
+plot(1000.0*tt, pexp, 'r')
+hold off;
+xlim([0, maxisi])
+title('Poisson');
+savefigpdf(gcf, 'poissonisihist.pdf', 10, 7);

pointprocesses/exercises/plotserialcorr.m

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

pointprocesses/exercises/pointprocesses-1.tex

@@ -15,6 +15,9 @@
 
 \begin{questions}
 
+  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
+  \question \qt{Read in chapter 10 on ``Spike-train analysis'' sections 10.1 -- 10.5!}\vspace{-3ex}
+
   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
   \question \qt{Statistics of interspike intervals}
   Download the files \code{poisson.mat},
@@ -61,11 +64,12 @@
     of the spike. When appropriate, the function should use milliseconds 
     for the time axis instead of seconds.
 
-    Use this function to plot the first second of the
+    Use this function to plot the first second of the spike rasters of
-    spike rasters of the three neurons.
+    the three neurons. Carefully look at the raster plots. How do the
+    spike rasters differ between the three neurons?
     \begin{solution}
-      \lstinputlisting{../code/rasterplot.m}
+      \lstinputlisting{rasterplot.m}
-      \lstinputlisting{../code/plotspikeraster.m}
+      \lstinputlisting{plotspikeraster.m}
       \mbox{}\\[-3ex]
       \colorbox{white}{\includegraphics[width=1\textwidth]{spikeraster}}
     \end{solution}
@@ -73,7 +77,7 @@
     \part Write a function that returns a single vector containing the
     interspike intervals of all trials of spike times.
     \begin{solution}
-      \lstinputlisting{../code/isis.m}
+      \lstinputlisting{isis.m}
     \end{solution}
     
     \part Write a function that computes an estimate of the
@@ -90,33 +94,51 @@
     the standard deviation and the coefficient of variation and
     display them in the plot as well.
 
-    Use this and the previous functions to compare the 
+    Use this and the previous functions to compare the interspike
-    interspike interval statistics of the three neurons.
+    interval statistics of the three neurons. How do the ISI
+    histograms differ between the three neurons?
     \begin{solution}
-      \lstinputlisting{../code/isihist.m}
+      \lstinputlisting{isihist.m}
-      \lstinputlisting{../code/plotisihist.m}
+      \lstinputlisting{plotisihist.m}
-      \lstinputlisting{../code/plotisihs.m}
+      \lstinputlisting{plotisihs.m}
       \mbox{}\\[-3ex]
       \colorbox{white}{\includegraphics[width=1\textwidth]{isihist}}
     \end{solution}
+
+    \part Compare the interspike interval histogram of the Poisson
+    spike train to the expected exponential distribution
+    \[ p(T) = \lambda e^{-\lambda T} \; .\] of a Poisson process with
+    firing rate $\lambda$. Estimate the firing rate from the inverse
+    of the mean interspike interval.
+    \begin{solution}
+      \lstinputlisting{plotpoissonisih.m}
+      %\mbox{}\\[-3ex]
+      %\colorbox{white}{\includegraphics[width=1\textwidth]{poissonisihist}}
+    \end{solution}
   
     % XXX Add return map!!! XXX
-    \part Write a function that computes and plots the serial
+    \part Write a function that computes the serial correlations of
-    correlations of interspike intervals for lags upto
+    interspike intervals for lags up to \code{maxlag}. The serial
-    \code{maxlag}. The serial correlations $\rho_k$ for lag $k$ of the
+    correlations $\rho_k$ for lag $k$ of the interspike intervals
-    interspike intervals $T_i$ are the correlation coefficients
+    $T_i$ are the correlation coefficients between interspike
-    between interspike intervals $T_i$ and the intervals $T_{i+k}$
+    intervals $T_i$ and the intervals $T_{i+k}$ that are shifted by
-    that are shifted by lag $k$:
+    lag $k$:
     \[ \rho_k = \frac{\langle (T_{i+k} - \langle T \rangle)(T_i -
       \langle T \rangle) \rangle}{\langle (T_i - \langle T
       \rangle)^2\rangle} = \frac{{\rm cov}(T_{i+k}, T_i)}{{\rm
-        var}(T_i)} = {\rm corr}(T_{i+k}, T_i) \] 
+        var}(T_i)} = {\rm corr}(T_{i+k}, T_i) \]
+
+    Write another function that plots the serial correlations.
 
-    Use this function to compare the serial correlations of the
+    Use the functions to compare the serial correlations of the
-    interspike intervals of the three neurons.
+    interspike intervals of the three neurons. What do the serial
+    correlations mean? How do they relate to what you see in the
+    raster plots? Do the serial correlations of the Poisson spike
+    train match the expectations?
     \begin{solution}
-      \lstinputlisting{../code/isiserialcorr.m}
+      \lstinputlisting{isiserialcorr.m}
-      \lstinputlisting{../code/plotserialcorr.m}
+      \lstinputlisting{plotisiserialcorr.m}
+      \lstinputlisting{plotserialcorr.m}
       \colorbox{white}{\includegraphics[width=1\textwidth]{serialcorr}}
     \end{solution}
     
@@ -127,16 +149,22 @@
   \question \qt{Statistics of spike counts}
   Now let's have a look at the statistics of the spike counts.
   \begin{parts}
-    \part Write a function that counts and returns the number of
+    \part Write a function that counts and returns a vector with the
-    spikes in windows of a given width $W$.
+    number of spikes in windows of a given width $W$.
 
     Use this function to generate a properly normalized histogram of
     spike counts for the data of the three types of neurons. Use
     100\,ms for the window width.
+
+    Compare the distributions with the Poisson distribution expected
+    for a Poisson spike train:
+    \[ P(k) = \frac{(\lambda W)^ke^{\lambda W}}{k!} \; , \] where
+    $\lambda$ is the rate of the spike train that you should estimate
+    from the data.
     \begin{solution}
-      \lstinputlisting{../code/spikecounts.m}
+      \lstinputlisting{counts.m}
       \newsolutionpage
-      \lstinputlisting{../code/spikecountshists.m}
+      \lstinputlisting{spikecounthists.m}
       \colorbox{white}{\includegraphics[width=1\textwidth]{spikecountshists}}
     \end{solution}
 
@@ -148,9 +176,9 @@
     (the mean spike count increases for larger window widths). The
     other plot showing the Fano factor as a function of window width.
     \begin{solution}
-      \lstinputlisting{../code/fanoplot.m}
+      \lstinputlisting{fanoplot.m}
       \newsolutionpage
-      \lstinputlisting{../code/fanoplots.m}
+      \lstinputlisting{fanoplots.m}
       \colorbox{white}{\includegraphics[width=1\textwidth]{fanoplotspoisson}}
       \colorbox{white}{\includegraphics[width=1\textwidth]{fanoplotspifou}}
       \colorbox{white}{\includegraphics[width=1\textwidth]{fanoplotslifadapt}}

pointprocesses/exercises/pointprocesses-2.tex

@@ -103,7 +103,7 @@ f\"ur gen\"ugend kleine $\Delta t$.
   soll. Hinweis: es gibt eine \code{matlab} Funktion, die die
   Fakult\"at $n!$ berechnet.
   \begin{solution}
-    \lstinputlisting{../code/counthist.m}
+    \lstinputlisting{counthist.m}
     \colorbox{white}{\includegraphics[width=0.48\textwidth]{poissoncounthistdist100hz10ms}}
     \colorbox{white}{\includegraphics[width=0.48\textwidth]{poissoncounthistdist100hz100ms}}
   \end{solution}
@@ -115,7 +115,7 @@ f\"ur gen\"ugend kleine $\Delta t$.
   den Ergebnissen angefertigt werden: (i) Varianz gegen Mittelwert der counts.
   (ii) Fano Faktor als Funktion der Fensterbreite.
   \begin{solution}
-    \lstinputlisting{../code/fano.m}
+    \lstinputlisting{fanoplot.m}
     \colorbox{white}{\includegraphics[width=0.98\textwidth]{poissonfano100hz}}
   \end{solution}
   

pointprocesses/exercises/rasterplot.m

@@ -0,0 +1,43 @@
+function rasterplot(spikes, tmax)
+% Display a spike raster of the spike times given in spikes.
+%
+% rasterplot(spikes, tmax)
+%
+% Arguments:
+%   spikes: a cell array of vectors of spike times in seconds
+%   tmax: plot spike raster up to tmax seconds
+    in_msecs = tmax < 1.5
+    spiketimes = [];
+    trials = [];
+    ntrials = length(spikes);
+    for k = 1:ntrials
+        times = spikes{k};
+        times = times(times<tmax);
+        if in_msecs
+            times = 1000.0*times;               % conversion to ms
+        end
+        % (x,y) pairs for start and stop of stroke
+        % plus nan separating strokes:
+        spiketimes = [spiketimes, ...
+                       [times(:)'; times(:)'; times(:)'*nan]];
+        trials = [trials, ...
+                   [ones(1, length(times)) * (k-0.4); ...
+                    ones(1, length(times)) * (k+0.4); ...
+                    ones(1, length(times)) * nan]];
+    end
+    % convert matrices into simple vectors of (x,y) pairs:
+    spiketimes = spiketimes(:);
+    trials = trials(:);
+    % plotting this is lightning fast:
+    plot(spiketimes, trials, 'k')
+    if in_msecs
+        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]);
+end
+

pointprocesses/exercises/spikecountshists.m

@@ -0,0 +1,22 @@
+w = 0.1;
+bins = 0:1:10;
+
+spikecounts{1} = counts(poissonspikes, w);
+spikecounts{2} = counts(pifouspikes, w);
+spikecounts{3} = counts(lifadaptspikes, w);
+titles = {'Poisson', 'PIF OU', 'LIF adapt'};
+for k = 1:3
+    subplot(1, 3, k);
+    [h, b] = hist(spikecounts{k}, bins);
+    bar(b, h/sum(h));
+    mu = mean(spikecounts{k});
+    ppois = poisspdf(bins, mu);
+    hold on;
+    scatter(bins, ppois, 'filled');
+    hold off;
+    title(titles{k})
+    xlabel('Spike count');
+    xlim([-0.5, 8.5]);
+    ylim([0.0 0.8]);
+end
+savefigpdf(gcf, 'spikecountshists.pdf', 20, 7);