teaching/scientificComputing
Archived
3
0
Fork 0
Code Issues Pull Requests Releases Wiki Activity

translation of point processes chapter and some tiny fixes elsewhere #1

Merged
jgrewe merged 19 commits from translation into master 2018-11-07 08:38:40 +00:00
Conversation 0 Commits 19 Files Changed 24 +1199 -626
24 changed files with 1199 additions and 626 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

4
bootstrap/lecture/bootstrapsem.py
View File

@ -15,13 +15,13 @@ x = rng.randn(nsamples)
# bootstrap the mean:
mus = []
for i in xrange(nresamples) :
for i in range(nresamples) :
mus.append(np.mean(x[rng.randint(0, nsamples, nsamples)]))
hmus, _ = np.histogram(mus, bins, density=True)
# many SRS:
musrs = []
for i in xrange(nresamples) :
for i in range(nresamples) :
musrs.append(np.mean(rng.randn(nsamples)))
hmusrs, _ = np.histogram(musrs, bins, density=True)

2
bootstrap/lecture/permutecorrelation.py
View File

@ -19,7 +19,7 @@ rd = np.corrcoef(x, y)[0, 1]
# permutation:
nperm = 1000
rs = []
for i in xrange(nperm) :
for i in range(nperm) :
xr=rng.permutation(x)
yr=rng.permutation(y)
rs.append( np.corrcoef(xr, yr)[0, 1] )

3
chapter.mk
View File

@ -8,7 +8,8 @@ PYPDFFILES=$(PYFILES:.py=.pdf)
pythonplots : $(PYPDFFILES)
$(PYPDFFILES) : %.pdf: %.py
python $<
echo $$(which python)
python3 $<
cleanpythonplots :
rm -f $(PYPDFFILES)

2
likelihood/lecture/mlemean.py
View File

@ -41,7 +41,7 @@ for mu in mus :
ax.text(mu-0.1, 0.04, '?', zorder=1, ha='right')
else :
ax.text(mu+0.1, 0.04, '?', zorder=1)
for k in xrange(len(mus)) :
for k in range(len(mus)) :
ax.plot(x, g[:,k], zorder=5)
ax.scatter(xd, 0.05*rng.rand(len(xd))+0.2, s=30, zorder=10)

8
plotting/code/simple_plot.m Normal file
View File

@ -0,0 +1,8 @@
frequency = 5; % frequency of the sine wave in Hz
time = 0.01:0.01:1.0; % the time axis
signal = sin(2 * pi * time * frequency);
plot(time, signal);
xlabel('time [s]');
ylabel('signal');
title('5Hz sine wave')

308
plotting/lecture/plotting.tex
View File

