From 20332e2013bdcb5781df804c73dc9c2bfac49068 Mon Sep 17 00:00:00 2001
From: Jan Benda <jan.benda@uni-tuebingen.de>
Date: Mon, 18 Jan 2021 13:13:35 +0100
Subject: [PATCH] [pointprocesses] improved chapter

---
 pointprocesses/code/counthist.m           |  25 +-
 pointprocesses/code/isihist.m             |   2 +-
 pointprocesses/code/isis.m                |   2 -
 pointprocesses/code/isiserialcorr.m       |  16 +-
 pointprocesses/code/plotisiserialcorr.m   |  16 ++
 pointprocesses/code/rasterplot.m          |  30 +-
 pointprocesses/lecture/fanoexamples.py    |  40 ++-
 pointprocesses/lecture/pointprocesses.tex | 330 +++++++++++++---------
 pointprocesses/lecture/rasterexamples.py  |   2 +-
 scientificcomputing-script.tex            |   3 +-
 10 files changed, 272 insertions(+), 194 deletions(-)
 create mode 100644 pointprocesses/code/plotisiserialcorr.m

diff --git a/pointprocesses/code/counthist.m b/pointprocesses/code/counthist.m
index 67b4cae..2f5b746 100644
--- a/pointprocesses/code/counthist.m
+++ b/pointprocesses/code/counthist.m
@@ -1,16 +1,11 @@
-function [counts, bins] = counthist(spikes, w)
-% Compute and plot histogram of spike counts.
+function counthist(spikes, w)
+% Plot histogram of spike counts.
 %
-% [counts, bins] = counthist(spikes, w)
+% counthist(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: the histogram of counts normalized to probabilities
-%   bins: the bin centers for the histogram
-
+%   w: duration of window in seconds for computing the counts
     % collect spike counts:
     tmax = spikes{1}(end);
     n = [];
@@ -21,12 +16,12 @@ function [counts, bins] = counthist(spikes, 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];
+%             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];
+        n = [n, nn];
     end
 
     % histogram of spike counts:
@@ -36,9 +31,7 @@ function [counts, bins] = counthist(spikes, w)
     counts = counts / sum(counts);
 
     % plot:
-    if nargout == 0
-        bar(bins, counts);
-        xlabel('counts k');
-        ylabel('P(k)');
-    end
+    bar(bins, counts);
+    xlabel('counts k');
+    ylabel('P(k)');
 end
diff --git a/pointprocesses/code/isihist.m b/pointprocesses/code/isihist.m
index 00b02c3..a24989c 100644
--- a/pointprocesses/code/isihist.m
+++ b/pointprocesses/code/isihist.m
@@ -8,7 +8,7 @@ function [pdf, centers] = isihist(isis, binwidth)
 %   binwidth: optional width in seconds to be used for the isi bins
 %
 % Returns:
-%   pdf: vector with probability density of interspike intervals in Hz
+%   pdf: vector with pdf of interspike intervals in Hertz
 %   centers: vector with centers of interspikeintervalls in seconds
 
     if nargin < 2
diff --git a/pointprocesses/code/isis.m b/pointprocesses/code/isis.m
index 65a346b..509dcc5 100644
--- a/pointprocesses/code/isis.m
+++ b/pointprocesses/code/isis.m
@@ -5,11 +5,9 @@ function isivec = isis(spikes)
 %
 % Arguments:
 %   spikes: a cell array of vectors of spike times in seconds
-%   isivec: a column vector with all the interspike intervalls
 %
 % Returns:
 %   isivec: a column vector with all the interspike intervalls
-
     isivec = [];
     for k = 1:length(spikes)
         difftimes = diff(spikes{k});
diff --git a/pointprocesses/code/isiserialcorr.m b/pointprocesses/code/isiserialcorr.m
index 21bb474..1b813fc 100644
--- a/pointprocesses/code/isiserialcorr.m
+++ b/pointprocesses/code/isiserialcorr.m
@@ -1,15 +1,15 @@
-function isicorr = isiserialcorr(isivec, maxlag)
+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 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)
@@ -21,14 +21,4 @@ function isicorr = isiserialcorr(isivec, maxlag)
             isicorr(k) = corr(isivec(1:end-lag), isivec(lag+1:end));
         end
     end
