From 7ec4cddd628f21ce9e49ad06368881a61b4bd2f3 Mon Sep 17 00:00:00 2001
From: Jan Benda <jan.benda@uni-tuebingen.de>
Date: Mon, 18 Jan 2021 23:03:30 +0100
Subject: [PATCH] [pointprocesses] improved spike count and fano factor
 exercises

---
 pointprocesses/code/counthist.m           | 37 ---------
 pointprocesses/code/fano.m                | 39 ----------
 pointprocesses/code/fanoplot.m            | 33 --------
 pointprocesses/code/fanoplots.m           |  9 ---
 pointprocesses/code/hompoissonspikes.m    |  3 -
 pointprocesses/code/isihist.m             | 14 +---
 pointprocesses/code/isis.m                |  4 +-
 pointprocesses/code/isiserialcorr.m       |  5 +-
 pointprocesses/code/plotcounthist.m       | 38 ++++-----
 pointprocesses/code/plotfanofactor.m      | 31 ++++++++
 pointprocesses/code/plotisihist.m         | 15 +---
 pointprocesses/code/plotisiserialcorr.m   |  2 -
 pointprocesses/code/rasterplot.m          |  4 +-
 pointprocesses/lecture/pointprocesses.tex | 95 ++++++++++++++---------
 14 files changed, 109 insertions(+), 220 deletions(-)
 delete mode 100644 pointprocesses/code/counthist.m
 delete mode 100644 pointprocesses/code/fano.m
 delete mode 100644 pointprocesses/code/fanoplot.m
 delete mode 100644 pointprocesses/code/fanoplots.m
 create mode 100644 pointprocesses/code/plotfanofactor.m

diff --git a/pointprocesses/code/counthist.m b/pointprocesses/code/counthist.m
deleted file mode 100644
index 2f5b746..0000000
--- a/pointprocesses/code/counthist.m
+++ /dev/null
@@ -1,37 +0,0 @@
-function counthist(spikes, w)
-% Plot histogram of spike counts.
-%
-% counthist(spikes, w)
-%
-% Arguments:
-%   spikes: a cell array of vectors of spike times in seconds
-%   w: duration of window in seconds for computing the counts
-    % collect spike counts:
-    tmax = spikes{1}(end);
-    n = [];
-    r = [];
-    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
-
-    % histogram of spike counts:
-    maxn = max(n);
-    [counts, bins] = hist(n, 0:1:maxn+10);
-    % normalize to probabilities:
-    counts = counts / sum(counts);
-
-    % plot:
-    bar(bins, counts);
-    xlabel('counts k');
-    ylabel('P(k)');
-end
diff --git a/pointprocesses/code/fano.m b/pointprocesses/code/fano.m
deleted file mode 100644
index a73cada..0000000
--- a/pointprocesses/code/fano.m
+++ /dev/null
@@ -1,39 +0,0 @@
-function fano( spikes )
-% computes fano factor as a function of window size
-% spikes: a cell array of vectors of spike times
-
-    tmax = spikes{1}(end);
-    windows = 0.01:0.05:0.01*tmax;
-    mc = windows;
-    vc = windows;
-    ff = windows;
-    fs = windows;
-    for j = 1:length(windows)
-        w = windows( j );
-        % collect counts:
-        n = [];
-        for k = 1:length(spikes)
-            for tk = 0:w:tmax-w
-                nn = sum( ( spikes{k} >= tk ) & ( spikes{k} < tk+w ) );
-                %nn = length( find( ( spikes{k} >= tk ) & ( spikes{k} < tk+w ) ) );
-                n = [ n nn ];
-            end
-        end
-        % statistics for current window:
-        mc(j) = mean( n );
-        vc(j) = var( n );
-        ff(j) = vc( j )/mc( j );
-        fs(j) = sqrt(vc( j )/mc( j ));
-    end
-    
-    subplot( 1, 2, 1 );
-    scatter( mc, vc, 'filled' );
-    xlabel( 'Mean count' );
-    ylabel( 'Count variance' );
-
-    subplot( 1, 2, 2 );
-    scatter( 1000.0*windows, fs, 'filled' );
-    xlabel( 'Window W [ms]' );
-    ylabel( 'Fano factor' );
-end
-
diff --git a/pointprocesses/code/fanoplot.m b/pointprocesses/code/fanoplot.m
deleted file mode 100644
index b53c42f..0000000
--- a/pointprocesses/code/fanoplot.m
+++ /dev/null
@@ -1,33 +0,0 @@
-function fanoplot(spikes, titles)
-% computes and plots fano factor as a function of window size
-% 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 = windows;
-    vc = windows;
-    ff = windows;
-    for j = 1:length(windows)
-        w = windows(j);
-        counts = spikecounts(spikes, w);
-        % statistics for current window:
-        mc(j) = mean(counts);
-        vc(j) = var(counts);
-        ff(j) = vc(j)/mc(j);
-    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, ff, '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
-
diff --git a/pointprocesses/code/fanoplots.m b/pointprocesses/code/fanoplots.m
deleted file mode 100644
index 949eb91..0000000
--- a/pointprocesses/code/fanoplots.m
+++ /dev/null
@@ -1,9 +0,0 @@
-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
diff --git a/pointprocesses/code/hompoissonspikes.m b/pointprocesses/code/hompoissonspikes.m
index e5caf62..9fe8446 100644
--- a/pointprocesses/code/hompoissonspikes.m
+++ b/pointprocesses/code/hompoissonspikes.m
@@ -2,8 +2,6 @@ function spikes = hompoissonspikes(rate, trials, tmax)
 % Generate spike times of a homogeneous poisson process
 % using the exponential interspike interval distribution.
 %