@ -11,123 +11,6 @@ results.
severe.}{\url{www.xkcd.com}}\label{xkcdplotting}
\end{figure}
\section{What makes a good plot?}
Plot should help/enable the interested reader to get a grasp of the
data and to understand the performed analysis and to critically assess
the presented results. The most important rule is the correct and
complete annotation of the plots. This starts with axis labels and
units and and extends to legends. Incomplete annotation can have
terrible consequences (\figref{xkcdplotting}).
The principle of \emph{ink minimization} may be used a a guiding
principle for appealing plots. It requires that the relation of amount
of ink spent on the data and that spent on other parts of the plot
should be strongly in favor of the data. Ornamental of otherwise
unnecessary gimicks should not be used in scientific contexts. An
exception can be made if the particular figure was designed for
didactic purposes and sometimes for presentations.
\begin{important}[Correct labeling of plots]
A data plot must be sufficiently labeled:
\begin{itemize}
\item Every axis must have a label and the correct unit, if it has
one.\\ (e.g. \code[xlabel()]{xlabel('Speed [m/s]'}).
\item When more than one line is plotted, they have to be labeled
using the figure legend, or similar \matlabfun{legend()}.
\item If using subplots that show similar information on the axes,
they should be scaled to show the same ranges to ease comparison
between plots. (e.g. \code[xlim()]{xlim([0 100])}.\\ If one
chooses to ignore this rule one should explicitly state this in
the figure caption and/or the descriptions in the text.
\item Labels must be large enough to be readable. In particular,
when using the figure in a presentation use large enough fonts.
\end{itemize}
\end{important}
\section{Things that should be avoided.}
When plotting scientific data we should take great care to avoid
suggestive or misleading presentations. Unnecessary additions and
fancy graphical effects make a plot frivolous and also violate the
\emph{ink minimization principle}. Illustrations in comic style
(\figref{comicexamplefig}) are not suited for scientific data in most
instances. For presentations or didactic purposes, however, using a
comic style may be helpful to indicate that the figure is a mere
sketch and the exact position of the data points is of no importance.
\begin{figure}[t]
\includegraphics[width=0.7\columnwidth]{outlier}\vspace{-3ex}
\titlecaption{Comic-like illustration.}{Obviously not suited to
present scientific data. In didactic or illustrative contexts they
can be helpful to focus on the important
aspects.}\label{comicexamplefig}
\end{figure}
The following figures show examples of misleading or suggestive
presentations of data. Several of the effects have been exaggerated to
make the point. A little more subtlety these methods are employed to
nudge the viewers experience into the desired direction. You can find
more examples on \url{https://en.wikipedia.org/wiki/Misleading_graph}.
\begin{figure}[p]
\includegraphics[width=0.35\textwidth]{misleading_pie}
\hspace{0.05\textwidth}
\includegraphics[width=0.35\textwidth]{sample_pie}
\titlecaption{Perspective distortion influences the perceived
size.}{By changing the perspective of the 3-D illustration the
highlighted segment \textbf{C} gains more weight than it should
have. In the left graph segments \textbf{A} and \textbf{C} appear
very similar. The 2-D plot on the right-hand side shows that this
is an
illusion. \url{https://en.wikipedia.org/wiki/Misleading_graph}}\label{misleadingpiefig}
\end{figure}
\begin{figure}[p]
\includegraphics[width=0.9\textwidth]{plot_scaling.pdf}
\titlecaption{Choosing the figure format and scaling of the axes
influences the perceived strength of a correlation.}{All subplots
show the same data. By choosing a certain figure size we can
pronounce or reduce the perceived strength of the correlation
in the data. Technically all three plots are correct.
}\label{misleadingscalingfig}
\end{figure}
\begin{figure}[p]
\begin{minipage}[t]{0.3\textwidth}
\includegraphics[width=0.8\textwidth]{improperly_scaled_graph}
\end{minipage}
\begin{minipage}[t]{0.3\textwidth}
\includegraphics[width=0.8\textwidth]{comparison_properly_improperly_graph}
\end{minipage}
\begin{minipage}[t]{0.3\textwidth}
\includegraphics[width=0.7\textwidth]{properly_scaled_graph}
\end{minipage}
\titlecaption{Scaling of markers and symbols.} {In these graphs
symbols have been used to illustrate the measurements made in two
categories. The measured value for category \textbf{B} is actually
three times the measured value for category \textbf{A}. In the
left graph the symbol for category \textbf{B} has been scaled to
triple height while maintaining the proportions. This appears just
fair and correct but leads to the effect that the covered surface
is not increased to the 3-fold but the 9-fold (center plot). The
plot on the right shows how it could have been done correctly.
\url{https://en.wikipedia.org/wiki/Misleading_graph}}\label{misleadingsymbolsfig}
\end{figure}
By using perspective effects in 3-D plot the perceived size can be
distorted into the desired direction. While the plot is correct in a
strict sense it is rather suggestive
(\figref{misleadingpiefig}). Similarly the choice of figure size and
proportions can lead to different interpretations of the
data. Stretching the y-extent of a graph leads to a stronger
impression of the correlation in the data. Compressing this axis will
lead to a much weaker perceived correlation
(\figref{misleadingscalingfig}). When using symbols to illustrate a
quantity we have to take care not to overrate of difference due to
symbol scaling (\figref{misleadingsymbolsfig}).
\section{The \matlab{} plotting system}
Plotting data in \matlab{} is rather straight forward for simple line
@ -191,19 +74,20 @@ the data was changed or the same kind of plot has to be created for a
number of datasets.
\begin{important}[Why manual editing should be avoided.]
On first glance the manual editing of a figure using common tools
like Corel draw, Illustrator, etc.\,appears some much more
convenient and less complex. This, however, is not entirely
true. What if the figure has to be re-drawn or updated? Then the
editing work starts all over again. Rather, there is a great risk
associated with this approach. Axes are shifted, fonts have not been
embedded into the final document, annotations have been copy pasted
between figures and are not valid. All of these mistakes can be
found in publications and then require an erratum, which is not
On first glance the manual editing of a figure using tools such as
Corel draw, Illustrator, etc.\,appears much more convenient and less
complex than coding everything into the analysis scripts. This,
however, is not entirely true. What if the figure has to be re-drawn
or updated? Then the editing work starts all over again. Rather,
there is a great risk associated with the manual editing
approach. Axes may be shifted, fonts have not been embedded into the
final document, annotations have been copy pasted between figures
and are not valid. All of these mistakes can be found in
publications and then require an erratum, which is not
desirable. Even if it appears more cumbersome in the beginning one
should always try to create publication-ready figures directly from
the data analysis tool using scripts or functions to properly
layout the plot.
the data analysis tool using scripts or functions to properly layout
the plot.
\end{important}
\subsection{Simple plotting}
@ -259,6 +143,12 @@ additional options consult the help.
\end{tabular}
\end{table}
The following listing shows a simple line plot with axis labeling and a title
\lstinputlisting[caption={A simple plot showing a sinewave.},
label=simpleplotlisting]{simple_plot.m}
\subsection{Changing properties of a line plot}
The properties of line plots can be changed by passing more arguments
@ -419,14 +309,14 @@ output format (box\,\ref{graphicsformatbox}).
\end{tabular}
\end{minipage}
It is often meaningful to store of data plots generated by \matlab{}
using a vector graphics format. When in doubt they can usually be
It is advisable to store of data plots generated by \matlab{}
using a vector graphics format. In doubt they can usually be
easily converted to a bitmap format. The way from a bitmap to a
vector graphic is not possible without a loss in quality. Storing a
plot that contains a very large set of graphical elements (e.g.\,a
plot that contains very large sets of graphical elements (e.g.\,a
raster-plot showing thousands of action potentials) may, on the
other hand, lead to very large files that can be hard to
handle. Saving such a plot using a bitmap format may be more
handle. Saving such plots using a bitmap format may be more
efficient.
\end{ibox}
@ -584,24 +474,24 @@ its properties. See the \matlab{} help for more information.
\begin{figure}[ht]
\includegraphics[width=0.9\linewidth]{errorbars}
\titlecaption{Adding error bars to a line plot}{\textbf{A}
\titlecaption{Indicating the estimation error in plots.}{\textbf{A}
symmetrical error around the mean (e.g.\ using the standard
deviation). \textbf{B} Errorbars of an asymmetrical distribution
of the data (note: the average value is now the median and the
errors are the lower and upper quartiles). \textbf{C} A shaded
area is used to illustrate the spread of the data. See
listing\,\ref{errorbarlisting}}\label{errorbarplot}
listing\,\ref{errorbarlisting} for A and C and listing\,\ref{errorbarlisting2} }\label{errorbarplot}
\end{figure}
\lstinputlisting[caption={Illustrating estimation errors. Script that
creates \figref{errorbarplot}.},
\lstinputlisting[caption={Illustrating estimation errors using error bars. Script that
creates \figref{errorbarplot}. A, B},
label=errorbarlisting, firstline=13, lastline=29,
basicstyle=\ttfamily\scriptsize]{errorbarplot.m}
\subsubsection{Fill}
For a few years now it has become fancy to illustrate the error not
using errorbars but by drawing a shaded area around the mean. Beside
their fancyness there is also a real argument in favor of using error
the fancyness there is also a real argument in favor of using error
areas instead of errorbars: In case you have a lot of data points with
respective errorbars such that they would merge in the figure it is
cleaner and probably easier to read and handle if one uses an error
@ -613,8 +503,8 @@ with the vertex points of the polygon. For each x-value we now have
two y-values (average minus error and average plus error). Further, we
want the vertices to be connected in a defined order. One can achieve
this by going back and forth on the x-axis; we append a reversed
version of the x-values to the original x-values using the \code{cat}
and inversion is done using the \code{fliplr} command (line 3 in
version of the x-values to the original x-values using \code{cat} and
\code{fliplr} for concatenation and inversion, respectively (line 3 in
listing \ref{errorbarlisting2}; Depending on the layout of your data
you may need concatenate along a different dimension of the data and
use \code{flipud} instead). The y-coordinates of the polygon vertices
@ -625,27 +515,26 @@ property defines the transparency (or rather the opaqueness) of the
area. The provided alpha value is a number between 0 and 1 with zero
leading to invisibility and a value of one to complete
opaqueness. Finally, we use the normal plot command to draw a line
connecting the average values.
connecting the average values (line 12).
\lstinputlisting[caption={Illustrating estimation errors. Script that
creates \figref{errorbarplot}.}, label=errorbarlisting2,
\lstinputlisting[caption={Illustrating estimation errors using a shaded area. Script that
creates \figref{errorbarplot} C.}, label=errorbarlisting2,
firstline=30,
basicstyle=\ttfamily\scriptsize]{errorbarplot.m}
\subsection{Annotations, text}
Sometimes want to highlight certain parts of a plot or simply add an
annotation that does not fit or belong to the legend. In these cases
we can use the \code[text()]{text()} or
\code[annotation()]{annotation()} function to add this information to
the plot. While \varcode{text} simply prints out the given text string
at the defined position (for example line in
The \code[text()]{text()} or \code[annotation()]{annotation()} are
used for highlighting certain parts of a plot or simply adding an
annotation that does not fit or does not belong into the legend.
While \varcode{text} simply prints out the given text string at the
defined position (for example line in
listing\,\ref{regularsubplotlisting}) the \varcode{annotation}
function allows to add some more advanced highlights like arrows,
lines, ellipses, or rectangles. Figure\,\ref{annotationsplot} shows
some examples, the respective code can be found in
listing\,\ref{annotationsplotlisting}. For more options consult the
documentation.
\matlab{} help.
\begin{figure}[ht]
\includegraphics[width=0.5\linewidth]{annotations}
@ -714,20 +603,133 @@ Lissajous figure. The basic steps are:
movie.}, label=animationlisting, firstline=3, lastline=33,
basicstyle=\ttfamily\scriptsize]{movie_example.m}
\section{What makes a good plot?}
Plot should help/enable the interested reader to get a grasp of the
data and to understand the performed analysis and to critically assess
the presented results. The most important rule is the correct and
complete annotation of the plots. This starts with axis labels and
units and and extends to legends. Incomplete annotation can have
terrible consequences (\figref{xkcdplotting}).
The principle of \emph{ink minimization} may be used a a guiding
principle for appealing plots. It requires that the relation of amount
of ink spent on the data and that spent on other parts of the plot
should be strongly in favor of the data. Ornamental of otherwise
unnecessary gimicks should not be used in scientific contexts. An
exception can be made if the particular figure was designed for
didactic purposes and sometimes for presentations.
\begin{important}[Correct labeling of plots]
A data plot must be sufficiently labeled:
\begin{itemize}
\item Every axis must have a label and the correct unit, if it has
one.\\ (e.g. \code[xlabel()]{xlabel('Speed [m/s]'}).
\item When more than one line is plotted, they have to be labeled
using the figure legend, or similar \matlabfun{legend()}.
\item If using subplots that show similar information on the axes,
they should be scaled to show the same ranges to ease comparison
between plots. (e.g. \code[xlim()]{xlim([0 100])}.\\ If one
chooses to ignore this rule one should explicitly state this in
the figure caption and/or the descriptions in the text.
\item Labels must be large enough to be readable. In particular,
when using the figure in a presentation use large enough fonts.
\end{itemize}
\end{important}
\section{Things that should be avoided.}
When plotting scientific data we should take great care to avoid
suggestive or misleading presentations. Unnecessary additions and
fancy graphical effects make a plot frivolous and also violate the
\emph{ink minimization principle}. Illustrations in comic style
(\figref{comicexamplefig}) are not suited for scientific data in most
instances. For presentations or didactic purposes, however, using a
comic style may be helpful to indicate that the figure is a mere
sketch and the exact position of the data points is of no importance.
\begin{figure}[t]
\includegraphics[width=0.7\columnwidth]{outlier}\vspace{-3ex}
\titlecaption{Comic-like illustration.}{Obviously not suited to
present scientific data. In didactic or illustrative contexts they
can be helpful to focus on the important
aspects.}\label{comicexamplefig}
\end{figure}
The following figures show examples of misleading or suggestive
presentations of data. Several of the effects have been exaggerated to
make the point. A little more subtlety these methods are employed to
nudge the viewers experience into the desired direction. You can find
more examples on \url{https://en.wikipedia.org/wiki/Misleading_graph}.
\begin{figure}[p]
\includegraphics[width=0.35\textwidth]{misleading_pie}
\hspace{0.05\textwidth}
\includegraphics[width=0.35\textwidth]{sample_pie}
\titlecaption{Perspective distortion influences the perceived
size.}{By changing the perspective of the 3-D illustration the
highlighted segment \textbf{C} gains more weight than it should
have. In the left graph segments \textbf{A} and \textbf{C} appear
very similar. The 2-D plot on the right-hand side shows that this
is an
illusion. \url{https://en.wikipedia.org/wiki/Misleading_graph}}\label{misleadingpiefig}
\end{figure}
\begin{figure}[p]
\includegraphics[width=0.9\textwidth]{plot_scaling.pdf}
\titlecaption{Choosing the figure format and scaling of the axes
influences the perceived strength of a correlation.}{All subplots
show the same data. By choosing a certain figure size we can
pronounce or reduce the perceived strength of the correlation
in the data. Technically all three plots are correct.
}\label{misleadingscalingfig}
\end{figure}
\begin{figure}[p]
\begin{minipage}[t]{0.3\textwidth}
\includegraphics[width=0.8\textwidth]{improperly_scaled_graph}
\end{minipage}
\begin{minipage}[t]{0.3\textwidth}
\includegraphics[width=0.8\textwidth]{comparison_properly_improperly_graph}
\end{minipage}
\begin{minipage}[t]{0.3\textwidth}
\includegraphics[width=0.7\textwidth]{properly_scaled_graph}
\end{minipage}
\titlecaption{Scaling of markers and symbols.} {In these graphs
symbols have been used to illustrate the measurements made in two
categories. The measured value for category \textbf{B} is actually
three times the measured value for category \textbf{A}. In the
left graph the symbol for category \textbf{B} has been scaled to
triple height while maintaining the proportions. This appears just
fair and correct but leads to the effect that the covered surface
is not increased to the 3-fold but the 9-fold (center plot). The
plot on the right shows how it could have been done correctly.
\url{https://en.wikipedia.org/wiki/Misleading_graph}}\label{misleadingsymbolsfig}
\end{figure}
By using perspective effects in 3-D plot the perceived size can be
distorted into the desired direction. While the plot is correct in a
strict sense it is rather suggestive
(\figref{misleadingpiefig}). Similarly the choice of figure size and
proportions can lead to different interpretations of the
data. Stretching the y-extent of a graph leads to a stronger
impression of the correlation in the data. Compressing this axis will
lead to a much weaker perceived correlation
(\figref{misleadingscalingfig}). When using symbols to illustrate a
quantity we have to take care not to overrate of difference due to
symbol scaling (\figref{misleadingsymbolsfig}).
\section{Summary}
A good plot of scientific data displays the data completely and
seriously without too many distractions. Misleading or suggestive
plots as may result from perspective presentations, inappropriate
scaling of axes of symbols should be avoided.
scaling of axes and symbols should be avoided.
\noindent When combining several line plots within the same figure one should
consider adapting color \textbf{and} line style (solid, dashed,
dotted. etc.) to make the distinguishable even in black-and-white
prints. Combinations of red and green are no good choice since they
prints. Combinations of red and green are not a good choice since they
cannot be distinguished by people with red-green blindness.
\vspace{2ex}
@ -737,6 +739,8 @@ Key ingredients for a good data plot:
\item Complete labeling.
\item Plotted lines and curves must be distinguishable.
\item No suggestive or misleading presentation.
\item The right balance of line width, font size and size of the figure.
\item The right balance of line width, font size and size of the
figure, this may depend on the purpose, for presentations slightly
thicker lines help.
\item Error bars wherever they are appropriate.
\end{itemize}

2
pointprocesses/lecture/firingrates.py
View File

@ -86,7 +86,7 @@ def plot_isi_rate(spike_times, max_t=30, dt=1e-4):
def get_binned_rate(times, bin_width=0.05, max_t=30., dt=1e-4):
time = np.arange(0., max_t, dt)
bins = np.arange(0., max_t, bin_width)
bin_indices = bins / dt
bin_indices = np.asarray(bins / dt, np.int)
hist, _ = sp.histogram(times, bins)
rate = np.zeros(time.shape)

6
pointprocesses/lecture/isihexamples.py
View File

@ -31,7 +31,7 @@ def pifspikes(input, trials, dt, D=0.1) :
times = []
v = vreset
noise = np.sqrt(2.0*D)*np.random.randn(len(input))/np.sqrt(dt)
for k in xrange(len(noise)) :
for k in range(len(noise)) :
v += (input[k]+noise[k])*dt/tau
if v >= vthresh :
v = vreset
@ -41,7 +41,7 @@ def pifspikes(input, trials, dt, D=0.1) :
def isis( spikes ) :
isi = []
for k in xrange(len(spikes)) :
for k in range(len(spikes)) :
isi.extend(np.diff(spikes[k]))
return isi
@ -76,7 +76,7 @@ rng = np.random.RandomState(54637281)
time = np.arange(0.0, duration, dt)
x = np.zeros(time.shape)+rate
n = rng.randn(len(time))*drate*tau/np.sqrt(dt)+rate
for k in xrange(1,len(x)) :
for k in range(1,len(x)) :
x[k] = x[k-1] + (n[k]-x[k-1])*dt/tau
x[x<0.0] = 0.0

672
pointprocesses/lecture/pointprocesses.tex
View File

@ -1,186 +1,186 @@
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\chapter{Analyse von Spiketrains}
\chapter{Spiketrain analysis}
\selectlanguage{ngerman}
\selectlanguage{english}
\determ[Aktionspotential]{Aktionspotentiale} (\enterm{spikes}) sind die Tr\"ager der
Information in Nervensystemen. Dabei ist in erster Linie nur der
Zeitpunkt des Auftretens eines Aktionspotentials von Bedeutung. Die
genaue Form des Aktionspotentials spielt keine oder nur eine
untergeordnete Rolle.
\enterm[Actionspotentials]{Actionspotentials} (\enterm{spikes}) are
the carriers of information in the nervous system. Thereby it is
mainly the time at which the spikes are generated that is of
importance. The waveform of the action potential is largely
stereotyped and does not carry information.
Nach etwas Vorverarbeitung haben elektrophysiologische Messungen
deshalb Listen von Spikezeitpunkten als Ergebniss --- sogenannte
\enterm{spiketrains}. Diese Messungen k\"onnen wiederholt werden und
es ergeben sich mehrere \enterm{trials} von Spiketrains
(\figref{rasterexamplesfig}).
The result of the processing of electrophysiological recordings are
series of spike times, which are then termed \enterm{spiketrains}. If
measurements are repeated we yield several \enterm{trials} of
spiketrains (\figref{rasterexamplesfig}).
Spiketrains sind Zeitpunkte von Ereignissen --- den Aktionspotentialen
--- und deren Analyse f\"allt daher in das Gebiet der Statistik von
sogenannten \determ[Punktprozess]{Punktprozessen}.
Spiketrains are times of events, the action potentials. The analysis
of these leads into the realm of the so called \enterm[point
process]{point processes}.
\begin{figure}[ht]
\includegraphics[width=1\textwidth]{rasterexamples}
\titlecaption{\label{rasterexamplesfig}Raster-Plot.}{Raster-Plot von
jeweils 10 Realisierungen eines station\"arenen Punktprozesses
(homogener Poisson Prozess mit Rate $\lambda=20$\;Hz, links) und
eines nicht-station\"aren Punktprozesses (perfect
integrate-and-fire Neuron getrieben mit Ohrnstein-Uhlenbeck
Rauschen mit Zeitkonstante $\tau=100$\,ms, rechts). Jeder
vertikale Strich markiert den Zeitpunkt eines Ereignisses.
Jede Zeile zeigt die Ereignisse eines trials.}
\titlecaption{\label{rasterexamplesfig}Raster-plot.}{Raster-plot of
ten realizations of a stationary point process (homogeneous point
process with a rate $\lambda=20$\;Hz, left) and an inhomogeneous
point process (perfect integrate-and-fire neuron dirven by
Ohrnstein-Uhlenbeck noise with a time-constant $\tau=100$\,ms,
right). Each vertical dash illustrates the time at which the
action potential was observed. Each line represents the event of
each trial.}
\end{figure}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Punktprozesse}
\section{Point processes}
Ein zeitlicher Punktprozess (\enterm{point process}) ist ein
stochastischer Prozess, der eine Abfolge von Ereignissen zu den Zeiten
$\{t_i\}$, $t_i \in \reZ$, generiert.
A temporal \enterm{point process} is a stochastic process that
generates a sequence of events at times $\{t_i\}$, $t_i \in
\reZ$.
\begin{ibox}{Beispiele von Punktprozessen}
Jeder Punktprozess wird durch einen sich in der Zeit kontinuierlich
entwickelnden Prozess generiert. Wann immer dieser Prozess eine
Schwelle \"uberschreitet wird ein Ereigniss des Punktprozesses
erzeugt. Zum Beispiel:
\begin{itemize}
\item Aktionspotentiale/Herzschlag: wird durch die Dynamik des
Membranpotentials eines Neurons/Herzzelle erzeugt.
\item Erdbeben: wird durch die Dynamik des Druckes zwischen
tektonischen Platten auf beiden Seiten einer geologischen Verwerfung
erzeugt.
\item Zeitpunkt eines Grillen/Frosch/Vogelgesangs: wird durch die
Dynamik des Nervensystems und des Muskelapparates erzeugt.
\end{itemize}
\begin{ibox}{Examples of point processes}
Every point process is generated by a temporally continuously
developing process. An event is generated whenever this process
reaches a certain threshold. For example:
\begin{itemize}
\item Action potentials/heart beat: created by the dynamics of the
neuron/sinoatrial node
\item Earthquake: defined by the dynamics of the pressure between
tectonical plates.
\item Evoked communication calls in crickets/frogs/birds: shaped by
the dynamics of nervous system and the muscle appartus.
\end{itemize}
\end{ibox}
\begin{figure}[t]
\texpicture{pointprocessscetch}
\titlecaption{\label{pointprocessscetchfig} Statistik von
Punktprozessen.}{Ein Punktprozess ist eine Abfolge von
Zeitpunkten $t_i$ die auch durch die Intervalle $T_i=t_{i+1}-t_i$
oder die Anzahl der Ereignisse $n_i$ beschrieben werden kann. }
\titlecaption{\label{pointprocessscetchfig} Statistics of point
processes.}{A point process is a sequence of instances in time
$t_i$ that can be characterized through the inter-event-intervals
$T_i=t_{i+1}-t_i$ and the number of events $n_i$. }
\end{figure}
F\"ur die Neurowissenschaften ist die Statistik der Punktprozesse
besonders wichtig, da die Zeitpunkte der Aktionspotentiale als
zeitlicher Punktprozess betrachtet werden k\"onnen und entscheidend
f\"ur die Informations\"ubertragung sind.
Bei Punktprozessen k\"onnen wir die Zeitpunkte $t_i$ ihres Auftretens,
die Intervalle zwischen diesen Zeitpunkten $T_i=t_{i+1}-t_i$, sowie
die Anzahl der Ereignisse $n_i$ bis zu einer bestimmten Zeit betrachten
(\figref{pointprocessscetchfig}).
In the neurosciences, the statistics of point processes is of
importance since the timing of the neuronal events (the action
potentials) is crucial for information transmission and can be treated
as such a process.
Zwei Punktprozesse mit verschiedenen Eigenschaften sind in
\figref{rasterexamplesfig} als Rasterplot dargestellt, bei dem die
Zeitpunkte der Ereignisse durch senkrechte Striche markiert werden.
Point processes can be described using the intervals between
successive events $T_i=t_{i+1}-t_i$ and the number of observed events
within a certain time window $n_i$ (\figref{pointprocessscetchfig}).
The events originating from a point process can be illustrated in form
of a scatter- or raster plot in which each vertical line indicates the
time of an event. The event from two different point processes are
shown in \figref{rasterexamplesfig}.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Intervallstatistik}
\section{Intervalstatistics}
Die Intervalle $T_i=t_{i+1}-t_i$ zwischen aufeinanderfolgenden
Ereignissen sind reelle, positive Zahlen. Bei Aktionspotentialen
heisen die Intervalle auch \determ{Interspikeintervalle}
(\enterm{interspike intervals}). Deren Statistik kann mit den
\"ublichen Gr\"o{\ss}en beschrieben werden.
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 \enterm{interspike intervals}. The statistics of these
are described using the common measures.
\begin{figure}[t]
\includegraphics[width=0.96\textwidth]{isihexamples}\vspace{-2ex}
\titlecaption{\label{isihexamplesfig}Interspikeintervall Histogramme}{der in
\figref{rasterexamplesfig} gezeigten Spikes.}
\titlecaption{\label{isihexamplesfig}Interspike interval
histogram}{of the spikes depicted in \figref{rasterexamplesfig}.}
\end{figure}
\begin{exercise}{isis.m}{}
Schreibe eine Funktion \code{isis()}, die aus mehreren trials von Spiketrains die
Interspikeintervalle bestimmt und diese in einem Vektor
zur\"uckgibt. Jeder trial der Spiketrains ist ein Vektor mit den
Spikezeiten gegeben in Sekunden als Element in einem \codeterm{cell-array}.
Implement a function \code{isis()} that calculates the interspike
intervals from several spike trains. The function should return a
single vector of intervals. The action potentials recorded in the
individual trials are stored as vectors of spike times within a
\codeterm{cell-array}. Spike times are given in seconds.
\end{exercise}
\subsection{Intervallstatistik erster Ordnung}
\subsection{First order interval statistics}
\begin{itemize}
\item Wahrscheinlichkeitsdichte $p(T)$ der Intervalle $T$
(\figref{isihexamplesfig}). Normiert auf $\int_0^{\infty} p(T) \; dT
\item Probability density $p(T)$ of the intervals $T$
(\figref{isihexamplesfig}). Normalized to $\int_0^{\infty} p(T) \; dT
= 1$.
\item Mittleres Intervall: $\mu_{ISI} = \langle T \rangle =
\item Average interval: $\mu_{ISI} = \langle T \rangle =
\frac{1}{n}\sum\limits_{i=1}^n T_i$.
\item Standardabweichung der Intervalle: $\sigma_{ISI} = \sqrt{\langle (T - \langle T
\rangle)^2 \rangle}$\vspace{1ex}
\item \determ{Variationskoeffizient} (\enterm{coefficient of variation}): $CV_{ISI} =
\item Standard deviation of the interspike intervals: $\sigma_{ISI} = \sqrt{\langle (T - \langle T
\rangle)^2 \rangle}$\vspace{1ex}
\item \enterm{Coefficient of variation}: $CV_{ISI} =
\frac{\sigma_{ISI}}{\mu_{ISI}}$.
\item \determ{Diffusionskoeffizient} (\enterm{diffusion coefficient}): $D_{ISI} =
\item \enterm{Diffusion coefficient}: $D_{ISI} =
\frac{\sigma_{ISI}^2}{2\mu_{ISI}^3}$.
\end{itemize}
\begin{exercise}{isihist.m}{}
Schreibe eine Funktion \code{isiHist()}, die einen Vektor mit Interspikeintervallen
entgegennimmt und daraus ein normiertes Histogramm der Interspikeintervalle
berechnet.
Implement a function \code{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
probability density as well as the centers of the bins.
\end{exercise}
\begin{exercise}{plotisihist.m}{}
Schreibe eine Funktion, die die Histogrammdaten der Funktion
\code{isiHist()} entgegennimmt, um das Histogramm zu plotten. Im
Plot sollen die Interspikeintervalle in Millisekunden aufgetragen
werden. Das Histogramm soll zus\"atzlich mit Mittelwert,
Standardabweichung und Variationskoeffizient der
Interspikeintervalle annotiert werden.
Implement a function that takes the return values of
\code{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.
\end{exercise}
\subsection{Korrelationen der Intervalle}
In \enterm{return maps} werden die um das \enterm{lag} $k$ verz\"ogerten
Intervalle $T_{i+k}$ gegen die Intervalle $T_i$ geplottet. Dies macht
m\"ogliche Abh\"angigkeiten von aufeinanderfolgenden Intervallen
sichtbar.
\subsection{Interval correlations}
So called \enterm{return maps} are used to illustrate
interdependencies between successive interspike intervals. The return
map plots the delayed interval $T_{i+k}$ against the interval
$T_i$. The parameter $k$ is called the \enterm{lag} $k$. Stationary
and non-stationary return maps are distinctly different
\figref{returnmapfig}.
\begin{figure}[t]
\includegraphics[width=1\textwidth]{returnmapexamples}
\includegraphics[width=1\textwidth]{serialcorrexamples}
\titlecaption{\label{returnmapfig}Interspikeintervall return maps und
serielle Korrelationen}{zwischen aufeinander folgenden Intervallen
im Abstand des Lags $k$.}
\titlecaption{\label{returnmapfig}Interspike interval analyses of a
stationary and a non-stationary pointprocess.}{Upper plots show the
return maps and the lower panels depict the serial correlation of
successive intervals separated by the lag $k$.}
\end{figure}
Solche Ab\"angigkeiten werden durch die \determ{serielle
Korrelationen} (\enterm{serial correlations}) der Intervalle
quantifiziert. Das ist der \determ{Korrelationskoeffizient} zwischen
aufeinander folgenden Intervallen getrennt durch lag $k$:
Such dependencies can be further quantified using the \enterm{serial
correlations} \figref{returnmapfig}. The serial correlation is the
correlation coefficient of the intervals $T_i$ and the intervals
delayed by the lag $T_{i+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) \]
\"Ublicherweise wird die Korrelation $\rho_k$ gegen den Lag $k$
aufgetragen (\figref{returnmapfig}). $\rho_0=1$ (Korrelation jedes
Intervalls mit sich selber).
= {\rm corr}(T_{i+k}, T_i) \] The resulting correlation coefficient
$\rho_k$ is usually plotted against the lag $k$
\figref{returnmapfig}. $\rho_0=1$ is the correlation of each interval
with itself and is always 1.
\begin{exercise}{isiserialcorr.m}{}
Schreibe eine Funktion \code{isiserialcorr()}, die einen Vektor mit Interspikeintervallen
entgegennimmt und daraus die seriellen Korrelationen berechnet und plottet.
\pagebreak[4]
Implement a function \code{isiserialcorr()} that takes a vector of
interspike intervals as input argument and calculates the serial
correlation. The function should further plot the serial
correlation. \pagebreak[4]
\end{exercise}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Z\"ahlstatistik}
\section{Count statistics}
% \begin{figure}[t]
% \includegraphics[width=0.48\textwidth]{poissoncounthist100hz10ms}\hfill
% \includegraphics[width=0.48\textwidth]{poissoncounthist100hz100ms}
% \titlecaption{\label{countstatsfig}Count Statistik.}{}
% \end{figure}
Die Anzahl der Ereignisse $n_i$ in Zeifenstern $i$ der
L\"ange $W$ ergeben ganzzahlige, positive Zufallsvariablen, die meist
durch folgende Sch\"atzer charakterisiert werden:
The number of events $n_i$ (counts) in a time window $i$ of the duration $W$
yields positive integer random numbers that are commonly quantified
using the following measures:
\begin{itemize}
\item Histogramm der counts $n_i$.
\item Mittlere Anzahl von Ereignissen: $\mu_N = \langle n \rangle$.
\item Varianz der Anzahl: $\sigma_n^2 = \langle (n - \langle n \rangle)^2 \rangle$.
\item \determ{Fano Faktor} (Varianz geteilt durch Mittelwert): $F = \frac{\sigma_n^2}{\mu_n}$.
\item Histogram of the counts $n_i$.
\item Average number of events: $\mu_N = \langle n \rangle$.
\item Variance of the counts: $\sigma_n^2 = \langle (n - \langle n \rangle)^2 \rangle$.
\item \determ{Fano Faktor} (The variance divided by the average): $F = \frac{\sigma_n^2}{\mu_n}$.
\end{itemize}
Insbesondere ist die mittlere Rate der Ereignisse $r$ (Spikes pro
Zeit, \determ{Feuerrate}) gemessen in Hertz \sindex[term]{Feuerrate!mittlere Rate}
And in particular the average firing rate $r$ (spike count per time interval
, \determ{Feuerrate}) that is given in Hertz \sindex[term]{Feuerrate!mittlere Rate}
\begin{equation}
\label{firingrate}
r = \frac{\langle n \rangle}{W} \; .
@ -200,315 +200,333 @@ Zeit, \determ{Feuerrate}) gemessen in Hertz \sindex[term]{Feuerrate!mittlere Rat
% \end{figure}
\begin{exercise}{counthist.m}{}
Schreibe eine Funktion \code{counthist()}, die aus mehreren trials
von Spiketrains die Verteilung der Anzahl der Spikes in Fenstern
einer der Funktion \"ubergegebenen Breite bestimmt, das Histogramm
plottet und zur\"uckgibt. Jeder trial der Spiketrains ist ein Vektor
mit den Spikezeiten gegeben in Sekunden als Element in einem
\codeterm{cell-array}.
Implement a function \code{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
\codeterm{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.
\end{exercise}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Homogener Poisson Prozess}
F\"ur kontinuierliche Me{\ss}gr\"o{\ss}en ist die Normalverteilung
u.a. wegen dem Zentralen Grenzwertsatz die Standardverteilung. Eine
\"ahnliche Rolle spielt bei Punktprozessen der \determ{Poisson
Prozess}.
Beim \determ[Poisson Prozess!homogener]{homogenen Poisson Prozess}
treten Ereignisse mit einer festen Rate $\lambda=\text{const.}$ auf
und sind unabh\"angig von der Zeit $t$ und unabh\"angig von den
Zeitpunkten fr\"uherer Ereignisse (\figref{hompoissonfig}). Die
Wahrscheinlichkeit zu irgendeiner Zeit ein Ereigniss in einem kleinen
Zeitfenster der Breite $\Delta t$ zu bekommen ist
\section{Homogeneous Poisson process}
The Gaussian distribution is, due to the central limit theorem, the
standard for continuous measures. The equivalent in the realm of point
processes is the \enterm{Poisson distribution}.
In a \enterm[Poisson process!homogeneous]{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 even within a
small time window of width $\Delta t$ is given by
\begin{equation}
\label{hompoissonprob}
P = \lambda \cdot \Delta t \; .
P = \lambda \cdot \Delta t \; .
\end{equation}
Beim \determ[Poisson Prozess!inhomogener]{inhomogenen Poisson Prozess}
h\"angt die Rate $\lambda$ von der Zeit ab: $\lambda = \lambda(t)$.
In an \enterm[Poisson process!inhomogeneous]{inhomogeneous Poisson
process}, however, the rate $\lambda$ depends on the time: $\lambda =
\lambda(t)$.
\begin{exercise}{poissonspikes.m}{}
Schreibe eine Funktion \code{poissonspikes()}, die die Spikezeiten
eines homogenen Poisson-Prozesses mit gegebener Rate in Hertz f\"ur
eine Anzahl von trials gegebener maximaler L\"ange in Sekunden in
einem \codeterm{cell-array} zur\"uckgibt. Benutze \eqnref{hompoissonprob}
um die Spikezeiten zu bestimmen.
Implement a function \code{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}
\begin{figure}[t]
\includegraphics[width=1\textwidth]{poissonraster100hz}
\titlecaption{\label{hompoissonfig}Rasterplot von Spikes eines homogenen
Poisson Prozesses mit $\lambda=100$\,Hz.}{}
\titlecaption{\label{hompoissonfig}Rasterplot of spikes of a
homogeneous Poisson process with a rate $\lambda=100$\,Hz.}{}
\end{figure}
\begin{figure}[t]
\includegraphics[width=0.45\textwidth]{poissonisihexp20hz}\hfill
\includegraphics[width=0.45\textwidth]{poissonisihexp100hz}
\titlecaption{\label{hompoissonisihfig}Interspikeintervallverteilungen
zweier Poissonprozesse.}{}
\titlecaption{\label{hompoissonisihfig}Distribution of interspike intervals of two Poisson processes.}{}
\end{figure}
Der homogene Poissonprozess hat folgende Eigenschaften:
The homogeneous Poisson process has the following properties:
\begin{itemize}
\item Die Intervalle $T$ sind exponentiell verteilt (\figref{hompoissonisihfig}):
\item Intervals $T$ are exponentially distributed (\figref{hompoissonisihfig}):
\begin{equation}
\label{poissonintervals}
p(T) = \lambda e^{-\lambda T} \; .
\end{equation}
\item Das mittlere Intervall ist $\mu_{ISI} = \frac{1}{\lambda}$ .
\item Die Varianz der Intervalle ist $\sigma_{ISI}^2 = \frac{1}{\lambda^2}$ .
\item Der Variationskoeffizient ist also immer $CV_{ISI} = 1$ .
\item Die \determ[serielle Korrelationen]{seriellen Korrelationen}
$\rho_k =0$ f\"ur $k>0$, da das Auftreten der Ereignisse
unabh\"angig von der Vorgeschichte ist. Ein solcher Prozess wird
auch \determ{Erneuerungsprozess} genannt (\enterm{renewal process}).
\item Die Anzahl der Ereignisse $k$ innerhalb eines Fensters der
L\"ange W ist \determ[Poisson-Verteilung]{Poissonverteilt}:
\item The average interval is $\mu_{ISI} = \frac{1}{\lambda}$ .
\item The variance of the intervals is $\sigma_{ISI}^2 = \frac{1}{\lambda^2}$ .
\item Thus, the coefficient of variation is always $CV_{ISI} = 1$ .
\item The serial correlation is $\rho_k =0$ for $k>0$, since the
occurrence of an event is independent of all previous events. Such a
process is also called a \enterm{renewal process}.
\item The number of events $k$ within a temporal window of duration
$W$ is Poisson distributed:
\[ P(k) = \frac{(\lambda W)^ke^{\lambda W}}{k!} \]
(\figref{hompoissoncountfig})
\item Der \determ{Fano Faktor} ist immer $F=1$ .
(\figref{hompoissoncountfig})
\item The Fano Faktor is always $F=1$ .
\end{itemize}
\begin{exercise}{hompoissonspikes.m}{}
Schreibe eine Funktion \code{hompoissonspikes()}, die die Spikezeiten
eines homogenen Poisson-Prozesses mit gegebener Rate in Hertz f\"ur
eine Anzahl von trials gegebener maximaler L\"ange in Sekunden in
einem \codeterm{cell-array} zur\"uckgibt. Benutze die exponentiell-verteilten
Interspikeintervalle \eqnref{poissonintervals}, um die Spikezeiten zu erzeugen.
Implement a function \code{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}
\begin{figure}[t]
\includegraphics[width=0.48\textwidth]{poissoncounthistdist100hz10ms}\hfill
\includegraphics[width=0.48\textwidth]{poissoncounthistdist100hz100ms}
\titlecaption{\label{hompoissoncountfig}Z\"ahlstatistik von Poisson Spikes.}{}
\titlecaption{\label{hompoissoncountfig}Count statistics of Poisson
spiketrains.}{}
\end{figure}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Zeitabh\"angige Feuerraten}
Bisher haben wir station\"are Spiketrains betrachtet, deren Statistik
sich innerhalb der Analysezeit nicht ver\"andert (station\"are
Punktprozesse). Meistens jedoch \"andert sich die Statistik der
Spiketrains eines Neurons mit der Zeit. Z.B. kann ein sensorisches
Neuron auf einen Reiz hin mit einer erh\"ohten Feuerrate antworten
(nichtstation\"arer Punktprozess).
Wie die mittlere Anzahl der Spikes sich mit der Zeit ver\"andert, die
\determ{Feuerrate} $r(t)$, ist die wichtigste Gr\"o{\ss}e bei
nicht-station\"aren Spiketrains. Die Einheit der Feuerrate ist Hertz,
also Anzahl Aktionspotentiale pro Sekunde. Es gibt verschiedene
Methoden diese zu bestimmen. Drei solcher Methoden sind in Abbildung
\ref{psthfig} dargestellt. Alle Methoden haben ihre Berechtigung und
ihre Vor- und Nachteile. Im folgenden werden die drei Methoden aus
Abbildung \ref{psthfig} n\"aher erl\"autert.
\section{Time-dependent firing rate}
So far we have discussed stationary spiketrains. The statistical properties
of these did not change within the observation time (stationary point
processes). Most commonly, however, this is not the case. A sensory
neuron, for example, might respond to a stimulus by modulating its
firing rate (non-stationary point process).
How the firing rate $r(t)$ changes over time is the most important
measure, when analyzing non-stationary spike trains. The unit of the
firing rate is Hertz, i.e. the number of action potentials per
second. There are different ways to estimate the firing rate and three
of these methods are illustrated in \figref{psthfig}. All of
these have their own justifications and pros- and cons. In the
following we will discuss these methods more
closely.
\begin{figure}[tp]
\includegraphics[width=\columnwidth]{firingrates}
\titlecaption{Bestimmung der zeitabh\"angigen
Feuerrate.}{\textbf{A)} Rasterplot eines Spiketrains. \textbf{B)}
Feurerrate aus der instantanen Feuerrate bestimmt. \textbf{C)}
klassisches PSTH mit der Binning Methode. \textbf{D)} Feuerrate
durch Faltung mit einem Gauss Kern bestimmt.}\label{psthfig}
\titlecaption{Estimating the time-dependent firing
rate.}{\textbf{A)} Rasterplot depicting the spiking activity of a
neuron. \textbf{B)} Firing rate calculated from the
\emph{instantaneous rate}. \textbf{C)} classical PSTH with the
\emph{binning} method. \textbf{D)} Firing rate estimated by
\emph{convolution} of the activity with a Gaussian
kernel.}\label{psthfig}
\end{figure}
\subsection{Instantane Feuerrate}
\subsection{Instantaneous firing rate}
\begin{figure}[tp]
\includegraphics[width=\columnwidth]{isimethod}
\titlecaption{Instantane Feuerrate.}{Skizze eines Spiketrains
(oben). Die Pfeile zwischen aufeinanderfolgenden
Aktionspotentialen mit den Zahlen in Millisekunden illustrieren
die Interspikeintervalle. Der Kehrwert des Interspikeintervalle
ergibt die instantane Feuerrate.}\label{instrate}
\titlecaption{Instantaneous firing rate.}{Sketch of the recorded
spiketrain (top). Arrows illustrate the interspike intervals and
number give the intervals in milliseconds. The inverse of the
interspike interval is the \emph{instantaneous firing
rate}.}\label{instratefig}
\end{figure}
Ein sehr einfacher Weg, die zeitabh\"angige Feuerrate zu bestimmen ist
die sogenannte \determ[Feuerrate!instantane]{instantane Feuerrate}
(\enterm[firing rate!instantaneous]{instantaneous firing rate}). Dabei
wird die Feuerrate aus dem Kehrwert der Interspikeintervalle, der Zeit
zwischen zwei aufeinander folgenden Aktionspotentialen
(\figref{instrate} A), bestimmt. Die abgesch\"atzte Feuerrate
(\figref{instrate} B) ist g\"ultig f\"ur das gesammte
Interspikeintervall. Diese Methode hat den Vorteil, dass sie sehr
einfach zu berechnen ist und keine Annahme \"uber eine relevante
Zeitskala (der Kodierung oder des Auslesemechanismus der
postsynaptischen Zelle) macht. $r(t)$ ist allerdings keine
kontinuierliche Funktion, die Spr\"unge in der Feuerrate k\"onnen
f\"ur manche Analysen nachteilig sein. Au{\ss}erdem wird die Feuerrate
nie gleich Null, auch wenn lange keine Aktionspotentiale generiert
wurden.
A very simple method for estimating the time-dependent firing rate is
the \enterm[firing rate!instantaneous]{instantaneous firing rate}. The
firing rate can be directly estimated as the inverse of the time
between successive spikes, the interspike-interval
(\figref{instratefig}).
\begin{equation}
\label{instantaneousrateeqn}
r_i = \frac{1}{T_i} .
\end{equation}
The instantaneous rate $r_i$ is valid for the whole interspike
interval. The method has the advantage of being extremely easy to
compute and that it does not make any assumptions about the relevant
timescale (of the encoding in the neuron or the decoding of a
postsynaptic neuron). The resulting $r(t)$, however, is no continuous
function, the firing rate jumps from one level to the next. Since the
interspike interval between successive spikes is never infinitely
long, the firing rate never reaches zero despite that the neuron may
not fire an action potential for a long time.
\begin{exercise}{instantaneousRate.m}{}
Implementiere die Absch\"atzung der Feuerrate auf Basis der
instantanen Feuerrate. Plotte die Feuerrate als Funktion der Zeit.
Implement a function that computes the instantaneous firing
rate. Plot the firing rate as a function of time.
\note{TODO: example data!!!}
\end{exercise}
\subsection{Peri-Stimulus-Zeit-Histogramm}
W\"ahrend die Instantane Rate den Kehrwert der Zeit von einem bis zum
n\"achsten Aktionspotential misst, sch\"atzt das sogenannte
\determ{Peri-Stimulus-Zeit-Histogramm} (\enterm{peri stimulus time
histogram}, \determ[PSTH|see{Peri-Stimulus-Zeit-Histogramm}]{PSTH})
die Wahrscheinlichkeit ab, zu einem Zeitpunkt Aktionspotentiale
anzutreffen. Es wird versucht die mittlere Rate \eqnref{firingrate} im
Grenzwert kleiner Beobachtungszeiten abzusch\"atzen:
\subsection{Peri-stimulus-time-histogram}
While the \emph{instantaneous firing rate} uses the interspike
interval, the \enterm{peri stimulus time histogram} (PSTH) uses the
spike count within observation windows of the duration $W$. It
estimates the probability of observing a spike within that observation
time. It tries to estimat the average rate in the limit of small
obersvation times \eqnref{psthrate}:
\begin{equation}
\label{psthrate}
r(t) = \lim_{W \to 0} \frac{\langle n \rangle}{W} \; ,
\end{equation}
wobei die Anzahl $n$ der Aktionspotentiale, die im Zeitintervall
$(t,t+W)$ aufgetreten sind, \"uber trials gemittelt wird. Eine solche
Rate enspricht der zeitabh\"angigen Rate $\lambda(t)$ des inhomogenen
Poisson-Prozesses.
where $\langle n \rangle$ is the across trial average number of action
potentials observed within the interval $(t, t+W)$. Such description
matches the time-dependent firing rate $\lambda(t)$ of an
inhomogeneous Poisson process.
Das PSTH \eqnref{psthrate} kann entweder \"uber die Binning-Methode
oder durch Verfaltung mit einem Kern bestimmt werden. Beiden Methoden
gemeinsam ist die Notwendigkeit der Wahl einer zus\"atzlichen Zeitskala,
die der Beobachtungszeit $W$ in \eqnref{psthrate} entspricht.
The firing probability can be estimated using the \emph{binning
method} or by using \emph{kernel density estimations}. Both methods
make an assumption about the relevant observation time-scale ($W$ in
\eqnref{psthrate}).
\subsubsection{Binning-Methode}
\subsubsection{Binning-method}
\begin{figure}[tp]
\includegraphics[width=\columnwidth]{binmethod}
\titlecaption{Bestimmung des PSTH mit der Binning Methode.}{Der
gleiche Spiketrain wie in \figref{instrate}. Die grauen Linien
markieren die Grenzen der Bins und die Zahlen geben die Anzahl der Spikes
in jedem Bin an (oben). Die Feuerrate ergibt sich aus dem
mit der Binbreite normierten Zeithistogramm (unten).}\label{binpsth}
\titlecaption{Estimating the PSTH using the binning method.}{ The
same spiketrain as shown in \figref{instratefig} is used
here. Vertical gray lines indicate the borders between adjacent
bins in which the number of action potentials is counted (red
numbers). The firing rate is then the histogram normalized to the
binwidth.}\label{binpsthfig}
\end{figure}
Bei der Binning-Methode wird die Zeitachse in gleichm\"aßige
Abschnitte (Bins) eingeteilt und die Anzahl Aktionspotentiale, die in
die jeweiligen Bins fallen, gez\"ahlt (\figref{binpsth} A). Um diese
Z\"ahlungen in die Feuerrate umzurechnen muss noch mit der Binweite
normiert werden. Das ist \"aquivalent zur Absch\"atzung einer
Wahrscheinlichkeitsdichte. Es kann auch die \code{hist()} Funktion zur
Bestimmung des PSTHs verwendet werden. \sindex[term]{Feuerrate!Binningmethode}
Die bestimmte Feuerrate gilt f\"ur das gesamte Bin (\figref{binpsth}
B). Das so berechnete PSTH hat wiederum eine stufige Form, die von der
Wahl der Binweite anh\"angt. $r(t)$ ist also keine stetige
Funktion. Die Binweite bestimmt die zeitliche Aufl\"osung der
Absch\"atzung. \"Anderungen in der Feuerrate, die innerhalb eines Bins
vorkommen k\"onnen nicht aufgl\"ost werden. Mit der Wahl der Binweite
wird somit eine Annahme \"uber die relevante Zeitskala des Spiketrains
gemacht.
The binning method separates the time axis into regular bins of the
bin width $W$ and counts for each bin the number of observed action
potentials (\figref{binpsthfig} top). The resulting histogram is then
normalized with the bin width $W$ to yield the firing rate shown in
the bottom trace of figure \ref{binpsthfig}. The above sketched
process is equivalent to estimating the probability density. It is
possible to estimate the PSTH using the \code{hist()} method
\sindex[term]{Feuerrate!Binningmethode}
The estimated firing rate is valid for the total duration of each
bin. This leads to the step-like plot shown in
\figref{binpsthfig}. $r(t)$ is thus not a contiunous function in
time. The binwidth defines the temporal resolution of the firing rate
estimation Changes that happen within a bin cannot be resolved. Thus
chosing a bin width implies an assumption about the relevant
time-scale.
\pagebreak[4]
\begin{exercise}{binnedRate.m}{}
Implementiere die Absch\"atzung der Feuerrate mit der ``binning''
Methode. Plotte das PSTH.
Implement a function that estimates the firing rate using the
``binning method''. The method should take the spike-times as an
input argument and returns the firing rate. Plot the PSTH.
\end{exercise}
\subsubsection{Faltungsmethode}
\subsubsection{Convolution method --- Kernel density estimation}
\begin{figure}[tp]
\includegraphics[width=\columnwidth]{convmethod}
\titlecaption{Bestimmung des PSTH mit der Faltungsmethode.}{Der
gleiche Spiketrain wie in \figref{instrate}. Bei der Verfaltung
des Spiketrains mit einem Faltungskern wird jeder Spike durch den
Faltungskern ersetzt (oben). Bei korrekter Normierung des
Kerns ergibt sich die Feuerrate direkt aus der \"Uberlagerung der
Kerne.}\label{convrate}
\titlecaption{Estimating the firing rate using the convolution
method.}{The same spiketrain as in \figref{instratefig}. The
convolution of the spiketrain with a convolution kernel basically
replaces each spike event with the kernel (top). A Gaussian kernel
is used here, but other kernels are also possible. If the kernel
is properly normalized the firing rate results directly form the
superposition of the kernels.}\label{convratefig}
\end{figure}
Bei der Faltungsmethode werden die harten Kanten der Bins der
Binning-Methode vermieden. Der Spiketrain wird mit einem Kern
verfaltet, d.h. jedes Aktionspotential wird durch den Kern ersetzt.
Zur Berechnung wird die Aktionspotentialfolge zun\"achst
``bin\"ar'' dargestellt. Dabei wird ein Spiketrain als
(Zeit-)Vektor dargestellt, in welchem die Zeitpunkte der
Aktionspotentiale als 1 notiert werden. Alle anderen Elemente des
Vektors sind 0. Anschlie{\ss}end wir dieser bin\"are Spiketrain mit
einem Gau{\ss}-Kern bestimmter Breite verfaltet:
\[r(t) = \int_{-\infty}^{\infty} \omega(\tau) \, \rho(t-\tau) \, {\rm d}\tau \; , \]
wobei $\omega(\tau)$ der Filterkern und $\rho(t)$ die bin\"are Antwort
ist. Bildlich geprochen wird jede 1 in $\rho(t)$ durch den Filterkern
ersetzt (Abbildung \ref{convrate} A). Wenn der Kern richtig normiert
wurde (Integral gleich Eins), ergibt sich die Feuerrate direkt aus der
\"Uberlagerung der Kerne (Abb. \ref{convrate} B). \sindex[term]{Feuerrate!Faltungsmethode}
Die Faltungsmethode f\"uhrt, anders als die anderen Methoden, zu einer
stetigen Funktion was insbesondere f\"ur spektrale Analysen von
Vorteil sein kann. Die Wahl der Kernbreite bestimmt, \"ahnlich zur
Binweite, die zeitliche Aufl\"osung von $r(t)$. Die Breite des Kerns
macht also auch wieder eine Annahme \"uber die relevante Zeitskala des
Spiketrains.
With the convolution method we avoid the sharp edges of the binning
method. The spiketrain is convolved with a \enterm{convolution
kernel}. Technically speaking we need to first create a binary
representation of the spike train. This binary representation is a
series of zeros and ones in which the ones denote the spike. Then this binary vector is convolved with a kernel of a certain width:
\[r(t) = \int_{-\infty}^{\infty} \omega(\tau) \, \rho(t-\tau) \, {\rm d}\tau \; , \]
where $\omega(\tau)$ represents the kernel and $\rho(t)$ the binary
representation of the response. The process of convolution can be
imagined as replacing each event of the spiketrain with the kernel
(figure \ref{convratefig} top). The superposition of the replaced
kernels is then the firing rate (if the kerel is correctly normalized
to an integral of one, figure \ref{convreffig}
bottom). \sindex[term]{Feuerrate!Faltungsmethode}
In contrast to the other methods the convolution methods leads to a
continuous function which is often desirable (in particular when
applying methods in the frequency domain). The choice of the kernel
width defines, similar to the bin width of the binning method, the
temporal resolution of the method and thus makes assumptions about the
relevate time-scale.
\pagebreak[4]
\begin{exercise}{convolutionRate.m}{}
Verwende die Faltungsmethode um die Feuerrate zu bestimmen. Plotte
das Ergebnis.
Implement the function that estimates the firing rate using the
convolution method. The method takes the spiketrain, the temporal
resolution of the recording (as the stepsize $dt$, in seconds) and
the width of the kernel (the standard deviation $\sigma$ of the
Gaussian kernel, in seconds) as input arguments. It returns the
firing rate. Plot the result.
\end{exercise}
\section{Spike-triggered Average}
Die graphischer Darstellung der Feuerrate allein reicht nicht aus um
den Zusammenhang zwischen neuronaler Antwort und einem Stimulus zu
analysieren. Eine Methode um mehr \"uber diesen Zusammenhang zu
erfahren, ist der \enterm{spike-triggered average}
(\enterm[STA|see{spike-triggered average}]{STA}). Der STA
The graphical representation of the neuronal activity alone is not
sufficient tot investigate the relation between the neuronal response
and a stimulus. One method to do this is the \enterm[STA|see
spike-triggered average]{spike-triggered average}. The STA
\begin{equation}
STA(\tau) = \langle s(t - \tau) \rangle = \frac{1}{N} \sum_{i=1}^{N} s(t_i - \tau)
\end{equation}
der $N$ Aktionspotentiale zu den Zeiten $t_i$ in Anwort auf den
Stimulus $s(t)$ ist der mittlere Stimulus, der zu einem
Aktionspotential in der neuronalen Antwort f\"uhrt.
Der STA l\"a{\ss}t sich relativ einfach berechnen, indem aus dem
Stimulus f\"ur jeden beobachteten Spike ein entsprechender Abschnitt
ausgeschnitten wird und diese dann gemittelt werde (\figref{stafig}).
of $N$ action potentials observed at the times $t_i$ in response to
the stimulus $s(t)$ is the average stimulus that led to a spike in the
neuron. The STA can be easily extracted by cutting snippets out of
$s(t)$ that surround the times of the spikes. The resulting stimulus
snippets are then averaged (\figref{stafig}).
\begin{figure}[t]
\includegraphics[width=\columnwidth]{sta}
\titlecaption{Spike-triggered Average eines P-Typ Elektrorezeptors
und Stimulusrekonstruktion.}{Der STA (links): der Rezeptor
wurde mit einem ``white-noise'' Stimulus getrieben. Zeitpunkt 0
ist der Zeitpunkt des beobachteten Aktionspotentials. Die Kurve
ergibt sich aus dem gemittelten Stimulus je 50\,ms vor und nach
einem Aktionspotential. Stimulusrekonstruktion mittels
STA (rechts). Die Zellantwort wird mit dem STA gefaltet um eine
Rekonstruktion des Stimulus zu erhalten.}\label{stafig}
\titlecaption{Spike-triggered average of a P-type electroreceptors
and the stimulus reconstruction.}{The neuron was driven by a
``white-noise'' stimulus (blue curve, right panel). The STA (left)
is the average stimulus that surrounds the times of the recorded
action potentials (40\,ms before and 20\,ms after the
spike). Using the STA as a convolution kernel for convolving the
spiketrain we can reconstruct the stimulus from the neuronal
response. In this way we can get an impression of the stimulus
features that are encoded in the neuronal response (right
panel).}\label{stafig}
\end{figure}
Aus dem STA k\"onnen verschiedene Informationen \"uber den
Zusammenhang zwischen Stimulus und neuronaler Antwort gewonnen
werden. Die Breite des STA repr\"asentiert die zeitliche Pr\"azision,
mit der Stimulus und Antwort zusammenh\"angen und wie weit das Neuron
zeitlich integriert. Die Amplitude des STA (gegeben in der gleichen
Einheit wie der Stimulus) deutet auf die Empfindlichkeit des Neurons
bez\"uglich des Stimulus hin. Eine hohe Amplitude besagt, dass es
einen starken Stimulus ben\"otigt, um ein Aktionspotential
hervorzurufen. Aus dem zeitlichen Versatz des STA kann die Zeit
abgelesen werden, die das System braucht, um auf den Stimulus zu
antworten.
Der STA kann auch dazu benutzt werden, aus den Antworten der Zelle den
Stimulus zu rekonstruieren (\figref{stafig} B). Bei der
\determ[invertierte Rekonstruktion]{invertierten Rekonstruktion} wird
die Zellantwort mit dem STA verfaltet.
From the STA we can extract several pieces of information about the
relation of stimulus and response. The width of the STA represents the
temporal precision with which the neuron encodes the stimulus waveform
and in which temporal window the neuron integrates the (sensory)
input. The amplitude of the STA tells something about the sensitivity
of the neuron. The STA is given in the same units as the stimulus and
a small amplitude indicates that the neuron needs only a small
stimulus amplitude to create a spike, a large amplitude on the
contrary suggests the opposite. The temporal delay between the STA and
the time of the spike is a consequence of the time the system (neuron)
needs to process the stimulus.
We can further use the STA to do a \enterm{reverse reconstruction} and
estimate the stimulus from the neuronal response (\figref{stafig},
right). For this, the spiketrain is convolved with the STA as a
kernel.
\begin{exercise}{spikeTriggeredAverage.m}{}
Implementiere eine Funktion, die den STA ermittelt. Verwende dazu
den Datensatz \file{sta\_data.mat}. Die Funktion sollte folgende
R\"uckgabewerte haben:
Implement a function that calculates the STA. Use the dataset
\file{sta\_data.mat}. The function expects the spike train, the
stimulus and the temporal resolution of the recording as input
arguments and should return the following information:
\vspace{-1ex}
\begin{itemize}
\setlength{\itemsep}{0ex}
\item den Spike-Triggered-Average.
\item die Standardabweichung der individuellen STAs.
\item die Anzahl Aktionspotentiale, die zur Berechnung des STA verwendet wurden.
\vspace{-2ex}
\item the spike-triggered average.
\item the standard deviation of the STA across the individual snippets.
\item The number of action potentials used to estimate the STA.
\end{itemize}
\end{exercise}
\begin{exercise}{reconstructStimulus.m}{}
Rekonstruiere den Stimulus mithilfe des STA und der Spike
Zeiten. Die Funktion soll Vektor als R\"uckgabewert haben, der
genauso gro{\ss} ist wie der Originalstimulus aus der Datei
Do the reverse reconstruction using the STA and the spike times. The
function should return the estimated stimulus in a vector that has
the same size as the original stimulus contained in file
\file{sta\_data.mat}.
\end{exercise}

515
pointprocesses/lecture/pointprocesses_de.tex Normal file
View File

@ -0,0 +1,515 @@
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\chapter{Analyse von Spiketrains}
\selectlanguage{ngerman}
\determ[Aktionspotential]{Aktionspotentiale} (\enterm{spikes}) sind die Tr\"ager der
Information in Nervensystemen. Dabei ist in erster Linie nur der
Zeitpunkt des Auftretens eines Aktionspotentials von Bedeutung. Die
genaue Form des Aktionspotentials spielt keine oder nur eine
untergeordnete Rolle.
Nach etwas Vorverarbeitung haben elektrophysiologische Messungen
deshalb Listen von Spikezeitpunkten als Ergebniss --- sogenannte
\enterm{spiketrains}. Diese Messungen k\"onnen wiederholt werden und
es ergeben sich mehrere \enterm{trials} von Spiketrains
(\figref{rasterexamplesfig}).
Spiketrains sind Zeitpunkte von Ereignissen --- den Aktionspotentialen
--- und deren Analyse f\"allt daher in das Gebiet der Statistik von
sogenannten \determ[Punktprozess]{Punktprozessen}.
\begin{figure}[ht]
\includegraphics[width=1\textwidth]{rasterexamples}
\titlecaption{\label{rasterexamplesfig}Raster-Plot.}{Raster-Plot von
jeweils 10 Realisierungen eines station\"arenen Punktprozesses
(homogener Poisson Prozess mit Rate $\lambda=20$\;Hz, links) und
eines nicht-station\"aren Punktprozesses (perfect
integrate-and-fire Neuron getrieben mit Ohrnstein-Uhlenbeck
Rauschen mit Zeitkonstante $\tau=100$\,ms, rechts). Jeder
vertikale Strich markiert den Zeitpunkt eines Ereignisses.
Jede Zeile zeigt die Ereignisse eines trials.}
\end{figure}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Punktprozesse}
Ein zeitlicher Punktprozess (\enterm{point process}) ist ein
stochastischer Prozess, der eine Abfolge von Ereignissen zu den Zeiten
$\{t_i\}$, $t_i \in \reZ$, generiert.
\begin{ibox}{Beispiele von Punktprozessen}
Jeder Punktprozess wird durch einen sich in der Zeit kontinuierlich
entwickelnden Prozess generiert. Wann immer dieser Prozess eine
Schwelle \"uberschreitet wird ein Ereigniss des Punktprozesses
erzeugt. Zum Beispiel:
\begin{itemize}
\item Aktionspotentiale/Herzschlag: wird durch die Dynamik des
Membranpotentials eines Neurons/Herzzelle erzeugt.
\item Erdbeben: wird durch die Dynamik des Druckes zwischen
tektonischen Platten auf beiden Seiten einer geologischen Verwerfung
erzeugt.
\item Zeitpunkt eines Grillen/Frosch/Vogelgesangs: wird durch die
Dynamik des Nervensystems und des Muskelapparates erzeugt.
\end{itemize}
\end{ibox}
\begin{figure}[t]
\texpicture{pointprocessscetch}
\titlecaption{\label{pointprocessscetchfig} Statistik von
Punktprozessen.}{Ein Punktprozess ist eine Abfolge von
Zeitpunkten $t_i$ die auch durch die Intervalle $T_i=t_{i+1}-t_i$
oder die Anzahl der Ereignisse $n_i$ beschrieben werden kann. }
\end{figure}
F\"ur die Neurowissenschaften ist die Statistik der Punktprozesse
besonders wichtig, da die Zeitpunkte der Aktionspotentiale als
zeitlicher Punktprozess betrachtet werden k\"onnen und entscheidend
f\"ur die Informations\"ubertragung sind.
Bei Punktprozessen k\"onnen wir die Zeitpunkte $t_i$ ihres Auftretens,
die Intervalle zwischen diesen Zeitpunkten $T_i=t_{i+1}-t_i$, sowie
die Anzahl der Ereignisse $n_i$ bis zu einer bestimmten Zeit betrachten
(\figref{pointprocessscetchfig}).
Zwei Punktprozesse mit verschiedenen Eigenschaften sind in
\figref{rasterexamplesfig} als Rasterplot dargestellt, bei dem die
Zeitpunkte der Ereignisse durch senkrechte Striche markiert werden.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Intervallstatistik}
Die Intervalle $T_i=t_{i+1}-t_i$ zwischen aufeinanderfolgenden
Ereignissen sind reelle, positive Zahlen. Bei Aktionspotentialen
heisen die Intervalle auch \determ{Interspikeintervalle}
(\enterm{interspike intervals}). Deren Statistik kann mit den
\"ublichen Gr\"o{\ss}en beschrieben werden.
\begin{figure}[t]
\includegraphics[width=0.96\textwidth]{isihexamples}\vspace{-2ex}
\titlecaption{\label{isihexamplesfig}Interspikeintervall Histogramme}{der in
\figref{rasterexamplesfig} gezeigten Spikes.}
\end{figure}
\begin{exercise}{isis.m}{}
Schreibe eine Funktion \code{isis()}, die aus mehreren trials von Spiketrains die
Interspikeintervalle bestimmt und diese in einem Vektor
zur\"uckgibt. Jeder trial der Spiketrains ist ein Vektor mit den
Spikezeiten gegeben in Sekunden als Element in einem \codeterm{cell-array}.
\end{exercise}
\subsection{Intervallstatistik erster Ordnung}
\begin{itemize}
\item Wahrscheinlichkeitsdichte $p(T)$ der Intervalle $T$
(\figref{isihexamplesfig}). Normiert auf $\int_0^{\infty} p(T) \; dT
= 1$.
\item Mittleres Intervall: $\mu_{ISI} = \langle T \rangle =
\frac{1}{n}\sum\limits_{i=1}^n T_i$.
\item Standardabweichung der Intervalle: $\sigma_{ISI} = \sqrt{\langle (T - \langle T
\rangle)^2 \rangle}$\vspace{1ex}
\item \determ{Variationskoeffizient} (\enterm{coefficient of variation}): $CV_{ISI} =
\frac{\sigma_{ISI}}{\mu_{ISI}}$.
\item \determ{Diffusionskoeffizient} (\enterm{diffusion coefficient}): $D_{ISI} =
\frac{\sigma_{ISI}^2}{2\mu_{ISI}^3}$.
\end{itemize}
\begin{exercise}{isihist.m}{}
Schreibe eine Funktion \code{isiHist()}, die einen Vektor mit Interspikeintervallen
entgegennimmt und daraus ein normiertes Histogramm der Interspikeintervalle
berechnet.
\end{exercise}
\begin{exercise}{plotisihist.m}{}
Schreibe eine Funktion, die die Histogrammdaten der Funktion
\code{isiHist()} entgegennimmt, um das Histogramm zu plotten. Im
Plot sollen die Interspikeintervalle in Millisekunden aufgetragen
werden. Das Histogramm soll zus\"atzlich mit Mittelwert,
Standardabweichung und Variationskoeffizient der
Interspikeintervalle annotiert werden.
\end{exercise}
\subsection{Korrelationen der Intervalle}
In \enterm{return maps} werden die um das \enterm{lag} $k$ verz\"ogerten
Intervalle $T_{i+k}$ gegen die Intervalle $T_i$ geplottet. Dies macht
m\"ogliche Abh\"angigkeiten von aufeinanderfolgenden Intervallen
sichtbar.
\begin{figure}[t]
\includegraphics[width=1\textwidth]{returnmapexamples}
\includegraphics[width=1\textwidth]{serialcorrexamples}
\titlecaption{\label{returnmapfig}Interspikeintervall return maps und
serielle Korrelationen}{zwischen aufeinander folgenden Intervallen
im Abstand des Lags $k$.}
\end{figure}
Solche Ab\"angigkeiten werden durch die \determ{serielle
Korrelationen} (\enterm{serial correlations}) der Intervalle
quantifiziert. Das ist der \determ{Korrelationskoeffizient} zwischen
aufeinander folgenden Intervallen getrennt durch 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) \]
\"Ublicherweise wird die Korrelation $\rho_k$ gegen den Lag $k$
aufgetragen (\figref{returnmapfig}). $\rho_0=1$ (Korrelation jedes
Intervalls mit sich selber).
\begin{exercise}{isiserialcorr.m}{}
Schreibe eine Funktion \code{isiserialcorr()}, die einen Vektor mit Interspikeintervallen
entgegennimmt und daraus die seriellen Korrelationen berechnet und plottet.
\pagebreak[4]
\end{exercise}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Z\"ahlstatistik}
% \begin{figure}[t]
% \includegraphics[width=0.48\textwidth]{poissoncounthist100hz10ms}\hfill
% \includegraphics[width=0.48\textwidth]{poissoncounthist100hz100ms}
% \titlecaption{\label{countstatsfig}Count Statistik.}{}
% \end{figure}
Die Anzahl der Ereignisse $n_i$ in Zeifenstern $i$ der
L\"ange $W$ ergeben ganzzahlige, positive Zufallsvariablen, die meist
durch folgende Sch\"atzer charakterisiert werden:
\begin{itemize}
\item Histogramm der counts $n_i$.
\item Mittlere Anzahl von Ereignissen: $\mu_N = \langle n \rangle$.
\item Varianz der Anzahl: $\sigma_n^2 = \langle (n - \langle n \rangle)^2 \rangle$.
\item \determ{Fano Faktor} (Varianz geteilt durch Mittelwert): $F = \frac{\sigma_n^2}{\mu_n}$.
\end{itemize}
Insbesondere ist die mittlere Rate der Ereignisse $r$ (Spikes pro
Zeit, \determ{Feuerrate}) gemessen in Hertz \sindex[term]{Feuerrate!mittlere Rate}
\begin{equation}
\label{firingrate}
r = \frac{\langle n \rangle}{W} \; .
\end{equation}
% \begin{figure}[t]
% \begin{minipage}[t]{0.49\textwidth}
% Poisson process $\lambda=100$\,Hz:\\
% \includegraphics[width=1\textwidth]{poissonfano100hz}
% \end{minipage}
% \hfill
% \begin{minipage}[t]{0.49\textwidth}
% LIF $I=10$, $\tau_{adapt}=100$\,ms:\\
% \includegraphics[width=1\textwidth]{lifadaptfano10-100ms}
% \end{minipage}
% \titlecaption{\label{fanofig}Fano factor.}{}
% \end{figure}
\begin{exercise}{counthist.m}{}
Schreibe eine Funktion \code{counthist()}, die aus mehreren trials
von Spiketrains die Verteilung der Anzahl der Spikes in Fenstern
einer der Funktion \"ubergegebenen Breite bestimmt, das Histogramm
plottet und zur\"uckgibt. Jeder trial der Spiketrains ist ein Vektor
mit den Spikezeiten gegeben in Sekunden als Element in einem
\codeterm{cell-array}.
\end{exercise}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Homogener Poisson Prozess}
F\"ur kontinuierliche Me{\ss}gr\"o{\ss}en ist die Normalverteilung
u.a. wegen dem Zentralen Grenzwertsatz die Standardverteilung. Eine
\"ahnliche Rolle spielt bei Punktprozessen der \determ{Poisson
Prozess}.
Beim \determ[Poisson Prozess!homogener]{homogenen Poisson Prozess}
treten Ereignisse mit einer festen Rate $\lambda=\text{const.}$ auf
und sind unabh\"angig von der Zeit $t$ und unabh\"angig von den
Zeitpunkten fr\"uherer Ereignisse (\figref{hompoissonfig}). Die
Wahrscheinlichkeit zu irgendeiner Zeit ein Ereigniss in einem kleinen
Zeitfenster der Breite $\Delta t$ zu bekommen ist
\begin{equation}
\label{hompoissonprob}
P = \lambda \cdot \Delta t \; .
\end{equation}
Beim \determ[Poisson Prozess!inhomogener]{inhomogenen Poisson Prozess}
h\"angt die Rate $\lambda$ von der Zeit ab: $\lambda = \lambda(t)$.
\begin{exercise}{poissonspikes.m}{}
Schreibe eine Funktion \code{poissonspikes()}, die die Spikezeiten
eines homogenen Poisson-Prozesses mit gegebener Rate in Hertz f\"ur
eine Anzahl von trials gegebener maximaler L\"ange in Sekunden in
einem \codeterm{cell-array} zur\"uckgibt. Benutze \eqnref{hompoissonprob}
um die Spikezeiten zu bestimmen.
\end{exercise}
\begin{figure}[t]
\includegraphics[width=1\textwidth]{poissonraster100hz}
\titlecaption{\label{hompoissonfig}Rasterplot von Spikes eines homogenen
Poisson Prozesses mit $\lambda=100$\,Hz.}{}
\end{figure}
\begin{figure}[t]
\includegraphics[width=0.45\textwidth]{poissonisihexp20hz}\hfill
\includegraphics[width=0.45\textwidth]{poissonisihexp100hz}
\titlecaption{\label{hompoissonisihfig}Interspikeintervallverteilungen
zweier Poissonprozesse.}{}
\end{figure}
Der homogene Poissonprozess hat folgende Eigenschaften:
\begin{itemize}
\item Die Intervalle $T$ sind exponentiell verteilt (\figref{hompoissonisihfig}):
\begin{equation}
\label{poissonintervals}
p(T) = \lambda e^{-\lambda T} \; .
\end{equation}
\item Das mittlere Intervall ist $\mu_{ISI} = \frac{1}{\lambda}$ .
\item Die Varianz der Intervalle ist $\sigma_{ISI}^2 = \frac{1}{\lambda^2}$ .
\item Der Variationskoeffizient ist also immer $CV_{ISI} = 1$ .
\item Die \determ[serielle Korrelationen]{seriellen Korrelationen}
$\rho_k =0$ f\"ur $k>0$, da das Auftreten der Ereignisse
unabh\"angig von der Vorgeschichte ist. Ein solcher Prozess wird
auch \determ{Erneuerungsprozess} genannt (\enterm{renewal process}).
\item Die Anzahl der Ereignisse $k$ innerhalb eines Fensters der
L\"ange W ist \determ[Poisson-Verteilung]{Poissonverteilt}:
\[ P(k) = \frac{(\lambda W)^ke^{\lambda W}}{k!} \]
(\figref{hompoissoncountfig})
\item Der \determ{Fano Faktor} ist immer $F=1$ .
\end{itemize}
\begin{exercise}{hompoissonspikes.m}{}
Schreibe eine Funktion \code{hompoissonspikes()}, die die Spikezeiten
eines homogenen Poisson-Prozesses mit gegebener Rate in Hertz f\"ur
eine Anzahl von trials gegebener maximaler L\"ange in Sekunden in
einem \codeterm{cell-array} zur\"uckgibt. Benutze die exponentiell-verteilten
Interspikeintervalle \eqnref{poissonintervals}, um die Spikezeiten zu erzeugen.
\end{exercise}
\begin{figure}[t]
\includegraphics[width=0.48\textwidth]{poissoncounthistdist100hz10ms}\hfill
\includegraphics[width=0.48\textwidth]{poissoncounthistdist100hz100ms}
\titlecaption{\label{hompoissoncountfig}Z\"ahlstatistik von Poisson Spikes.}{}
\end{figure}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Zeitabh\"angige Feuerraten}
Bisher haben wir station\"are Spiketrains betrachtet, deren Statistik
sich innerhalb der Analysezeit nicht ver\"andert (station\"are
Punktprozesse). Meistens jedoch \"andert sich die Statistik der
Spiketrains eines Neurons mit der Zeit. Z.B. kann ein sensorisches
Neuron auf einen Reiz hin mit einer erh\"ohten Feuerrate antworten
(nichtstation\"arer Punktprozess).
Wie die mittlere Anzahl der Spikes sich mit der Zeit ver\"andert, die
\determ{Feuerrate} $r(t)$, ist die wichtigste Gr\"o{\ss}e bei
nicht-station\"aren Spiketrains. Die Einheit der Feuerrate ist Hertz,
also Anzahl Aktionspotentiale pro Sekunde. Es gibt verschiedene
Methoden diese zu bestimmen. Drei solcher Methoden sind in Abbildung
\ref{psthfig} dargestellt. Alle Methoden haben ihre Berechtigung und
ihre Vor- und Nachteile. Im folgenden werden die drei Methoden aus
Abbildung \ref{psthfig} n\"aher erl\"autert.
\begin{figure}[tp]
\includegraphics[width=\columnwidth]{firingrates}
\titlecaption{Bestimmung der zeitabh\"angigen
Feuerrate.}{\textbf{A)} Rasterplot eines Spiketrains. \textbf{B)}
Feurerrate aus der instantanen Feuerrate bestimmt. \textbf{C)}
klassisches PSTH mit der Binning Methode. \textbf{D)} Feuerrate
durch Faltung mit einem Gauss Kern bestimmt.}\label{psthfig}
\end{figure}
\subsection{Instantane Feuerrate}
\begin{figure}[tp]
\includegraphics[width=\columnwidth]{isimethod}
\titlecaption{Instantane Feuerrate.}{Skizze eines Spiketrains
(oben). Die Pfeile zwischen aufeinanderfolgenden
Aktionspotentialen mit den Zahlen in Millisekunden illustrieren
die Interspikeintervalle. Der Kehrwert des Interspikeintervalle
ergibt die instantane Feuerrate.}\label{instrate}
\end{figure}
Ein sehr einfacher Weg, die zeitabh\"angige Feuerrate zu bestimmen ist
die sogenannte \determ[Feuerrate!instantane]{instantane Feuerrate}
(\enterm[firing rate!instantaneous]{instantaneous firing rate}). Dabei
wird die Feuerrate aus dem Kehrwert der Interspikeintervalle, der Zeit
zwischen zwei aufeinander folgenden Aktionspotentialen
(\figref{instrate} A), bestimmt. Die abgesch\"atzte Feuerrate
(\figref{instrate} B) ist g\"ultig f\"ur das gesammte
Interspikeintervall. Diese Methode hat den Vorteil, dass sie sehr
einfach zu berechnen ist und keine Annahme \"uber eine relevante
Zeitskala (der Kodierung oder des Auslesemechanismus der
postsynaptischen Zelle) macht. $r(t)$ ist allerdings keine
kontinuierliche Funktion, die Spr\"unge in der Feuerrate k\"onnen
f\"ur manche Analysen nachteilig sein. Au{\ss}erdem wird die Feuerrate
nie gleich Null, auch wenn lange keine Aktionspotentiale generiert
wurden.
\begin{exercise}{instantaneousRate.m}{}
Implementiere die Absch\"atzung der Feuerrate auf Basis der
instantanen Feuerrate. Plotte die Feuerrate als Funktion der Zeit.
\end{exercise}
\subsection{Peri-Stimulus-Zeit-Histogramm}
W\"ahrend die Instantane Rate den Kehrwert der Zeit von einem bis zum
n\"achsten Aktionspotential misst, sch\"atzt das sogenannte
\determ{Peri-Stimulus-Zeit-Histogramm} (\enterm{peri stimulus time
histogram}, \determ[PSTH|see{Peri-Stimulus-Zeit-Histogramm}]{PSTH})
die Wahrscheinlichkeit ab, zu einem Zeitpunkt Aktionspotentiale
anzutreffen. Es wird versucht die mittlere Rate \eqnref{firingrate} im
Grenzwert kleiner Beobachtungszeiten abzusch\"atzen:
\begin{equation}
\label{psthrate}
r(t) = \lim_{W \to 0} \frac{\langle n \rangle}{W} \; ,
\end{equation}
wobei die Anzahl $n$ der Aktionspotentiale, die im Zeitintervall
$(t,t+W)$ aufgetreten sind, \"uber trials gemittelt wird. Eine solche
Rate enspricht der zeitabh\"angigen Rate $\lambda(t)$ des inhomogenen
Poisson-Prozesses.
Das PSTH \eqnref{psthrate} kann entweder \"uber die Binning-Methode
oder durch Verfaltung mit einem Kern bestimmt werden. Beiden Methoden
gemeinsam ist die Notwendigkeit der Wahl einer zus\"atzlichen Zeitskala,
die der Beobachtungszeit $W$ in \eqnref{psthrate} entspricht.
\subsubsection{Binning-Methode}
\begin{figure}[tp]
\includegraphics[width=\columnwidth]{binmethod}
\titlecaption{Bestimmung des PSTH mit der Binning Methode.}{Der
gleiche Spiketrain wie in \figref{instrate}. Die grauen Linien
markieren die Grenzen der Bins und die Zahlen geben die Anzahl der Spikes
in jedem Bin an (oben). Die Feuerrate ergibt sich aus dem
mit der Binbreite normierten Zeithistogramm (unten).}\label{binpsth}
\end{figure}
Bei der Binning-Methode wird die Zeitachse in gleichm\"aßige
Abschnitte (Bins) eingeteilt und die Anzahl Aktionspotentiale, die in
die jeweiligen Bins fallen, gez\"ahlt (\figref{binpsth} A). Um diese
Z\"ahlungen in die Feuerrate umzurechnen muss noch mit der Binweite
normiert werden. Das ist \"aquivalent zur Absch\"atzung einer
Wahrscheinlichkeitsdichte. Es kann auch die \code{hist()} Funktion zur
Bestimmung des PSTHs verwendet werden. \sindex[term]{Feuerrate!Binningmethode}
Die bestimmte Feuerrate gilt f\"ur das gesamte Bin (\figref{binpsth}
B). Das so berechnete PSTH hat wiederum eine stufige Form, die von der
Wahl der Binweite anh\"angt. $r(t)$ ist also keine stetige
Funktion. Die Binweite bestimmt die zeitliche Aufl\"osung der
Absch\"atzung. \"Anderungen in der Feuerrate, die innerhalb eines Bins
vorkommen k\"onnen nicht aufgl\"ost werden. Mit der Wahl der Binweite
wird somit eine Annahme \"uber die relevante Zeitskala des Spiketrains
gemacht.
\pagebreak[4]
\begin{exercise}{binnedRate.m}{}
Implementiere die Absch\"atzung der Feuerrate mit der ``binning''
Methode. Plotte das PSTH.
\end{exercise}
\subsubsection{Faltungsmethode}
\begin{figure}[tp]
\includegraphics[width=\columnwidth]{convmethod}
\titlecaption{Bestimmung des PSTH mit der Faltungsmethode.}{Der
gleiche Spiketrain wie in \figref{instrate}. Bei der Verfaltung
des Spiketrains mit einem Faltungskern wird jeder Spike durch den
Faltungskern ersetzt (oben). Bei korrekter Normierung des
Kerns ergibt sich die Feuerrate direkt aus der \"Uberlagerung der
Kerne.}\label{convrate}
\end{figure}
Bei der Faltungsmethode werden die harten Kanten der Bins der
Binning-Methode vermieden. Der Spiketrain wird mit einem Kern
verfaltet, d.h. jedes Aktionspotential wird durch den Kern ersetzt.
Zur Berechnung wird die Aktionspotentialfolge zun\"achst
``bin\"ar'' dargestellt. Dabei wird ein Spiketrain als
(Zeit-)Vektor dargestellt, in welchem die Zeitpunkte der
Aktionspotentiale als 1 notiert werden. Alle anderen Elemente des
Vektors sind 0. Anschlie{\ss}end wir dieser bin\"are Spiketrain mit
einem Gau{\ss}-Kern bestimmter Breite verfaltet:
\[r(t) = \int_{-\infty}^{\infty} \omega(\tau) \, \rho(t-\tau) \, {\rm d}\tau \; , \]
wobei $\omega(\tau)$ der Filterkern und $\rho(t)$ die bin\"are Antwort
ist. Bildlich geprochen wird jede 1 in $\rho(t)$ durch den Filterkern
ersetzt (Abbildung \ref{convrate} A). Wenn der Kern richtig normiert
wurde (Integral gleich Eins), ergibt sich die Feuerrate direkt aus der
\"Uberlagerung der Kerne (Abb. \ref{convrate} B). \sindex[term]{Feuerrate!Faltungsmethode}
Die Faltungsmethode f\"uhrt, anders als die anderen Methoden, zu einer
stetigen Funktion was insbesondere f\"ur spektrale Analysen von
Vorteil sein kann. Die Wahl der Kernbreite bestimmt, \"ahnlich zur
Binweite, die zeitliche Aufl\"osung von $r(t)$. Die Breite des Kerns
macht also auch wieder eine Annahme \"uber die relevante Zeitskala des
Spiketrains.
\pagebreak[4]
\begin{exercise}{convolutionRate.m}{}
Verwende die Faltungsmethode um die Feuerrate zu bestimmen. Plotte
das Ergebnis.
\end{exercise}
\section{Spike-triggered Average}
Die graphischer Darstellung der Feuerrate allein reicht nicht aus um
den Zusammenhang zwischen neuronaler Antwort und einem Stimulus zu
analysieren. Eine Methode um mehr \"uber diesen Zusammenhang zu
erfahren, ist der \enterm{spike-triggered average}
(\enterm[STA|see{spike-triggered average}]{STA}). Der STA
\begin{equation}
STA(\tau) = \langle s(t - \tau) \rangle = \frac{1}{N} \sum_{i=1}^{N} s(t_i - \tau)
\end{equation}
der $N$ Aktionspotentiale zu den Zeiten $t_i$ in Anwort auf den
Stimulus $s(t)$ ist der mittlere Stimulus, der zu einem
Aktionspotential in der neuronalen Antwort f\"uhrt.
Der STA l\"a{\ss}t sich relativ einfach berechnen, indem aus dem
Stimulus f\"ur jeden beobachteten Spike ein entsprechender Abschnitt
ausgeschnitten wird und diese dann gemittelt werde (\figref{stafig}).
\begin{figure}[t]
\includegraphics[width=\columnwidth]{sta}
\titlecaption{Spike-triggered Average eines P-Typ Elektrorezeptors
und Stimulusrekonstruktion.}{Der STA (links): der Rezeptor
wurde mit einem ``white-noise'' Stimulus getrieben. Zeitpunkt 0
ist der Zeitpunkt des beobachteten Aktionspotentials. Die Kurve
ergibt sich aus dem gemittelten Stimulus je 50\,ms vor und nach
einem Aktionspotential. Stimulusrekonstruktion mittels
STA (rechts). Die Zellantwort wird mit dem STA gefaltet um eine
Rekonstruktion des Stimulus zu erhalten.}\label{stafig}
\end{figure}
Aus dem STA k\"onnen verschiedene Informationen \"uber den
Zusammenhang zwischen Stimulus und neuronaler Antwort gewonnen
werden. Die Breite des STA repr\"asentiert die zeitliche Pr\"azision,
mit der Stimulus und Antwort zusammenh\"angen und wie weit das Neuron
zeitlich integriert. Die Amplitude des STA (gegeben in der gleichen
Einheit wie der Stimulus) deutet auf die Empfindlichkeit des Neurons
bez\"uglich des Stimulus hin. Eine hohe Amplitude besagt, dass es
einen starken Stimulus ben\"otigt, um ein Aktionspotential
hervorzurufen. Aus dem zeitlichen Versatz des STA kann die Zeit
abgelesen werden, die das System braucht, um auf den Stimulus zu
antworten.
Der STA kann auch dazu benutzt werden, aus den Antworten der Zelle den
Stimulus zu rekonstruieren (\figref{stafig} B). Bei der
\determ[invertierte Rekonstruktion]{invertierten Rekonstruktion} wird
die Zellantwort mit dem STA verfaltet.
\begin{exercise}{spikeTriggeredAverage.m}{}
Implementiere eine Funktion, die den STA ermittelt. Verwende dazu
den Datensatz \file{sta\_data.mat}. Die Funktion sollte folgende
R\"uckgabewerte haben:
\vspace{-1ex}
\begin{itemize}
\setlength{\itemsep}{0ex}
\item den Spike-Triggered-Average.
\item die Standardabweichung der individuellen STAs.
\item die Anzahl Aktionspotentiale, die zur Berechnung des STA verwendet wurden.
\vspace{-2ex}
\end{itemize}
\end{exercise}
\begin{exercise}{reconstructStimulus.m}{}
Rekonstruiere den Stimulus mithilfe des STA und der Spike
Zeiten. Die Funktion soll Vektor als R\"uckgabewert haben, der
genauso gro{\ss} ist wie der Originalstimulus aus der Datei
\file{sta\_data.mat}.
\end{exercise}
\selectlanguage{english}