-    
-    if nargout == 0
-        % plot:
-        plot(lags, isicorr, '-b');
-        hold on;
-        scatter(lags, isicorr, 50.0, 'b', 'filled');
-        hold off;
-        xlabel('Lag k')
-        ylabel('\rho_k')
-    end
 end
diff --git a/pointprocesses/code/plotisiserialcorr.m b/pointprocesses/code/plotisiserialcorr.m
new file mode 100644
index 0000000..3f21418
--- /dev/null
+++ b/pointprocesses/code/plotisiserialcorr.m
@@ -0,0 +1,16 @@
+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
+    [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
diff --git a/pointprocesses/code/rasterplot.m b/pointprocesses/code/rasterplot.m
index 7b28820..ccbfe60 100644
--- a/pointprocesses/code/rasterplot.m
+++ b/pointprocesses/code/rasterplot.m
@@ -2,21 +2,35 @@ 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 upto tmax 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 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'); 
+        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
-    if tmax < 1.5
+    % 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
diff --git a/pointprocesses/lecture/fanoexamples.py b/pointprocesses/lecture/fanoexamples.py
index 228b473..2e23382 100644
--- a/pointprocesses/lecture/fanoexamples.py
+++ b/pointprocesses/lecture/fanoexamples.py
@@ -91,22 +91,38 @@ def plot_count_fano(ax1, ax2, spikes):
     ax2.xaxis.set_minor_locator(mpt.NullLocator())
     ax2.set_yticks(np.arange(0.0, 1.2, 0.5))
 
+
+def plot_fano(ax, spikes):
+    wins = np.logspace(-2, 0.0, 200)
+    mean_counts = np.zeros(len(wins))
+    var_counts = np.zeros(len(wins))
+    for k, win in enumerate(wins):
+        counts = []
+        for times in spikes:
+            c, _ = np.histogram(times, np.arange(0.0, duration, win))
+            counts.extend(c)
+        mean_counts[k] = np.mean(counts)
+        var_counts[k] = np.var(counts)
+    ax.plot(1000.0*wins, var_counts/mean_counts, **lsB)
+    ax.set_xlabel('Window', 'ms')
+    ax.set_ylim(0.0, 1.2)
+    ax.set_xscale('log')
+    ax.set_xticks([10, 100, 1000])
+    ax.set_xticklabels(['10', '100', '1000'])
+    ax.xaxis.set_minor_locator(mpt.NullLocator())
+    ax.set_yticks(np.arange(0.0, 1.2, 0.5))
+
     
 if __name__ == "__main__":
     homspikes = hompoisson(rate, trials, duration)
     inhspikes = oupifspikes(rate, trials, duration, dt, 0.3, drate, tau)
-    fig, axs = plt.subplots(2, 2)
+    fig, (ax1, ax2) = plt.subplots(1, 2)
     fig.subplots_adjust(**adjust_fs(fig, top=0.5, right=2.0))
-    plot_count_fano(axs[0,0], axs[0,1], homspikes)
-    axs[0,0].text(0.1, 0.95, 'Poisson', transform=axs[0,0].transAxes)
-    axs[0,0].set_xlabel('')
-    axs[0,1].set_xlabel('')
-    axs[0,0].xaxis.set_major_formatter(mpt.NullFormatter())
-    axs[0,1].xaxis.set_major_formatter(mpt.NullFormatter())
-    plot_count_fano(axs[1,0], axs[1,1], inhspikes)
-    axs[1,1].axhline(1.0, **lsGrid)
-    axs[1,0].text(0.1, 0.95, 'OU noise', transform=axs[1,0].transAxes)
-    fig.text(0.01, 0.58, 'Count variance', va='center', rotation='vertical')
-    fig.text(0.51, 0.58, 'Fano factor', va='center', rotation='vertical')
+    plot_fano(ax1, homspikes)
+    ax1.set_ylabel('Fano factor')
+    ax1.text(0.1, 0.95, 'Poisson', transform=ax1.transAxes)
+    plot_fano(ax2, inhspikes)
+    ax2.axhline(1.0, **lsGrid)
+    ax2.text(0.1, 0.95, 'OU noise', transform=ax2.transAxes)
     plt.savefig('fanoexamples.pdf')
     plt.close()
diff --git a/pointprocesses/lecture/pointprocesses.tex b/pointprocesses/lecture/pointprocesses.tex
index f1a57f9..8a297ee 100644
--- a/pointprocesses/lecture/pointprocesses.tex
+++ b/pointprocesses/lecture/pointprocesses.tex
@@ -13,9 +13,9 @@ information. Analyzing the statistics of spike times and their
 relation to sensory stimuli or motor actions is central to
 neuroscientific research. With multi-electrode arrays it is nowadays
 possible to record from hundreds or even thousands of neurons
-simultaneously. The open challenge is how to analyze such data sets in
-order to understand how neural systems work. Let's start with the
-basics in this chapter.
+simultaneously. The open challenge is how to analyze such huge data
+sets in smart ways in order to gain insights into the way neural
+systems work. Let's start with the basics in this chapter.
 
 The result of the pre-processing of electrophysiological recordings
 are series of spike times, which are termed \enterm[spike train]{spike
@@ -27,15 +27,16 @@ process]{Punktprozess}{point processes}.
 
 \begin{figure}[bt]
   \includegraphics[width=1\textwidth]{rasterexamples}
-  \titlecaption{\label{rasterexamplesfig}Raster plot.}{Raster plots of
-    ten trials of data illustrating the times of action
-    potentials. Each vertical stroke illustrates the time at which an
-    action potential was observed. Each row displays the events of one
-    trial. Shown is a stationary point process (homogeneous point
-    process with a rate $\lambda=20$\;Hz, left) and an non-stationary
-    point process with a rate that varies in time (noisy perfect
-    integrate-and-fire neuron driven by Ornstein-Uhlenbeck noise with
-    a time-constant $\tau=100$\,ms, right).}
+  \titlecaption{\label{rasterexamplesfig}Raster plots of spike
+    trains.}{Raster plots of ten trials of data illustrating the times
+    of action potentials. Each vertical stroke illustrates the time at
+    which an action potential was observed. Each row displays the
+    events of one trial. Shown is a stationary point process
+    (homogeneous point process with a rate $\lambda=20$\;Hz, left) and
+    an non-stationary point process with a rate that varies in time
+    (noisy perfect integrate-and-fire neuron driven by
+    Ornstein-Uhlenbeck noise with a time-constant $\tau=100$\,ms,
+    right).}
 \end{figure}
 
 
@@ -59,9 +60,10 @@ process]{Punktprozess}{point processes}.
 \begin{figure}[tb]
   \includegraphics{pointprocesssketch}
   \titlecaption{\label{pointprocesssketchfig} Statistics of point
-    processes.}{A point process is a sequence of instances in time
-    $t_i$ that can be also characterized by inter-event intervals
-    $T_i=t_{i+1}-t_i$ and event counts $n_i$.}
+    processes.}{A temporal point process is a sequence of events in
+    time, $t_i$, that can be also characterized by the corresponding
+    inter-event intervals $T_i=t_{i+1}-t_i$ and event counts $n_i$,
+    i.e. the number of events that occurred so far.}
 \end{figure}
 
 \noindent