-% spikes = hompoissonspikes(rate, trials, tmax)
-%
 % Arguments:
 %   rate: the rate of the Poisson process in Hertz
 %   trials: number of trials that should be generated
@@ -11,7 +9,6 @@ function spikes = hompoissonspikes(rate, trials, tmax)
 %
 % Returns:
 %   spikes: a cell array of vectors of spike times in seconds
-
     spikes = cell(trials, 1);
     mu = 1.0/rate;
     nintervals = 2*round(tmax/mu);
diff --git a/pointprocesses/code/isihist.m b/pointprocesses/code/isihist.m
index a24989c..6c89b38 100644
--- a/pointprocesses/code/isihist.m
+++ b/pointprocesses/code/isihist.m
@@ -1,25 +1,13 @@
 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
+%   binwidth: 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);
diff --git a/pointprocesses/code/isis.m b/pointprocesses/code/isis.m
index 509dcc5..d5eb7d6 100644
--- a/pointprocesses/code/isis.m
+++ b/pointprocesses/code/isis.m
@@ -1,13 +1,11 @@
 function isivec = isis(spikes)
 % returns a single list of isis computed from all trials in spikes
 %
-% isivec = isis(spikes)
-%
 % Arguments:
 %   spikes: a cell array of vectors of spike times in seconds
 %
 % Returns:
-%   isivec: a column vector with all the interspike intervalls
+%   isivec: a column vector with all the interspike intervals
     isivec = [];
     for k = 1:length(spikes)
         difftimes = diff(spikes{k});
diff --git a/pointprocesses/code/isiserialcorr.m b/pointprocesses/code/isiserialcorr.m
index 1b813fc..fc5ec17 100644
--- a/pointprocesses/code/isiserialcorr.m
+++ b/pointprocesses/code/isiserialcorr.m
@@ -1,8 +1,6 @@
 function [isicorr, lags] = isiserialcorr(isivec, maxlag)
 % serial correlation of interspike intervals
 %
-% isicorr = isiserialcorr(isivec, maxlag)
-%
 % Arguments:
 %   isivec: vector of interspike intervals in seconds
 %   maxlag: the maximum lag
@@ -16,8 +14,7 @@ function [isicorr, lags] = isiserialcorr(isivec, maxlag)
         lag = lags(k);
         if length(isivec) > lag+10        % ensure "enough" data
 	    % NOTE: the arguments to corr must be column vectors!
-            % We insure this in the isis() function that 
-            % generates the isivec.
+            % We insure this already in the isis() function.
             isicorr(k) = corr(isivec(1:end-lag), isivec(lag+1:end));
         end
     end