4
pointprocesses/lecture/rasterexamples.py
View File

@ -31,7 +31,7 @@ def pifspikes(input, trials, dt, D=0.1) :
times = []
v = vreset
noise = np.sqrt(2.0*D)*np.random.randn(len(input))/np.sqrt(dt)
for k in xrange(len(noise)) :
for k in range(len(noise)) :
v += (input[k]+noise[k])*dt/tau
if v >= vthresh :
v = vreset
@ -55,7 +55,7 @@ rng = np.random.RandomState(54637281)
time = np.arange(0.0, duration, dt)
x = np.zeros(time.shape)+rate
n = rng.randn(len(time))*drate*tau/np.sqrt(dt)+rate
for k in xrange(1,len(x)) :
for k in range(1,len(x)) :
x[k] = x[k-1] + (n[k]-x[k-1])*dt/tau
x[x<0.0] = 0.0

6
pointprocesses/lecture/returnmapexamples.py
View File

@ -31,7 +31,7 @@ def pifspikes(input, trials, dt, D=0.1) :
times = []
v = vreset
noise = np.sqrt(2.0*D)*np.random.randn(len(input))/np.sqrt(dt)
for k in xrange(len(noise)) :
for k in range(len(noise)) :
v += (input[k]+noise[k])*dt/tau
if v >= vthresh :
v = vreset
@ -41,7 +41,7 @@ def pifspikes(input, trials, dt, D=0.1) :
def isis( spikes ) :
isi = []
for k in xrange(len(spikes)) :
for k in range(len(spikes)) :
isi.extend(np.diff(spikes[k]))
return np.array( isi )
@ -84,7 +84,7 @@ rng = np.random.RandomState(54637281)
time = np.arange(0.0, duration, dt)
x = np.zeros(time.shape)+rate
n = rng.randn(len(time))*drate*tau/np.sqrt(dt)+rate
for k in xrange(1,len(x)) :
for k in range(1,len(x)) :
x[k] = x[k-1] + (n[k]-x[k-1])*dt/tau
x[x<0.0] = 0.0