@@ -77,75 +79,108 @@ plot in which each vertical line indicates the time of an event. The
 event from two different point processes are shown in
 \figref{rasterexamplesfig}.  In addition to the event times, point
 processes can be described using the intervals $T_i=t_{i+1}-t_i$
-between successive events or the number of observed events within a
-certain time window $n_i$ (\figref{pointprocesssketchfig}).
-
-\begin{exercise}{rasterplot.m}{}
-  Implement a function \varcode{rasterplot()} that displays the times of
-  action potentials within the first \varcode{tmax} seconds in a raster
-  plot. The spike times (in seconds) recorded in the individual trials
-  are stored as vectors of times within a cell array.
-\end{exercise}
+between successive events or the number of observed events $n_i$
+within a certain time window (\figref{pointprocesssketchfig}).
+
+In \enterm[point process!stationary]{stationary} point processes the
+statistics does not change over time. In particular, the rate of the
+process, the number of events per time, is constant. The homogeneous
+point process shown in \figref{rasterexamplesfig} on the left is an
+example of a stationary point process. Although locally within each
+trial there are regions with many events per time and others with long
+intervals between events, the average number of events within a small
+time window over trials does not change in time. In the first sections
+of this chapter we introduce various statistics for characterizing
+stationary point processes.
+
+On the other hand, the example shown in \figref{rasterexamplesfig} on
+the right is a non-stationary point process. The rate of the events in
+all trials first decreases and then increases again. This common
+change in rate may have been caused by some sensory input to the
+neuron. How to estimate temporal changes in event rates (firing rates)
+is covered in section~\ref{nonstationarysec}.
 
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 \section{Homogeneous Poisson process}
 