diff --git a/pointprocesses/code/plotcounthist.m b/pointprocesses/code/plotcounthist.m
index 7a0a51d..c487e9a 100644
--- a/pointprocesses/code/plotcounthist.m
+++ b/pointprocesses/code/plotcounthist.m
@@ -1,24 +1,14 @@
-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);
\ No newline at end of file
+function counthist(spikes, w)
+% Plot histogram of spike counts.
+%
+% Arguments:
+%   spikes: a cell array of vectors of spike times in seconds
+%   w: duration of window in seconds for computing the counts
+    n = counts(spikes, w);
+    maxn = max(n);
+    [counts, bins] = hist(n, 0:1:maxn+10);
+    counts = counts / sum(counts);
+    bar(bins, counts);
+    xlabel('counts k');
+    ylabel('P(k)');
+end
diff --git a/pointprocesses/code/plotfanofactor.m b/pointprocesses/code/plotfanofactor.m
new file mode 100644
index 0000000..7335bdf
--- /dev/null
+++ b/pointprocesses/code/plotfanofactor.m
@@ -0,0 +1,31 @@
+function plotfanofactor(spikes, wmin, wmax)
+% Compute and plot Fano factor as a function of window size.
+%
+% Arguments:
+%   spikes: a cell array of vectors of spike times in seconds
+%   wmin: minimum window size in seconds
+%   wmax: maximum window size in seconds
+    windows = logspace(log10(wmin), log10(wmax), 100);
+    mc = zeros(1, length(windows));
+    vc = zeros(1, length(windows));
+    for k = 1:length(windows)
+        w = windows(k);
+        n = counts(spikes, w);
+        mc(k) = mean(n);
+        vc(k) = var(n);
+    end
+    
+    subplot(1, 2, 1);
+    scatter(mc, vc, 'filled');
+    xlabel('Mean count');
+    ylabel('Count variance');
+
+    subplot(1, 2, 2);
+    scatter(1000.0*windows, vc ./ mc, 'filled');
+    xlabel('Window [ms]');
+    ylabel('Fano factor');
+    xlim(1000.0*[windows(1) windows(end)])
+    ylim([0.0 1.1]);
+    set(gca, 'XScale', 'log');
+end
+
diff --git a/pointprocesses/code/plotisihist.m b/pointprocesses/code/plotisihist.m
index 2f3cf93..bfdf179 100644
--- a/pointprocesses/code/plotisihist.m
+++ b/pointprocesses/code/plotisihist.m
@@ -1,24 +1,13 @@
 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:
+%   binwidth: width in seconds to be used for the ISI bins
+    [pdf, centers] = isihist(isis, binwidth);
     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), ...
diff --git a/pointprocesses/code/plotisiserialcorr.m b/pointprocesses/code/plotisiserialcorr.m
index 3f21418..ce806e3 100644
--- a/pointprocesses/code/plotisiserialcorr.m
+++ b/pointprocesses/code/plotisiserialcorr.m
@@ -1,8 +1,6 @@
 function isicorr = plotisiserialcorr(isivec, maxlag)
 % plot serial correlation of interspike intervals
 %
-% plotisiserialcorr(isivec, maxlag)
-%
 % Arguments:
 %   isivec: vector of interspike intervals in seconds
 %   maxlag: the maximum lag
diff --git a/pointprocesses/code/rasterplot.m b/pointprocesses/code/rasterplot.m
index beb1970..7551552 100644
--- a/pointprocesses/code/rasterplot.m
+++ b/pointprocesses/code/rasterplot.m
@@ -1,8 +1,6 @@
 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
@@ -21,7 +19,7 @@ function rasterplot(spikes, tmax)
                     ones(1, length(times)) * (k+0.4); ...
                     ones(1, length(times)) * nan]];
     end
-    % convert matrices into simple vectors of (x,y) pairs:
+    % convert matrices into column vectors of (x,y) pairs:
     spiketimes = spiketimes(:);
     trials = trials(:);
     % plotting this is lightning fast:
diff --git a/pointprocesses/lecture/pointprocesses.tex b/pointprocesses/lecture/pointprocesses.tex
index 8a297ee..f6796d2 100644
--- a/pointprocesses/lecture/pointprocesses.tex
+++ b/pointprocesses/lecture/pointprocesses.tex
@@ -239,19 +239,18 @@ describing univariate data sets of real numbers:
 
 \begin{exercise}{isihist.m}{}
   Implement a function \varcode{isihist()} that calculates the
-  normalized interspike interval histogram. The function should take
-  two input arguments; (i) a vector of interspike intervals and (ii)
-  the width of the bins used for the histogram. It returns the
+  normalized interspike interval histogram. The function should take a
+  vector of interspike intervals and the width of the bins to be used
+  for the histogram as input arguments. The function returns the
   probability density as well as the centers of the bins.
 \end{exercise}
 
 \begin{exercise}{plotisihist.m}{}
-  Implement a function that takes the returned values of
-  \varcode{isihist()} as input arguments and then plots the data. The
-  plot should show the histogram with the x-axis scaled to
-  milliseconds and should be annotated with the average ISI, the
-  standard deviation, and the coefficient of variation of the ISIs
-  (\figref{isihexamplesfig}).
+  Implement a function that uses the \varcode{isihist()} function from
+  the previous exercise to plot an ISI histogram. The plot shows the
+  histogram with the x-axis scaled to milliseconds, annotated with the
+  average ISI, the standard deviation, and the coefficient of
+  variation of the ISIs (\figref{isihexamplesfig}).
 \end{exercise}
 
 \subsection{Interval correlations}