6
pointprocesses/lecture/serialcorrexamples.py
View File

@ -31,7 +31,7 @@ def pifspikes(input, trials, dt, D=0.1) :
times = []
v = vreset
noise = np.sqrt(2.0*D)*np.random.randn(len(input))/np.sqrt(dt)
for k in xrange(len(noise)) :
for k in range(len(noise)) :
v += (input[k]+noise[k])*dt/tau
if v >= vthresh :
v = vreset
@ -41,7 +41,7 @@ def pifspikes(input, trials, dt, D=0.1) :
def isis( spikes ) :
isi = []
for k in xrange(len(spikes)) :
for k in range(len(spikes)) :
isi.extend(np.diff(spikes[k]))
return np.array( isi )
@ -95,7 +95,7 @@ rng = np.random.RandomState(54637281)
time = np.arange(0.0, duration, dt)
x = np.zeros(time.shape)+rate
n = rng.randn(len(time))*drate*tau/np.sqrt(dt)+rate
for k in xrange(1,len(x)) :
for k in range(1,len(x)) :
x[k] = x[k-1] + (n[k]-x[k-1])*dt/tau
x[x<0.0] = 0.0

2
pointprocesses/lecture/sta.py
View File

@ -6,7 +6,7 @@ from IPython import embed
def plot_sta(times, stim, dt, t_min=-0.1, t_max=.1):
count = 0
sta = np.zeros((abs(t_min) + abs(t_max))/dt)
sta = np.zeros(int((abs(t_min) + abs(t_max))/dt))
time = np.arange(t_min, t_max, dt)
if len(stim.shape) > 1 and stim.shape[1] > 1:
stim = stim[:,1]