-The Gaussian distribution is, because of the central limit theorem,
-the standard distribution for continuous measures. The equivalent in
-the realm of point processes is the
-\entermde[distribution!Poisson]{Verteilung!Poisson-}{Poisson distribution}.
+Before we are able to start analyzing point processes, we need some
+data. As before, because we are working with computers, we can easily
+simulate them. To get us started we here briefly introduce the
+homogeneous Poisson process. The Poisson process is to point processes
+what the Gaussian distribution is to the statistics of real-valued
+data. It is a standard against which everything is compared.
 
-In a \entermde[Poisson process!homogeneous]{Poissonprozess!homogener}{homogeneous Poisson
-  process} the events occur at a fixed rate $\lambda=\text{const}$ and
-are independent of both the time $t$ and occurrence of previous events
-(\figref{hompoissonfig}). The probability of observing an event within
-a small time window of width $\Delta t$ is given by
-\begin{equation}
-  \label{hompoissonprob}
-  P = \lambda \cdot \Delta t \; .
-\end{equation}
+In a \entermde[Poisson
+process!homogeneous]{Poissonprozess!homogener}{homogeneous Poisson
+  process} events at every given time, that is within every small time
+window of width $\Delta t$ occur with the same constant
+probability. This probability is independent of absolute time and
+independent of any events occurring before. To observe an event right
+after an event is as likely as to observe an event at some specific
+time later on.
 
-In an \entermde[Poisson process!inhomogeneous]{Poissonprozess!inhomogener}{inhomogeneous Poisson
-  process}, however, the rate $\lambda$ depends on time: $\lambda =
-\lambda(t)$.
-
-\begin{exercise}{poissonspikes.m}{}
-  Implement a function \varcode{poissonspikes()} that uses a homogeneous
-  Poisson process to generate events at a given rate for a certain
-  duration and a number of trials. The rate should be given in Hertz
-  and the duration of the trials is given in seconds. The function
-  should return the event times in a cell-array. Each entry in this
-  array represents the events observed in one trial. Apply
-  \eqnref{hompoissonprob} to generate the event times.
-\end{exercise}
+The simplest way to simulate events of a homogeneous Poisson process
+is based on the exponential distribution \eqref{hompoissonexponential}
+of event intervals. We randomly draw intervals from this distribution
+and then sum them up to convert them to event times. The only
+parameter of a homogeneous Poisson process is its rate. It defines how
+many events per time are expected.
+
+Here is a function that generates several trials of a homogeneous
+Poisson process:
+
+\lstinputlisting[caption={hompoissonspikes.m}]{hompoissonspikes.m}
+
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
+\section{Raster plot}
+
+Let's generate some events with this function and display them in a
+raster plot as in \figref{rasterexamplesfig}. For the raster plot we
+need to draw for each event a line at the corresponding time of the
+event and at the height of the corresponding trial. You can go with a
+one for-loop through the trials and then with another for-loop through
+the event times and every time plot a two-point line. This, however,
+is very slow. The fastest way is to concatenate the coordinates of all
+strokes into large vectors, separate the events by \varcode{nan}
+entries, and pass this to a single call of \varcode{plot()}:
+
+\pageinputlisting[caption={rasterplot.m}]{rasterplot.m}
+
+Adapt this function to your needs and use it where ever possible to
+illustrate your spike train data to your readers. They appreciate
+seeing your raw data and being able to judge the data for themselves
+before you go on analyzing firing rates, interspike-interval
+correlations, etc.
 