@@ -332,8 +331,8 @@ characterized by the non-zero serial correlations
 
 \begin{exercise}{isiserialcorr.m}{}
   Implement a function \varcode{isiserialcorr()} that takes a vector
-  of interspike intervals as input argument and calculates the serial
-  correlations up to some maximum lag.
+  of interspike intervals as input and computes serial correlations up
+  to some maximum lag.
 \end{exercise}
 
 \begin{exercise}{plotisiserialcorr.m}{}
@@ -398,27 +397,10 @@ Because spike counts are unitless and positive numbers the
   \end{equation}
   is a commonly used measure for quantifying the variability of event
   counts relative to the mean number of events. In particular for
-  homogeneous Poisson processes the Fano factor equals one,
-  independently of the time window $W$.
+  homogeneous Poisson processes the Fano factor equals exactly one and
+  is independent of the time window $W$.
 \end{itemize}
 
-Note that all of these statistics depend in general on the chosen
-window length $W$. The average spike count, for example, grows
-linearly with $W$ for sufficiently large time windows: $\langle n
-\rangle = r W$, \eqnref{firingrate}. Doubling the counting window
-doubles the spike count. As does the spike-count variance
-(\figref{fanofig}). At smaller time windows the statistics of the
-event counts might depend on the particular duration of the counting
-window. There might be an optimal time window for which the variance
-of the spike count is minimal. The Fano factor plotted as a function
-of the time window illustrates such properties of point processes
-(\figref{fanofig}).
-
-This also has consequences for information transmission in neural
-systems. The lower the variance in spike count relative to the
-averaged count, the higher the signal-to-noise ratio at which
-information encoded in the mean spike count is transmitted.
-
 \begin{figure}[t]
   \includegraphics{fanoexamples}
   \titlecaption{\label{fanofig} Fano factor.}{Counting events in time
@@ -435,13 +417,52 @@ information encoded in the mean spike count is transmitted.
     interspike intervals (right).}
 \end{figure}
 
-\begin{exercise}{counthist.m}{}
-  Implement a function \varcode{counthist()} that calculates and plots
-  the distribution of spike counts observed in a certain time
-  window. The function should take two input arguments: a cell-array
-  of vectors containing the spike times in seconds observed in a
-  number of trials, and the duration of the time window that is used
-  to evaluate the counts.
+Note that all of these statistics depend in general on the chosen
+window length $W$ used for counting the events. The average spike
+count, for example, grows linearly with $W$ for sufficiently large
+time windows: $\langle n \rangle = r W$, \eqnref{firingrate}. Doubling
+the counting window doubles the spike count. As does the spike-count
+variance (\figref{fanofig}). At smaller time windows the statistics of
+the event counts might depend on the particular duration of the
+counting window. There could be an optimal time window for which the
+variance of the spike count is minimal. The Fano factor plotted as a
+function of the time window illustrates such properties of point
+processes in a single graph (\figref{fanofig}).
+
+This also has consequences for information transmission in neural
+systems. The membrane time constant of a post-synaptic neuron defines
+a counting window. If input spikes arrive within about one time
+constant they are added up. The lower the variance in spike count
+relative to the averaged count, the higher the signal-to-noise ratio
+at which information encoded in the mean spike count is
+transmitted. If the membrane time constant of the target neuron
+matches the counting window where the Fano factor is minimal, then
+information is potentially transmitted at the highest possible
+signal-to-noise ratio. For this reason, the Fano factor is used in the
+Neurosciences to quantify and analyze reliability of neuronal
+responses.
+
+\begin{exercise}{counts.m}{}
+  Write a function \varcode{counts()} that counts the number of spikes
+  in windows of given duration and returns the counts in a single
+  vector. Spike times are passed as a cell-array of vectors,
+  containing the spike times in seconds observed in a number of
+  trials, to the function.
+\end{exercise}
+
+\begin{exercise}{plotcounthist.m}{}
+  Write a function that takes a cell-array with spike times as input
+  and plots a normalized histogram of the spike counts counted in
+  windows of a given duration.
+\end{exercise}
+
+\begin{exercise}{plotfanofactor.m}{}
+  Write a function that takes a cell-array with spike times as input
+  and plots in one plot count variances a function of the
+  corresponding mean counts and in a second plot the Fano factor as a
+  function of the duration of the count window in logarithmic
+  scale. Two arguments of the function take the minimum and maximum
+  duration of the count window.
 \end{exercise}