21
programming/exercises/boolean_logical_indexing.tex
View File

@ -62,9 +62,11 @@ following pattern: ``variables\_datatypes\_\{lastname\}.m''
\part Execute and explain: \verb+bitand(10, 8)+
\part Execute and explain: \verb+bitor(10, 8)+
\end{parts}
\item Implement the following Boolean expressions. Test using randomly selected integer values for \verb+x+ and \verb+y+.
\item Implement the following Boolean expressions. Test your
implementations using random integer
numbers for \verb+x+ and \verb+y+ (\verb+randi+).
\begin{parts}
\part The result should be \verb+true+ if \verb+x+ greater than \verb+y+ and the sum of \verb+x+ and \verb+y+ is not less than 100.
\part The result should be \verb+true+ if \verb+x+ is greater than \verb+y+ and the sum of \verb+x+ and \verb+y+ is not less than 100.
\part The result should be \verb+true+ if \verb+x+ and \verb+y+ are not equal zero or \verb+x+ and \verb+y+ are equal.
\end{parts}
\end{questions}
@ -82,7 +84,7 @@ matrices that match in certain criteria. This process is called
the following commands and explain.
\begin{parts}
\part \verb+x < 5+
\part \verb+x( x < 5) )+
\part \verb+x( (x < 5) )+
\part \verb+x( (y <= 2) )+
\part \verb+x( (x > 2) | (y < 8) )+
\part \verb+x( (x == 0) & (y == 0) )+
@ -94,6 +96,19 @@ matrices that match in certain criteria. This process is called
\verb+x >= 33 and x < 66+ to 1 and all \verb+x >= 66+ to 2.
\part Count the number of elements in each class using Boolean expressions (\verb+sum+ can be used to count the matches).
\end{parts}
\question Plotting a periodic signal.
\begin{parts}
\part Load the file ``signal.mat'' into the workspace (use
\verb+load('signal.mat')+ or use the UI for it). It contains two
variables \verb+signal+ and \verb+time+. The signal has a period of
0.2\,s
\part What is the size of the data?
\part What is the temporal resolution of the time axis?
\part Plot the full data. Make sure that the axes are properly labeled (script chapter 3, or the matlab documentation).
\part Use logical indexing to select the periods individually and plot them into the same plot (first period starts at time 0.0\,s, the next at 0.2\,s and so on).
\end{parts}
\end{questions}
\end{document}