-\begin{exercise}{hompoissonspikes.m}{}
-  Implement a function \varcode{hompoissonspikes()} that uses a
-  homogeneous Poisson process to generate spike events at a given rate
-  for a certain duration and a number of trials. The rate should be
-  given in Hertz and the duration of the trials is given in
-  seconds. The function should return the event times in a
-  cell-array. Each entry in this array represents the events observed
-  in one trial. Apply \eqnref{poissonintervals} to generate the event
-  times.
-\end{exercise}
 
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
 \section{Interval statistics}
 
-The intervals $T_i=t_{i+1}-t_i$ between successive events are real
-positive numbers. In the context of action potentials they are
-referred to as \entermde[interspike
+Let's start with the interval statistics of stationary point
+processes.  The intervals $T_i=t_{i+1}-t_i$ between successive events
+are real positive numbers. In the context of action potentials they
+are referred to as \entermde[interspike
 interval]{Interspikeintervall}{interspike intervals}, in short
 \entermde[ISI|see{interspike
-  interval}]{ISI|see{Interspikeintervall}}{ISI}s. The statistics of
-interspike intervals are described using common measures for
-describing the statistics of real-valued variables:
+  interval}]{ISI|see{Interspikeintervall}}{ISI}s.  For analyzing event
+intervals we can use all the usual statistics that we know from
+describing univariate data sets of real numbers:
 
 \begin{figure}[t]
-  \includegraphics[width=0.96\textwidth]{isihexamples}\vspace{-2ex}
+  \includegraphics[width=0.96\textwidth]{isihexamples}
   \titlecaption{\label{isihexamplesfig}Interspike-interval
-    histograms}{of the spike trains shown in \figref{rasterexamplesfig}.}
+    histograms}{of the spike trains shown in
+    \figref{rasterexamplesfig}.  The intervals of the homogeneous
+    Poisson process (left) are exponentially distributed according to
+    \eqnref{hompoissonexponential} (red). Typically for many sensory
+    neurons that fire more regularly is an ISI histogram like the one
+    shown on the right. There is a preferred interval where the
+    distribution peaks. The distribution falls off quickly towards
+    smaller intervals, and very small intervals are absent, probably
+    because of refractoriness of the spike generator. The tail of the
+    distribution again approaches an exponential distribution as for
+    the Poisson process.}
 \end{figure}
 
 \begin{exercise}{isis.m}{}
@@ -168,19 +203,25 @@ describing the statistics of real-valued variables:
     p_{exp}(T) = \lambda e^{-\lambda T}
   \end{equation}
   of a homogeneous Poisson spike train with rate $\lambda$.
-\item Mean interval: $\mu_{ISI} = \langle T \rangle =
-  \frac{1}{n}\sum\limits_{i=1}^n T_i$. The average time it takes from
-  one event to the next. The inverse of the mean interval is identical
-  with the mean rate $\lambda$ (number of events per time, see below)
-  of the process.
-\item Standard deviation of intervals: $\sigma_{ISI} = \sqrt{\langle
-    (T - \langle T \rangle)^2 \rangle}$. Periodically spiking neurons
-  have little variability in their intervals, whereas many cortical
-  neurons cover a wide range with their intervals. The standard
-  deviation of homogeneous Poisson spike trains, $\sigma_{ISI} =
-  \frac{1}{\lambda}$, also equals the inverse rate.  Whether the
-  standard deviation of intervals is low or high, however, is better
-  quantified by the
+\item Mean interval
+  \begin{equation}
+    \label{meanisi}
+    \mu_{ISI} = \langle T \rangle = \frac{1}{n}\sum_{i=1}^n T_i
+  \end{equation}
+  The average time it takes from one event to the next. For stationary
+  point processes the inverse of the mean interval is identical with
+  the mean rate $\lambda$ (number of events per time, see below) of
+  the process.
+\item Standard deviation of intervals
+  \begin{equation}
+    \label{stdisi}
+    \sigma_{ISI} = \sqrt{\langle (T - \langle T \rangle)^2 \rangle}
+  \end{equation}
+  Periodically spiking neurons have little variability in their
+  intervals, whereas many cortical neurons cover a wide range with
+  their intervals. The standard deviation of homogeneous Poisson spike
+  trains also equals the inverse rate.  Whether the standard deviation
+  of intervals is low or high, however, is better quantified by the
 \item \entermde[coefficient of
   variation]{Variationskoeffizient}{Coefficient of variation}, the
   standard deviation of the ISIs relative to their mean:
@@ -188,28 +229,29 @@ describing the statistics of real-valued variables:
     \label{cvisi}
     CV_{ISI} = \frac{\sigma_{ISI}}{\mu_{ISI}}
   \end{equation}
-  Homogeneous Poisson spike trains have an CV of exactly one. The
-  lower the CV the more regularly firing a neuron is firing. CVs
-  larger than one are also possible in spike trains with small
-  intervals separated by really long ones.
+  Homogeneous Poisson spike trains have an $CV$ of exactly one. The
+  lower the $CV$ the more regularly a neuron is firing. $CV$s larger than
+  one are also possible in spike trains with small intervals separated
+  by really long ones.
 %\item \entermde[diffusion coefficient]{Diffusionskoeffizient}{Diffusion coefficient}: $D_{ISI} =
 %  \frac{\sigma_{ISI}^2}{2\mu_{ISI}^3}$.
 \end{itemize}
 
 \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 further returns the
+  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
   probability density as well as the centers of the bins.
 \end{exercise}
 
 \begin{exercise}{plotisihist.m}{}
-  Implement a function that takes the return values of
-  \varcode{isiHist()} as input arguments and then plots the data. The
+  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.
+  standard deviation, and the coefficient of variation of the ISIs
+  (\figref{isihexamplesfig}).
 \end{exercise}
 
 \subsection{Interval correlations}
@@ -227,19 +269,20 @@ correlation coefficients. We form $(x,y)$ data pairs by taking the
 series of intervals $T_i$ as $x$-data values and pairing them with the
 $k$-th next intervals $T_{i+k}$ as $y$-data values. The parameter $k$
 is called \enterm{lag} (\determ{Verz\"ogerung}). For lag one we pair
-each interval with the next one. A \entermde[return map]{return
-  map}{Return map} illustrates dependencies between successive
-intervals by simply plotting the intervals $T_{i+k}$ against the
-intervals $T_i$ in a scatter plot (\figref{returnmapfig}). For Poisson
-spike trains there is no structure beyond the one expected from the
-exponential interspike interval distribution, hinting at neighboring
-interspike intervals being independent of each other. For the spike
-train based on an Ornstein-Uhlenbeck process the return map is more
-clustered along the diagonal, hinting at a positive correlation
-between succeeding intervals. That is, short intervals are more likely
-to be followed by short ones and long intervals more likely by long
-ones. This temporal structure was already clearly visible in the spike
-raster shown in \figref{rasterexamplesfig}.
+each interval with the next one. A \entermde{return map}{return map}
+illustrates dependencies between successive intervals by simply
+plotting the intervals $T_{i+k}$ against the intervals $T_i$ in a
+scatter plot (\figref{returnmapfig}). For Poisson spike trains there
+is no structure beyond the one expected from the exponential
+interspike interval distribution, hinting at neighboring interspike
+intervals being independent of each other. For the spike train based
+on an Ornstein-Uhlenbeck process the return map is more clustered
+along the diagonal, hinting at a positive correlation between
+succeeding intervals. That is, short intervals are more likely to be
+followed by short ones and long intervals more likely by long
+ones. This temporal structure was already clearly visible in each
+trial of the spike raster shown in \figref{rasterexamplesfig} on the
+right.
 
 \begin{figure}[tp]
   \includegraphics[width=1\textwidth]{serialcorrexamples}
@@ -259,17 +302,17 @@ raster shown in \figref{rasterexamplesfig}.
 
 Such dependencies can be further quantified by
 \entermde[correlation!serial]{Korrelation!serielle}{serial
-  correlations}. These are the correlation coefficients between the
-intervals $T_{i+k}$ and $T_i$ in dependence on lag $k$:
+  correlations}. These quantify the correlations between successive
+intervals by Pearson's correlation coefficients between the intervals
+$T_{i+k}$ and $T_i$ in dependence on lag $k$:
 \begin{equation}
   \label{serialcorrelation}
   \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)
 \end{equation}
 The serial correlations $\rho_k$ are usually plotted against the lag