47
programming/exercises/control_flow.tex
View File

@ -14,7 +14,7 @@
%%%%% text size %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\usepackage[left=20mm,right=20mm,top=25mm,bottom=25mm]{geometry}
\pagestyle{headandfoot}
\header{{\bfseries\large Exercise 4}}{{\bfseries\large Control Flow}}{{\bfseries\large 30. October, 2017}}
\header{{\bfseries\large Exercise 5}}{{\bfseries\large Control Flow}}{{\bfseries\large 30. October, 2018}}
\firstpagefooter{Dr. Jan Grewe}{Phone: 29 74588}{Email:
jan.grewe@uni-tuebingen.de}
\runningfooter{}{\thepage}{}
@ -33,39 +33,44 @@
\begin{center}
\textbf{\Large Introduction to scientific computing}\\[1ex]
{\large Jan Grewe, Jan Benda}\\[-3ex]
Abteilung Neuroethologie \hfill --- \hfill Institut f\"ur Neurobiologie \hfill --- \hfill \includegraphics[width=0.28\textwidth]{UT_WBMW_Black_RGB} \\
Neuroethologie \hfill --- \hfill Institute for Neurobiology \hfill --- \hfill \includegraphics[width=0.28\textwidth]{UT_WBMW_Black_RGB} \\
\end{center}
The exercises are meant for self-monitoring, revision of the lecture
topic. You should try to solve them on your own. Your solution should
be submitted as a single script (m-file) in the Ilias system. Each
task should be solved in its own ``cell''. Each cell must be
executable on its own. The file should be named according to the following pattern:
``variables\_datatypes\_\{lastname\}.m'' benannt werden
The exercises are meant for self-monitoring and revision of the
lecture. You should try to solve them on your own. Your solution
should be submitted as a single script (m-file) in the Ilias
system. Each task should be solved in its own ``cell''. Each cell must
be executable on its own. The file should be named according to the
following pattern: ``variables\_datatypes\_\{lastname\}.m''
(e.g. variables\_datentypes\_mueller.m).
\begin{questions}
\question Implement \code{for} loops in which the \emph{running variable}:
\begin{parts}
\part ... assumes values from 0 to 10. Print the value of the running variable for each iteration step.
\begin{solution}
for i = 1:10
disp(i);
end;
\end{solution}
\part ... assumes values from 10 to 0. Print out the value for each iteration.
\part ... assumes values from 0 to 1 in steps of 0.1. Print out the value for each iteration.
\end{parts}
\question Indexing in vectros:
\question Indexing in vectors:
\begin{parts}
\part Define a vector \code{x} that contains values from 1:100.
\part Define a vector \code{x} that contains values ranging from 101 to 200.
\part Use a \code{for} loop to print each element of \code{x}. Use the running variable for the indexing.
\begin{solution}
\code{for i = 1:length(x); disp(x(i)); end}
\end{solution}
\part Do the same without using the running variable for indexing.
\part Do the same without use the running variable directly (not for indexing).
\begin{solution}
\code{for i : x; disp(i); end;}
\end{solution}
\end{parts}
\question Create a vector \verb+x+ that contains 50 random numbers in the range 0 - 10.
\question Create a vector \verb+x+ that contains 50 random numbers in the value range 0 - 10.
\begin{parts}
\part Use a loop to calculate the arithmetic mean. The mean is defined as:
$\overline{x}=\frac{1}{n}\sum\limits_{i=0}^{n}x_i $.
@ -74,20 +79,20 @@ executable on its own. The file should be named according to the following patte
\part Search the MATLAB help for functions that do the calculations for you :-).
\end{parts}
\question Implement a\code{while} loop
\question Implement a \code{while} loop
\begin{parts}
\part ... that iterates 100 times. Print out the current iteration number.
\part ... that iterates 100 times. Print out the current iteration number.
\part ... iterates endlessly. You can interrupt the execution with the command \code{Strg + C}.
\end{parts}
\question Use an endless \code{while} loop to individually print out the elements of a vector that contains 10 elements.
\question Use an endless \code{while} loop to individually print out the elements of a vector that contains 10 elements. Stop the execution after all elements have been printed.
\question Use the endless \verb+while+ loop to draw random numbers
(\verb+randn+) until you it is larger than 1.33.
(\verb+randn+) until you hit one that is larger than 1.33.
\begin{parts}
\part Count the number of required iterations.
\part Use a \code{for} loop to run the previous test 1000 times. Remeber all counts and calculate the mean of it.
\part Create a plot that schows for each of the 1000 trials, how often you had to draw.
\part Use a \code{for} loop to run the previous test 1000 times. Remember all counts and calculate the mean of it.
\part Create a plot that shows for each of the 1000 trials, how often you had to draw.
\part Play around with the threshold. What happens?
\end{parts}
@ -100,19 +105,19 @@ executable on its own. The file should be named according to the following patte
\end{parts}
\question Test the random number generator! To do so count the
number of elements that fall in the classes defined by the edges
number of elements that fall into the classes defined by the edges
[0.0, 0.2, 0.4, 0.6, 0.8, 1.0]. Store the results in a vector. Use a
loop to draw 1000 random numbers with \verb+rand()+ (see help). What
would be the expectation? What is the result?
\question String parsing: It is quite common to use the names of datasets to encode the conditions under which they were recorded. To find out the required information we have to \emph{parse} he name.
\question String parsing: It is quite common to use the names of datasets to encode the conditions under which they were recorded. To find out the required information we have to \emph{parse} the name.
\begin{parts}
\part Create a variable
\verb+filename = '2015-10-12_100Hz_1.25V.dat'+. Obviously, the
underscore was used as a delimiter.
\part Use a \verb+for+ loop to find the positions of the underscores and store these in a vector.
\part Use a second loop to iterate over the previously filled `position' vector and use this information to
cut filename appropriately.
cut \verb+filename+ appropriately.
\part Print out the individual parts.
\end{parts}

12
programming/exercises/scripts_functions.tex
View File

@ -16,7 +16,7 @@
%%%%% text size %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\usepackage[left=20mm,right=20mm,top=25mm,bottom=25mm]{geometry}
\pagestyle{headandfoot}
\header{{\bfseries\large \"Ubung 5}}{{\bfseries\large Scripts and functions}}{{\bfseries\large 07. November, 2017}}
\header{{\bfseries\large Exercise 6}}{{\bfseries\large Scripts and functions}}{{\bfseries\large 06. November, 2018}}
\firstpagefooter{Prof. Jan Benda}{Phone: 29 74 573}{Email:
jan.benda@uni-tuebingen.de}
\runningfooter{}{\thepage}{}
@ -46,7 +46,6 @@
aboveskip=10pt
}
\newcommand{\code}[1]{\texttt{#1}}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
@ -60,7 +59,7 @@
\end{center}
The exercises are meant for self-monitoring and revision of the
lecture topic. You should try to solve them on your own. In contrast
lecture. You should try to solve them on your own. In contrast
to previous exercises, the solutions can not be saved in a single file
but each question needs an individual file. Combine the files into a
single zip archive and submit it via ILIAS. Name the archive according
@ -108,7 +107,7 @@ to the pattern: ``scripts\_functions\_\{surname\}.zip''.
\part{}
Improve the function that it takes the duration of the time axis,
the amplitude and the frequency as input arguments. The
the amplitude, and the frequency as input arguments. The
calculation should use a temporal stepsize that is 0.01 of the
frequency.
\begin{solution}
@ -144,7 +143,8 @@ to the pattern: ``scripts\_functions\_\{surname\}.zip''.
\question Random Walk. In a 1-D random walk an \emph{agent} walks
randomly either in the one ($+1$) or the other ($-1$)
direction. With each simulation step one direction is chosen and the
position is updated accordingly.
position is updated accordingly. \textbf{There are some differences
to the solution discussed in the lecture!}
\begin{parts}
\part{}
@ -172,7 +172,7 @@ to the pattern: ``scripts\_functions\_\{surname\}.zip''.
Now we want to know how the probability of $p_{+1}$ (the
probability to walk into the $+1$ direction) impacts the random
walk. Vary $p_{+1}$ in the range $0.5 \le p_{+1} < 0.8$. Do 10
random walks for four probabilities (apply the same thresholds for
random walks for the four probabilities (apply the same thresholds for
stopping the simulations as before).
\begin{solution}
\lstinputlisting{randomwalkscriptc.m}

2
programming/exercises/variables_types.tex
View File

@ -1,4 +1,4 @@
\documentclass[12pt,a4paper,pdftex]{exam}
\documentclass[12pt,a4paper,pdftex, answers]{exam}
\usepackage[german]{babel}
\usepackage{natbib}

24
programming/lecture/programming.tex
View File

@ -40,14 +40,12 @@ variable.
\begin{figure}
\centering
\begin{subfigure}{.5\textwidth}
\includegraphics[width=0.8\textwidth]{variable}
\label{variable:a}
\includegraphics[width=0.8\textwidth]{variable}\label{variable:a}
\end{subfigure}%
\begin{subfigure}{.5\textwidth}
\includegraphics[width=.8\textwidth]{variableB}
\label{variable:b}
\includegraphics[width=.8\textwidth]{variableB}\label{variable:b}
\end{subfigure}
\titlecaption{Variables} point to a memory
\titlecaption{Variables}{ point to a memory
address. They further are described by their name and
data type. The variable's value is stored as a pattern of binary
values (0 or 1). When reading the variable this pattern is
@ -630,7 +628,7 @@ matrix). The function \code{cat()} allows to concatenate n-dimensional
matrices.
To request the length of a vector we used the function
\code{length()}. This function is \tetbf{not} suited to request
\code{length()}. This function is \textbf{not} suited to request
information about the size of a matrix. As mentioned above,
\code{length()} would return the length of the largest dimension. The
function \code{size()} however, returns the length in each dimension
@ -1031,7 +1029,7 @@ segment of data of a certain time span (the stimulus was on,
\varcode{x} represent the measured data at the times in
\varcode{t}.
\item Use logical indexing to select those values that have been
recorded in the time span form 5--6\,s.
recorded in the time span from 5--6\,s.
\end{itemize}
\end{exercise}
@ -1095,7 +1093,7 @@ ans =
Generally, a program is executed line by line from top to
bottom. Sometimes this is not the desired behavior, or the other way
round, it is needed to skip certain parts or execute others
repeatedly. High-level programming languages like \matlab{} offer
repeatedly. High-level programming languages such as \matlab{} offer
statements that allow to manipulate the control flow. There are two
major classes of such statements:
@ -1120,7 +1118,7 @@ x =
120
\end{lstlisting}
This kind of program is fine but it is rather repetitive. The only
This kind of program solves the taks but it is rather repetitive. The only
thing that changes is the increasing factor. The repetition of such
very similar lines of code is bad programming style. This is not only
a matter of taste but there are severe drawbacks to this style:
@ -1136,7 +1134,7 @@ a matter of taste but there are severe drawbacks to this style:
repetitions. It is easy to forget a single change.
\item Readability: repetitive code is terrible to read and to
understand. (I) one tends to skip repetitions (its the same,
anyways) and misses the essential change. (II), the duplication of
anyway) and misses the essential change. (II), the duplication of
code leads to long and hard-to-parse programs.
\end{enumerate}
All imperative programming languages offer a solution: the loop. It is
@ -1146,7 +1144,7 @@ used whenever the same commands have to be repeated.
\subsubsection{The \code{for} --- loop}
The most common type of loop is the \codeterm{for-loop}. It
consists of a \codeterm[Loop!head]{head} and the
\codeterm[Loop!body]{body}. The head defines how often the code of the
\codeterm[Loop!body]{body}. The head defines how often the code in the
body is executed. In \matlab{} the head begins with the keyword
\code{for} which is followed by the \codeterm{running variable}. In
\matlab{} a for-loop always operates on vectors. With each
@ -1191,7 +1189,7 @@ while x == true % head with a Boolean expression
% execute this code if the expression yields true
end
\end{lstlisting}
\begin{exercise}{factorialWhileLoop.m}{}
Implement the factorial of a number \varcode{n} using a \code{while}
-- loop.
@ -1201,7 +1199,7 @@ end
\begin{exercise}{neverendingWhile.m}{}
Implement a \code{while}--loop that is never-ending. Hint: the body
is executed as long as the Boolean expression in the head is
true. You can escape the loop by pressing \keycode{Ctrl+C}.
\code{true}. You can escape the loop by pressing \keycode{Ctrl+C}.
\end{exercise}

8
programmingstyle/code/calculateSines.m
View File

@ -1,5 +1,7 @@
function sines = calculateSines(x, amplitudes, frequencies)
% Function calculates sinewaves with all combinations of
% sines = calculateSines(x, amplitudes, frequencies)
%
% Function calculates sinewaves with all combinations of
% given amplitudes and frequencies.
% Arguments: x, a vector of radiants for which the sine should be
% computed.
@ -8,11 +10,11 @@ function sines = calculateSines(x, amplitudes, frequencies)
%
% Returns: a 3-D Matrix of sinewaves, 2nd dimension represents
% the amplitudes, 3rd the frequencies.
sines = zeros(length(x), length(amplitudes), length(frequencies));
for i = 1:length(amplitudes)
sines(:,i,:) = sinesWithFrequencies(x, amplitudes(i), frequencies);
sines(:,i,:) = sinesWithFrequencies(x, amplitudes(i), frequencies);
end
end

157
programmingstyle/lecture/programmingstyle.tex
View File

@ -6,20 +6,19 @@
%\selectlanguage{ngerman}
Cultivating a good code style not a matter of good taste but is
a key ingredient for understandability, maintainability and in the end
a key ingredient for understandability, maintainability and, in the end,
facilitates reproducibility of scientific results. Programs should be
written and structured in a way that supports outsiders as well the
author himself --- a few weeks or months after it was written --- to
understand the programs' rationale. Clean code pays off for the
original author as well as others that are supposed to use the code.
Clean code addresses several issues:
\begin{enumerate}
\item The programs' structure.
\item Naming of scripts and functions.
\item Naming of variables and constants.
\item Application of indentation empty lines to define blocks.
\item Application of indentation and empty lines to define blocks.
\item Use of comments and inline documentation.
\item Delegation of repeated code to functions and dedicated
subroutines.
@ -29,13 +28,12 @@ Clean code addresses several issues:
While introducing scripts and functions we suggested a typical program
layout (box\,\ref{whenscriptsbox}). The idea is to create a single
entry point by having one script that controls the rest of program by
managing data and results and calling functions that work on the data
and produce the results. Applying this structure makes it easy to
understand the flow of the program but two questions remain: (i) How
to organize the files on the file system and (ii) how to name them
that the controlling script is easily identified among the other
\codeterm{m-files}.
entry point by having one script that controls the rest of the program
by calling functions that work on the data and managing the
results. Applying this structure makes it easy to understand the flow
of the program but two questions remain: (i) How to organize the files
on the file system and (ii) how to name them that the controlling
script is easily identified among the other \codeterm{m-files}.
Upon installation ``MATLAB'' creates a folder called \emph{MATLAB} in
the user space (Windows: My files, Linux: Documents, MacOS:
@ -110,24 +108,11 @@ Box~\ref{matlabpathbox}).
\end{lstlisting}
\end{ibox}
\section{Naming scripts and functions}
\matlab{} will search the search path (Box \ref{matlabpathbox})
exclusively by name. It is case-sensitive this implies that the files
\file{test\_function.m} and \file{Test\_function.m} are two different
things. It is self-evident that choosing such names is nonsensical
because the name contains no cue about the difference between the two
and it further tells close to nothing about the purpose. Finding good
names is not trivial sometimes it is harder than the programming
itself. Expressive names, however, pay off! Expressive means that the
name provides information about the purpose.
\begin{important}[Naming scripts and functions]
Function and script names should be expressive in the sense that the
name provides information about the function's purpose
(\file{estimate\_firingrate.m} tells much more than
\file{exercise1.m}). Choosing a good name replaces large parts of
the documentation.
\end{important}
\section{Naming things}
The dictum of good code style is: ``Program code must be readable.''
Expressive names are extraordinarily important in this respect. Even
if it is tricky to find expressive names that are not overly long,
naming should be taken seriously.
\matlab{} has a few rules about names: Names must not start with a
number, they must not contain blanks or other special characters like
@ -144,26 +129,41 @@ patterns:
There are other common patterns such as the \emph{camelCase} in which
the first character of compound words is capitalized. Other
conventions use the underscore to separate the individual words
\emph{snake\_case}. A function that counts the number of action
conventions use the underscore to separate the individual words (
\emph{snake\_case}). A function that counts the number of action
potentials could be named \file{spikeCount.m} or
\file{spike\_count.m}.
The same naming rules apply for scripts and functions as well as
variables and constants.
\section{Naming variables and constants}
\subsection{Naming scripts and functions}
\matlab{} will search the search path (Box \ref{matlabpathbox})
exclusively by name. This search is case-sensitive which implies that
the files \file{test\_function.m} and \file{Test\_function.m} are two
different things. It is self-evident that choosing such names is
nonsensical because the tiny difference in the name contains no cue
about the difference between the two versions and the function names
themselves tell close to nothing about the purpose. Finding good names
is not trivial. Sometimes it is harder than the programming
itself. Choosing \emph{expressive names} that provide information about a
function's purpose, however, pays off!
\begin{important}[Naming scripts and functions]
Names of functions and scripts should be expressive in the sense
that the name provides information about the function's purpose.
(\file{estimate\_firingrate.m} tells much more than
\file{exercise1.m}). Choosing a good name replaces large parts of
the documentation.
\end{important}
\matlab{} applies the same rules for naming variables and constants as
for the naming of scripts and functions. The dictum of good
code style is: ``Program code must be readable.'' Expressive
names are extraordinarily important in this respect. Even if it is
tricky to find expressive names that are not overly long, naming
should be taken seriously.
\subsection{Naming variables and constants}
While the names of scripts and functions describe the purpose, names
of variables describe the stored content. A variable storing the
average number of actions potentials could be called
average number of actions potentials could be called\\
\varcode{average\_spike\_count}. If this variable is meant to store
multiple spike counts the plural form would be appropriate
multiple spike counts the plural form would be appropriate\\
(\varcode{average\_spike\_counts}).
The control variables used in the head of a \code{for} loop are often
@ -190,7 +190,7 @@ to comprehend. Even though the \matlab{} language (as many others)
does not enforce indentation, indentation is very powerful for
defining coherent blocks. The \matlab{} editor supports this by an
auto-indentation mechanism. A selected section of the code and be
automatically indented by pressing the \keycode{Ctrl-I} combination.
automatically indented by pressing \keycode{Ctrl-I}.
Interspersing empty lines is very helpful to separate regions in the
code that belong together. Too many empty lines, however lead to
@ -198,8 +198,11 @@ hard-to-read code because it might require more space than a granted
by the screen and thus takes overview.
The following two listings show basically the same implementation of a
random walk once in a rather chaotic version (listing
\ref{chaoticcode}) then in cleaner way (listing \ref{cleancode})
random walk\footnote{A random walk is a simple simulation of Brownian
motion. In each simulation step an agent takes a step into a
randomly chosen direction.} once in a rather chaotic version
(listing \ref{chaoticcode}) then in cleaner way (listing
\ref{cleancode})
\begin{lstlisting}[label=chaoticcode, caption={Chaotic implementation of the random-walk.}]
num_runs = 10; max_steps = 1000;
@ -245,25 +248,24 @@ end
\section{Using comments}
It is common to provide extra information about the meaning of program
code by adding comments to it. In \matlab{} comments are indicated by
the percent character \code{\%}. Anything that is written in the
respective line following the percent is ignored and considered a
comment. When used sparsely comments can immensely important for
understanding. Comments are short sentences that describe the meaning
of the (following) lines in the program code. During the initial
implementation of a function they can be used to guide the development
but have the tendency to blow up the code and decrease readability. By
choosing expressive variable and function names, most lines should be
self-explanatory.
For example stating the obvious does not really help:\\
\varcode{ x = x + 2; \% add two to x}\\
code by adding comments. In \matlab{} comments are indicated by the
percent character \code{\%}. Anything that follows the percent
character in a line is ignored and considered a comment. When used
sparsely comments can be immensely helpful. Comments
are short sentences that describe the meaning of the (following) lines
in the program code. During the initial implementation of a function
they can be used to guide the development but have the tendency to
blow up the code and decrease readability. By choosing expressive
variable and function names, most lines should be self-explanatory.
For example stating the obvious does not really help and should be
avoided:\\ \varcode{ x = x + 2; \% add two to x}\\
\begin{important}[Using comments]
\begin{itemize}
\item Comments describe the rationale of the respective code block.
\item Comments are good and helpful --- they have to be true, however!
\item A wrong comment is worse than a non-existent comment!
\item Comments are good and helpful --- they must be true, however!
\item A wrong comment is worse than a non-existent one!
\item Comments must be maintained just as the code. Otherwise they
may become wrong and worse than meaningless!
\end{itemize}
@ -303,7 +305,7 @@ well documented function.
Comments and empty lines are used to organize code into logical blocks
and to briefly explain what they do. Whenever one feels tempted to do
this, one could also consider to delegate the respective task to a
function. In most cases this is preferable.
function. In most cases this is preferable.
Not delegating the tasks leads to very long \codeterm{m-files} which
can be confusing. Sometimes such a code is called ``spaghetti
@ -325,9 +327,9 @@ Generally, functions live in their own \codeterm{m-files} that have
the same name as the function itself. Delegating tasks to functions
thus leads to a large set of \codeterm{m-files} which increases
complexity and may lead to confusion. If the delegated functionality
is used in multiple instances, it is advisable to do so. On the other
hand, when the delegated functionality is only used within the context
of another function \matlab{} allows to define
is used in multiple instances, it is still advisable to do so. On the
other hand, when the delegated functionality is only used within the
context of another function \matlab{} allows to define
\codeterm[function!local]{local functions} and
\codeterm[function!nested]{nested functions} within the same
file. Listing \ref{localfunctions} shows an example of a local
@ -336,17 +338,19 @@ function definition.
\pagebreak[3] \lstinputlisting[label=localfunctions, caption={Example
for local functions.}]{calculateSines.m}
Local function live in the same \codeterm{m-file} as the main function
and are only available in this context. Each local function has its
own \codeterm{scope}, that is, the local function can not access (read
or write) variables of the calling function.
\emph{Local function} live in the same \codeterm{m-file} as the main
function and are only available in this context. Each local function
has its own \codeterm{scope}, that is, the local function can not
access (read or write) variables of the calling function. Interaction
with the local function requires to pass all required arguments and to
take care of the return values of the function.
This is different in so called \codeterm[function!nested]{nested
functions}. These are defined within the body of the parent function
(between the keywords \code{function} and \code{end}) and have full
access to all variables defined in the parent function. Working (in
particular changing) the parent's variables is handy on the one side,
but is also risky. One should take care when defining nested functions.
\emph{Nested functions} are different in this respect. They are
defined within the body of the parent function (between the keywords
\code{function} and \code{end}) and have full access to all variables
defined in the parent function. Working (in particular changing) the
parent's variables is handy on the one side, but is also risky. One
should take care when defining nested functions.
\section{Specifics when using scripts}
@ -355,7 +359,7 @@ A similar problem as with nested function arises when using scripts
become available in the global \codeterm{Workspace}. There is the risk
of name conflicts, that is, a called sub-script redefines or uses the
same variable name and may \emph{silently} change its content. The
user will not be notified by this change and the calling script may
user will not be notified about this change and the calling script may
expect a completely different content. Bugs that are based on such
mistakes are hard to find since the program itself looks perfectly
fine.
@ -379,12 +383,15 @@ provides important information to track and fix the bug.
\item Scripts should work independently of existing variables in the
global workspace.
\item It is advisable to start a script with deleting variables
(\code{clear}) from the workspace and most of the times it is also
good to close all open figures (\code{close all}).
\item Often it is advisable to start a script with deleting
variables (\code{clear}) from the workspace and most of the times
it is also good to close all open figures (\code{close all}). Be
careful if a the respective script has been called by another one.
\item Clean up the workspace at the end of a script. Delete
(\code{clear}) all variables that are no longer needed.
\item Consider to write functions instead of scripts.
\end{itemize}
\end{important}

8
regression/lecture/regression.tex
View File

@ -86,7 +86,7 @@ large deviations.
$f_{cost}(\{(x_i, y_i)\}|\{y^{est}_i\})$ is a so called
\enterm{objective function} or \enterm{cost function}. We aim to adapt
the model parameters to minimize the error (mean square error) and
thus the \emph{objective function}. In Chapter~\ref{maximumlikelihood}
thus the \emph{objective function}. In Chapter~\ref{maximumlikelihoodchapter}
we will show that the minimization of the mean square error is
equivalent to maximizing the likelihood that the observations
originate from the model (assuming a normal distribution of the data
@ -270,7 +270,7 @@ The gradient is given by partial derivatives
(Box~\ref{partialderivativebox}) with respect to the parameters $m$
and $b$ of the linear equation. There is no need to calculate it
analytically but it can be estimated from the partial derivatives
using the difference quotient (Box~\ref{differentialquotient}) for
using the difference quotient (Box~\ref{differentialquotientbox}) for
small steps $\Delta m$ und $\Delta b$. For example the partial
derivative with respect to $m$:
@ -383,7 +383,7 @@ introduce how this can be done without using the gradient descent
Problems that involve nonlinear computations on parameters, e.g. the
rate $\lambda$ in the exponential function $f(x;\lambda) =
\exp(\lambda x)$, do not have an analytical solution. To find minima
e^{\lambda x}$, do not have an analytical solution. To find minima
in such functions numerical methods such as the gradient descent have
to be applied.
@ -396,7 +396,7 @@ objective functions while more specialized functions are specifically
designed for optimizations in the least square error sense
\matlabfun{lsqcurvefit()}.
\newpage
%\newpage
\begin{important}[Beware of secondary minima!]
Finding the absolute minimum is not always as easy as in the case of
the linear equation. Often, the error surface has secondary or local

2
scientificcomputing-script.tex
View File

@ -67,7 +67,7 @@
\lstset{inputpath=regression/code}
\include{regression/lecture/regression}
\setboolean{showexercisesolutions}{true}
\setboolean{showexercisesolutions}{false}
\graphicspath{{likelihood/lecture/}{likelihood/lecture/figures/}}
\lstset{inputpath=likelihood/code}
\include{likelihood/lecture/likelihood}

4
statistics/lecture/displayunivariatedata.py
View File

@ -53,7 +53,7 @@ ax.annotate('mean plus\nstd. dev.',
arrowprops=dict(arrowstyle="->", relpos=(0.5,0.0),
connectionstyle="angle3,angleA=-60,angleB=80") )
ax = fig.add_axes([xpos, ypos, width, height], axis_bgcolor='none')
ax = fig.add_axes([xpos, ypos, width, height])
ax.spines['right'].set_visible(False)
ax.spines['top'].set_visible(False)
ax.spines['left'].set_visible(False)
@ -92,7 +92,7 @@ ax.annotate('median',
arrowprops=dict(arrowstyle="->", relpos=(0.8,0.0),
connectionstyle="angle3,angleA=-60,angleB=20") )
ax = fig.add_axes([xpos+width+0.03, ypos, 0.98-(xpos+width+0.03), height], axis_bgcolor='none')
ax = fig.add_axes([xpos+width+0.03, ypos, 0.98-(xpos+width+0.03), height])
ax.spines['right'].set_visible(False)
ax.spines['top'].set_visible(False)
ax.xaxis.set_ticks_position('bottom')