-$k$ for a range small range of lags
-(\figref{returnmapfig}). $\rho_0=1$ is the correlation of each
-interval with itself and always equals one.
+$k$ for a small range of lags (\figref{returnmapfig}). $\rho_0=1$ is
+the correlation of each interval with itself and always equals one.
 
 If the serial correlations all equal zero, $\rho_k =0$ for $k>0$, then
 the length of an interval is independent of all the previous
@@ -282,15 +325,21 @@ thus spike trains may approximate renewal processes.
 
 However, other variables like the intracellular calcium concentration
 or the states of slowly switching ion channels may carry information
-from one interspike interval to the next and thus introducing
-correlations. Such non-renewal dynamics can then be described by the
-non-zero serial correlations (\figref{returnmapfig}).
+from one interspike interval to the next and thus introduce
+correlations between intervals. Such non-renewal dynamics is then
+characterized by the non-zero serial correlations
+(\figref{returnmapfig}).
 
 \begin{exercise}{isiserialcorr.m}{}
-  Implement a function \varcode{isiserialcorr()} that takes a vector of
-  interspike intervals as input argument and calculates the serial
-  correlation. The function should further plot the serial
-  correlation.
+  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.
+\end{exercise}
+
+\begin{exercise}{plotisiserialcorr.m}{}
+  Implement a function \varcode{plotisiserialcorr()} that takes a
+  vector of interspike intervals as input argument and generates a
+  plot of the serial correlations.
 \end{exercise}
 
 
@@ -326,10 +375,11 @@ split into many segments $i$, each of duration $W$, and the number of
 events $n_i$ in each of the segments can be counted. The integer event
 counts can be quantified in the usual ways:
 \begin{itemize}
-\item Histogram of the counts $n_i$. For homogeneous Poisson spike
-  trains with rate $\lambda$ the resulting probability distributions
-  follow a Poisson distribution (\figref{countstatsfig}), where the
-  probability of counting $k$ events within a time window $W$ is given by
+\item Histogram of the counts $n_i$ appropriately normalized to
+  probability distributions. For homogeneous Poisson spike trains with
+  rate $\lambda$ the resulting probability distributions follow a
+  Poisson distribution (\figref{countstatsfig}), where the probability
+  of counting $k$ events within a time window $W$ is given by
   \begin{equation}
     \label{poissondist}
     P(k) = \frac{(\lambda W)^k e^{\lambda W}}{k!}
@@ -338,7 +388,7 @@ counts can be quantified in the usual ways:
 \item Variance of counts:
   $\sigma_n^2 = \langle (n - \langle n \rangle)^2 \rangle$.
 \end{itemize}
-Because spike counts are unitless and positive numbers, the
+Because spike counts are unitless and positive numbers the
 \begin{itemize}
 \item \entermde{Fano Faktor}{Fano factor} (variance of counts divided
   by average count)
@@ -351,7 +401,6 @@ Because spike counts are unitless and positive numbers, the
   homogeneous Poisson processes the Fano factor equals one,
   independently of the time window $W$.
 \end{itemize}
-is an additional measure quantifying event counts.
 
 Note that all of these statistics depend in general on the chosen
 window length $W$. The average spike count, for example, grows
@@ -372,33 +421,33 @@ information encoded in the mean spike count is transmitted.
 
 \begin{figure}[t]
   \includegraphics{fanoexamples}
-  \titlecaption{\label{fanofig}
-    Count variance and Fano factor.}{Variance of event counts as a
-    function of mean counts obtained by varying the duration of the
-    count window (left). Dividing the count variance by the respective
-    mean results in the Fano factor that can be plotted as a function
-    of the count window (right). For Poisson spike trains the variance
-    always equals the mean counts and consequently the Fano factor
-    equals one irrespective of the count window (top). A spike train
-    with positive correlations between interspike intervals (caused by
-    an Ornstein-Uhlenbeck process) has a minimum in the Fano factor,
-    that is an analysis window for which the relative count variance
-    is minimal somewhere close to the correlation time scale of the
-    interspike intervals (bottom).}
+  \titlecaption{\label{fanofig} Fano factor.}{Counting events in time
+    windows of given duration and then dividing the variance of the
+    counts by their mean results in the Fano factor. Here, the Fano
+    factor is plotted as a function of the duration of the window used
+    to count events.  For Poisson spike trains the variance always
+    equals the mean counts and consequently the Fano factor equals one
+    irrespective of the count window (left). A spike train with
+    positive correlations between interspike intervals (caused by an
+    Ornstein-Uhlenbeck process) has a minimum in the Fano factor, that
+    is an analysis window for which the relative count variance is
+    minimal somewhere close to the correlation time scale of the
+    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: (i) a
-  cell-array of vectors containing the spike times in seconds observed
-  in a number of trials, and (ii) the duration of the time window that
-  is used to evaluate the counts.
+  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.
 \end{exercise}
 
 
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 \section{Time-dependent firing rate}
+\label{nonstationarysec}
 
 So far we have discussed stationary spike trains. The statistical properties
 of these did not change within the observation time (stationary point
@@ -566,6 +615,7 @@ relevate time-scale.
 \end{exercise}
 
 \section{Spike-triggered Average}
+\label{stasec}
 
 The graphical representation of the neuronal activity alone is not
 sufficient tot investigate the relation between the neuronal response
diff --git a/pointprocesses/lecture/rasterexamples.py b/pointprocesses/lecture/rasterexamples.py
index 07b8628..7a9f4f0 100644
--- a/pointprocesses/lecture/rasterexamples.py
+++ b/pointprocesses/lecture/rasterexamples.py
@@ -85,7 +85,7 @@ def plot_inhomogeneous_spikes(ax):
 
 
 if __name__ == "__main__":
-    fig, (ax1, ax2) = plt.subplots(1, 2, figsize=cm_size(figure_width, 0.5*figure_width))
+    fig, (ax1, ax2) = plt.subplots(1, 2)
     fig.subplots_adjust(**adjust_fs(fig, left=4.0, right=1.0, top=1.2))
     plot_homogeneous_spikes(ax1)
     plot_inhomogeneous_spikes(ax2)
diff --git a/scientificcomputing-script.tex b/scientificcomputing-script.tex
index 13ca240..016a6a4 100644
--- a/scientificcomputing-script.tex
+++ b/scientificcomputing-script.tex
@@ -116,9 +116,10 @@
 \includechapter{pointprocesses}
 
 % add a chapter on simulating point-processes/spike trains
-% (Poisson spike trains, integrate-and-fire models)
+% (Poisson spike trains, LMP models, integrate-and-fire models, full HH like models)
 
 % add chapter on information theory, mutual information, stimulus reconstruction, coherence
+% move STA here!
 
 % add a chapter on Bayesian inference (the Neuroscience of it and a
 % bit of application for statistical problems).