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: # bootstrap the mean:
mus = [] mus = []
for i in xrange(nresamples) : for i in range(nresamples) :
mus.append(np.mean(x[rng.randint(0, nsamples, nsamples)])) mus.append(np.mean(x[rng.randint(0, nsamples, nsamples)]))
hmus, _ = np.histogram(mus, bins, density=True) hmus, _ = np.histogram(mus, bins, density=True)
# many SRS: # many SRS:
musrs = [] musrs = []
for i in xrange(nresamples) : for i in range(nresamples) :
musrs.append(np.mean(rng.randn(nsamples))) musrs.append(np.mean(rng.randn(nsamples)))
hmusrs, _ = np.histogram(musrs, bins, density=True) 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: # permutation:
nperm = 1000 nperm = 1000
rs = [] rs = []
for i in xrange(nperm) : for i in range(nperm) :
xr=rng.permutation(x) xr=rng.permutation(x)
yr=rng.permutation(y) yr=rng.permutation(y)
rs.append( np.corrcoef(xr, yr)[0, 1] ) rs.append( np.corrcoef(xr, yr)[0, 1] )

3
chapter.mk
View File

@ -8,7 +8,8 @@ PYPDFFILES=$(PYFILES:.py=.pdf)
pythonplots : $(PYPDFFILES) pythonplots : $(PYPDFFILES)
$(PYPDFFILES) : %.pdf: %.py $(PYPDFFILES) : %.pdf: %.py
python $< echo $$(which python)
python3 $<
cleanpythonplots : cleanpythonplots :
rm -f $(PYPDFFILES) 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') ax.text(mu-0.1, 0.04, '?', zorder=1, ha='right')
else : else :
ax.text(mu+0.1, 0.04, '?', zorder=1) 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.plot(x, g[:,k], zorder=5)
ax.scatter(xd, 0.05*rng.rand(len(xd))+0.2, s=30, zorder=10) 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} severe.}{\url{www.xkcd.com}}\label{xkcdplotting}
\end{figure} \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} \section{The \matlab{} plotting system}
Plotting data in \matlab{} is rather straight forward for simple line 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. number of datasets.
\begin{important}[Why manual editing should be avoided.] \begin{important}[Why manual editing should be avoided.]
On first glance the manual editing of a figure using common tools On first glance the manual editing of a figure using tools such as
like Corel draw, Illustrator, etc.\,appears some much more Corel draw, Illustrator, etc.\,appears much more convenient and less
convenient and less complex. This, however, is not entirely complex than coding everything into the analysis scripts. This,
true. What if the figure has to be re-drawn or updated? Then the however, is not entirely true. What if the figure has to be re-drawn
editing work starts all over again. Rather, there is a great risk or updated? Then the editing work starts all over again. Rather,
associated with this approach. Axes are shifted, fonts have not been there is a great risk associated with the manual editing
embedded into the final document, annotations have been copy pasted approach. Axes may be shifted, fonts have not been embedded into the
between figures and are not valid. All of these mistakes can be final document, annotations have been copy pasted between figures
found in publications and then require an erratum, which is not 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 desirable. Even if it appears more cumbersome in the beginning one
should always try to create publication-ready figures directly from should always try to create publication-ready figures directly from
the data analysis tool using scripts or functions to properly the data analysis tool using scripts or functions to properly layout
layout the plot. the plot.
\end{important} \end{important}
\subsection{Simple plotting} \subsection{Simple plotting}
@ -259,6 +143,12 @@ additional options consult the help.
\end{tabular} \end{tabular}
\end{table} \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} \subsection{Changing properties of a line plot}
The properties of line plots can be changed by passing more arguments The properties of line plots can be changed by passing more arguments
@ -419,14 +309,14 @@ output format (box\,\ref{graphicsformatbox}).
\end{tabular} \end{tabular}
\end{minipage} \end{minipage}
It is often meaningful to store of data plots generated by \matlab{} It is advisable to store of data plots generated by \matlab{}
using a vector graphics format. When in doubt they can usually be using a vector graphics format. In doubt they can usually be
easily converted to a bitmap format. The way from a bitmap to a 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 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 raster-plot showing thousands of action potentials) may, on the
other hand, lead to very large files that can be hard to 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. efficient.
\end{ibox} \end{ibox}
@ -584,24 +474,24 @@ its properties. See the \matlab{} help for more information.
\begin{figure}[ht] \begin{figure}[ht]
\includegraphics[width=0.9\linewidth]{errorbars} \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 symmetrical error around the mean (e.g.\ using the standard
deviation). \textbf{B} Errorbars of an asymmetrical distribution deviation). \textbf{B} Errorbars of an asymmetrical distribution
of the data (note: the average value is now the median and the of the data (note: the average value is now the median and the
errors are the lower and upper quartiles). \textbf{C} A shaded errors are the lower and upper quartiles). \textbf{C} A shaded
area is used to illustrate the spread of the data. See 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} \end{figure}
\lstinputlisting[caption={Illustrating estimation errors. Script that \lstinputlisting[caption={Illustrating estimation errors using error bars. Script that
creates \figref{errorbarplot}.}, creates \figref{errorbarplot}. A, B},
label=errorbarlisting, firstline=13, lastline=29, label=errorbarlisting, firstline=13, lastline=29,
basicstyle=\ttfamily\scriptsize]{errorbarplot.m} basicstyle=\ttfamily\scriptsize]{errorbarplot.m}
\subsubsection{Fill} \subsubsection{Fill}
For a few years now it has become fancy to illustrate the error not 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 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 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 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 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 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 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 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} version of the x-values to the original x-values using \code{cat} and
and inversion is done using the \code{fliplr} command (line 3 in \code{fliplr} for concatenation and inversion, respectively (line 3 in
listing \ref{errorbarlisting2}; Depending on the layout of your data listing \ref{errorbarlisting2}; Depending on the layout of your data
you may need concatenate along a different dimension of the data and you may need concatenate along a different dimension of the data and
use \code{flipud} instead). The y-coordinates of the polygon vertices 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 area. The provided alpha value is a number between 0 and 1 with zero
leading to invisibility and a value of one to complete leading to invisibility and a value of one to complete
opaqueness. Finally, we use the normal plot command to draw a line 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 \lstinputlisting[caption={Illustrating estimation errors using a shaded area. Script that
creates \figref{errorbarplot}.}, label=errorbarlisting2, creates \figref{errorbarplot} C.}, label=errorbarlisting2,
firstline=30, firstline=30,
basicstyle=\ttfamily\scriptsize]{errorbarplot.m} basicstyle=\ttfamily\scriptsize]{errorbarplot.m}
\subsection{Annotations, text} \subsection{Annotations, text}
Sometimes want to highlight certain parts of a plot or simply add an The \code[text()]{text()} or \code[annotation()]{annotation()} are
annotation that does not fit or belong to the legend. In these cases used for highlighting certain parts of a plot or simply adding an
we can use the \code[text()]{text()} or annotation that does not fit or does not belong into the legend.
\code[annotation()]{annotation()} function to add this information to While \varcode{text} simply prints out the given text string at the
the plot. While \varcode{text} simply prints out the given text string defined position (for example line in
at the defined position (for example line in
listing\,\ref{regularsubplotlisting}) the \varcode{annotation} listing\,\ref{regularsubplotlisting}) the \varcode{annotation}
function allows to add some more advanced highlights like arrows, function allows to add some more advanced highlights like arrows,
lines, ellipses, or rectangles. Figure\,\ref{annotationsplot} shows lines, ellipses, or rectangles. Figure\,\ref{annotationsplot} shows
some examples, the respective code can be found in some examples, the respective code can be found in
listing\,\ref{annotationsplotlisting}. For more options consult the listing\,\ref{annotationsplotlisting}. For more options consult the
documentation. \matlab{} help.
\begin{figure}[ht] \begin{figure}[ht]
\includegraphics[width=0.5\linewidth]{annotations} \includegraphics[width=0.5\linewidth]{annotations}
@ -714,20 +603,133 @@ Lissajous figure. The basic steps are:
movie.}, label=animationlisting, firstline=3, lastline=33, movie.}, label=animationlisting, firstline=3, lastline=33,
basicstyle=\ttfamily\scriptsize]{movie_example.m} 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} \section{Summary}
A good plot of scientific data displays the data completely and A good plot of scientific data displays the data completely and
seriously without too many distractions. Misleading or suggestive seriously without too many distractions. Misleading or suggestive
plots as may result from perspective presentations, inappropriate 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 \noindent When combining several line plots within the same figure one should
consider adapting color \textbf{and} line style (solid, dashed, consider adapting color \textbf{and} line style (solid, dashed,
dotted. etc.) to make the distinguishable even in black-and-white 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. cannot be distinguished by people with red-green blindness.
\vspace{2ex} \vspace{2ex}
@ -737,6 +739,8 @@ Key ingredients for a good data plot:
\item Complete labeling. \item Complete labeling.
\item Plotted lines and curves must be distinguishable. \item Plotted lines and curves must be distinguishable.
\item No suggestive or misleading presentation. \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. \item Error bars wherever they are appropriate.
\end{itemize} \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): def get_binned_rate(times, bin_width=0.05, max_t=30., dt=1e-4):
time = np.arange(0., max_t, dt) time = np.arange(0., max_t, dt)
bins = np.arange(0., max_t, bin_width) 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) hist, _ = sp.histogram(times, bins)
rate = np.zeros(time.shape) 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 = [] times = []
v = vreset v = vreset
noise = np.sqrt(2.0*D)*np.random.randn(len(input))/np.sqrt(dt) 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 v += (input[k]+noise[k])*dt/tau
if v >= vthresh : if v >= vthresh :
v = vreset v = vreset
@ -41,7 +41,7 @@ def pifspikes(input, trials, dt, D=0.1) :
def isis( spikes ) : def isis( spikes ) :
isi = [] isi = []
for k in xrange(len(spikes)) : for k in range(len(spikes)) :
isi.extend(np.diff(spikes[k])) isi.extend(np.diff(spikes[k]))
return isi return isi
@ -76,7 +76,7 @@ rng = np.random.RandomState(54637281)
time = np.arange(0.0, duration, dt) time = np.arange(0.0, duration, dt)
x = np.zeros(time.shape)+rate x = np.zeros(time.shape)+rate
n = rng.randn(len(time))*drate*tau/np.sqrt(dt)+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[k] = x[k-1] + (n[k]-x[k-1])*dt/tau
x[x<0.0] = 0.0 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 \enterm[Actionspotentials]{Actionspotentials} (\enterm{spikes}) are
Information in Nervensystemen. Dabei ist in erster Linie nur der the carriers of information in the nervous system. Thereby it is
Zeitpunkt des Auftretens eines Aktionspotentials von Bedeutung. Die mainly the time at which the spikes are generated that is of
genaue Form des Aktionspotentials spielt keine oder nur eine importance. The waveform of the action potential is largely
untergeordnete Rolle. stereotyped and does not carry information.
Nach etwas Vorverarbeitung haben elektrophysiologische Messungen The result of the processing of electrophysiological recordings are
deshalb Listen von Spikezeitpunkten als Ergebniss --- sogenannte series of spike times, which are then termed \enterm{spiketrains}. If
\enterm{spiketrains}. Diese Messungen k\"onnen wiederholt werden und measurements are repeated we yield several \enterm{trials} of
es ergeben sich mehrere \enterm{trials} von Spiketrains spiketrains (\figref{rasterexamplesfig}).
(\figref{rasterexamplesfig}).
Spiketrains sind Zeitpunkte von Ereignissen --- den Aktionspotentialen Spiketrains are times of events, the action potentials. The analysis
--- und deren Analyse f\"allt daher in das Gebiet der Statistik von of these leads into the realm of the so called \enterm[point
sogenannten \determ[Punktprozess]{Punktprozessen}. process]{point processes}.
\begin{figure}[ht] \begin{figure}[ht]
\includegraphics[width=1\textwidth]{rasterexamples} \includegraphics[width=1\textwidth]{rasterexamples}
\titlecaption{\label{rasterexamplesfig}Raster-Plot.}{Raster-Plot von \titlecaption{\label{rasterexamplesfig}Raster-plot.}{Raster-plot of
jeweils 10 Realisierungen eines station\"arenen Punktprozesses ten realizations of a stationary point process (homogeneous point
(homogener Poisson Prozess mit Rate $\lambda=20$\;Hz, links) und process with a rate $\lambda=20$\;Hz, left) and an inhomogeneous
eines nicht-station\"aren Punktprozesses (perfect point process (perfect integrate-and-fire neuron dirven by
integrate-and-fire Neuron getrieben mit Ohrnstein-Uhlenbeck Ohrnstein-Uhlenbeck noise with a time-constant $\tau=100$\,ms,
Rauschen mit Zeitkonstante $\tau=100$\,ms, rechts). Jeder right). Each vertical dash illustrates the time at which the
vertikale Strich markiert den Zeitpunkt eines Ereignisses. action potential was observed. Each line represents the event of
Jede Zeile zeigt die Ereignisse eines trials.} each trial.}
\end{figure} \end{figure}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Punktprozesse} \section{Point processes}
Ein zeitlicher Punktprozess (\enterm{point process}) ist ein A temporal \enterm{point process} is a stochastic process that
stochastischer Prozess, der eine Abfolge von Ereignissen zu den Zeiten generates a sequence of events at times $\{t_i\}$, $t_i \in
$\{t_i\}$, $t_i \in \reZ$, generiert. \reZ$.
\begin{ibox}{Beispiele von Punktprozessen} \begin{ibox}{Examples of point processes}
Jeder Punktprozess wird durch einen sich in der Zeit kontinuierlich Every point process is generated by a temporally continuously
entwickelnden Prozess generiert. Wann immer dieser Prozess eine developing process. An event is generated whenever this process
Schwelle \"uberschreitet wird ein Ereigniss des Punktprozesses reaches a certain threshold. For example:
erzeugt. Zum Beispiel: \begin{itemize}
\begin{itemize} \item Action potentials/heart beat: created by the dynamics of the
\item Aktionspotentiale/Herzschlag: wird durch die Dynamik des neuron/sinoatrial node
Membranpotentials eines Neurons/Herzzelle erzeugt. \item Earthquake: defined by the dynamics of the pressure between
\item Erdbeben: wird durch die Dynamik des Druckes zwischen tectonical plates.
tektonischen Platten auf beiden Seiten einer geologischen Verwerfung \item Evoked communication calls in crickets/frogs/birds: shaped by
erzeugt. the dynamics of nervous system and the muscle appartus.
\item Zeitpunkt eines Grillen/Frosch/Vogelgesangs: wird durch die \end{itemize}
Dynamik des Nervensystems und des Muskelapparates erzeugt.
\end{itemize}
\end{ibox} \end{ibox}
\begin{figure}[t] \begin{figure}[t]
\texpicture{pointprocessscetch} \texpicture{pointprocessscetch}
\titlecaption{\label{pointprocessscetchfig} Statistik von \titlecaption{\label{pointprocessscetchfig} Statistics of point
Punktprozessen.}{Ein Punktprozess ist eine Abfolge von processes.}{A point process is a sequence of instances in time
Zeitpunkten $t_i$ die auch durch die Intervalle $T_i=t_{i+1}-t_i$ $t_i$ that can be characterized through the inter-event-intervals
oder die Anzahl der Ereignisse $n_i$ beschrieben werden kann. } $T_i=t_{i+1}-t_i$ and the number of events $n_i$. }
\end{figure} \end{figure}
F\"ur die Neurowissenschaften ist die Statistik der Punktprozesse In the neurosciences, the statistics of point processes is of
besonders wichtig, da die Zeitpunkte der Aktionspotentiale als importance since the timing of the neuronal events (the action
zeitlicher Punktprozess betrachtet werden k\"onnen und entscheidend potentials) is crucial for information transmission and can be treated
f\"ur die Informations\"ubertragung sind. as such a process.
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 Point processes can be described using the intervals between
\figref{rasterexamplesfig} als Rasterplot dargestellt, bei dem die successive events $T_i=t_{i+1}-t_i$ and the number of observed events
Zeitpunkte der Ereignisse durch senkrechte Striche markiert werden. 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 The intervals $T_i=t_{i+1}-t_i$ between successive events are real
Ereignissen sind reelle, positive Zahlen. Bei Aktionspotentialen positive numbers. In the context of action potentials they are
heisen die Intervalle auch \determ{Interspikeintervalle} referred to as \enterm{interspike intervals}. The statistics of these
(\enterm{interspike intervals}). Deren Statistik kann mit den are described using the common measures.
\"ublichen Gr\"o{\ss}en beschrieben werden.
\begin{figure}[t] \begin{figure}[t]
\includegraphics[width=0.96\textwidth]{isihexamples}\vspace{-2ex} \includegraphics[width=0.96\textwidth]{isihexamples}\vspace{-2ex}
\titlecaption{\label{isihexamplesfig}Interspikeintervall Histogramme}{der in \titlecaption{\label{isihexamplesfig}Interspike interval
\figref{rasterexamplesfig} gezeigten Spikes.} histogram}{of the spikes depicted in \figref{rasterexamplesfig}.}
\end{figure} \end{figure}
\begin{exercise}{isis.m}{} \begin{exercise}{isis.m}{}
Schreibe eine Funktion \code{isis()}, die aus mehreren trials von Spiketrains die Implement a function \code{isis()} that calculates the interspike
Interspikeintervalle bestimmt und diese in einem Vektor intervals from several spike trains. The function should return a
zur\"uckgibt. Jeder trial der Spiketrains ist ein Vektor mit den single vector of intervals. The action potentials recorded in the
Spikezeiten gegeben in Sekunden als Element in einem \codeterm{cell-array}. individual trials are stored as vectors of spike times within a
\codeterm{cell-array}. Spike times are given in seconds.
\end{exercise} \end{exercise}
\subsection{Intervallstatistik erster Ordnung} \subsection{First order interval statistics}
\begin{itemize} \begin{itemize}
\item Wahrscheinlichkeitsdichte $p(T)$ der Intervalle $T$ \item Probability density $p(T)$ of the intervals $T$
(\figref{isihexamplesfig}). Normiert auf $\int_0^{\infty} p(T) \; dT (\figref{isihexamplesfig}). Normalized to $\int_0^{\infty} p(T) \; dT
= 1$. = 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$. \frac{1}{n}\sum\limits_{i=1}^n T_i$.
\item Standardabweichung der Intervalle: $\sigma_{ISI} = \sqrt{\langle (T - \langle T \item Standard deviation of the interspike intervals: $\sigma_{ISI} = \sqrt{\langle (T - \langle T
\rangle)^2 \rangle}$\vspace{1ex} \rangle)^2 \rangle}$\vspace{1ex}
\item \determ{Variationskoeffizient} (\enterm{coefficient of variation}): $CV_{ISI} = \item \enterm{Coefficient of variation}: $CV_{ISI} =
\frac{\sigma_{ISI}}{\mu_{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}$. \frac{\sigma_{ISI}^2}{2\mu_{ISI}^3}$.
\end{itemize} \end{itemize}
\begin{exercise}{isihist.m}{} \begin{exercise}{isihist.m}{}
Schreibe eine Funktion \code{isiHist()}, die einen Vektor mit Interspikeintervallen Implement a function \code{isiHist()} that calculates the normalized
entgegennimmt und daraus ein normiertes Histogramm der Interspikeintervalle interspike interval histogram. The function should take two input
berechnet. 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} \end{exercise}
\begin{exercise}{plotisihist.m}{} \begin{exercise}{plotisihist.m}{}
Schreibe eine Funktion, die die Histogrammdaten der Funktion Implement a function that takes the return values of
\code{isiHist()} entgegennimmt, um das Histogramm zu plotten. Im \code{isiHist()} as input arguments and then plots the data. The
Plot sollen die Interspikeintervalle in Millisekunden aufgetragen plot should show the histogram with the x-axis scaled to
werden. Das Histogramm soll zus\"atzlich mit Mittelwert, milliseconds and should be annotated with the average ISI, the
Standardabweichung und Variationskoeffizient der standard deviation and the coefficient of variation.
Interspikeintervalle annotiert werden.
\end{exercise} \end{exercise}
\subsection{Korrelationen der Intervalle} \subsection{Interval correlations}
In \enterm{return maps} werden die um das \enterm{lag} $k$ verz\"ogerten So called \enterm{return maps} are used to illustrate
Intervalle $T_{i+k}$ gegen die Intervalle $T_i$ geplottet. Dies macht interdependencies between successive interspike intervals. The return
m\"ogliche Abh\"angigkeiten von aufeinanderfolgenden Intervallen map plots the delayed interval $T_{i+k}$ against the interval
sichtbar. $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] \begin{figure}[t]
\includegraphics[width=1\textwidth]{returnmapexamples} \includegraphics[width=1\textwidth]{returnmapexamples}
\includegraphics[width=1\textwidth]{serialcorrexamples} \includegraphics[width=1\textwidth]{serialcorrexamples}
\titlecaption{\label{returnmapfig}Interspikeintervall return maps und \titlecaption{\label{returnmapfig}Interspike interval analyses of a
serielle Korrelationen}{zwischen aufeinander folgenden Intervallen stationary and a non-stationary pointprocess.}{Upper plots show the
im Abstand des Lags $k$.} return maps and the lower panels depict the serial correlation of
successive intervals separated by the lag $k$.}
\end{figure} \end{figure}
Solche Ab\"angigkeiten werden durch die \determ{serielle Such dependencies can be further quantified using the \enterm{serial
Korrelationen} (\enterm{serial correlations}) der Intervalle correlations} \figref{returnmapfig}. The serial correlation is the
quantifiziert. Das ist der \determ{Korrelationskoeffizient} zwischen correlation coefficient of the intervals $T_i$ and the intervals
aufeinander folgenden Intervallen getrennt durch lag $k$: 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)} \[ \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) \] = {\rm corr}(T_{i+k}, T_i) \] The resulting correlation coefficient
\"Ublicherweise wird die Korrelation $\rho_k$ gegen den Lag $k$ $\rho_k$ is usually plotted against the lag $k$
aufgetragen (\figref{returnmapfig}). $\rho_0=1$ (Korrelation jedes \figref{returnmapfig}. $\rho_0=1$ is the correlation of each interval
Intervalls mit sich selber). with itself and is always 1.
\begin{exercise}{isiserialcorr.m}{} \begin{exercise}{isiserialcorr.m}{}
Schreibe eine Funktion \code{isiserialcorr()}, die einen Vektor mit Interspikeintervallen Implement a function \code{isiserialcorr()} that takes a vector of
entgegennimmt und daraus die seriellen Korrelationen berechnet und plottet. interspike intervals as input argument and calculates the serial
\pagebreak[4] correlation. The function should further plot the serial
correlation. \pagebreak[4]
\end{exercise} \end{exercise}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Z\"ahlstatistik} \section{Count statistics}
% \begin{figure}[t] % \begin{figure}[t]
% \includegraphics[width=0.48\textwidth]{poissoncounthist100hz10ms}\hfill % \includegraphics[width=0.48\textwidth]{poissoncounthist100hz10ms}\hfill
% \includegraphics[width=0.48\textwidth]{poissoncounthist100hz100ms} % \includegraphics[width=0.48\textwidth]{poissoncounthist100hz100ms}
% \titlecaption{\label{countstatsfig}Count Statistik.}{} % \titlecaption{\label{countstatsfig}Count Statistik.}{}
% \end{figure} % \end{figure}
The number of events $n_i$ (counts) in a time window $i$ of the duration $W$
Die Anzahl der Ereignisse $n_i$ in Zeifenstern $i$ der yields positive integer random numbers that are commonly quantified
L\"ange $W$ ergeben ganzzahlige, positive Zufallsvariablen, die meist using the following measures:
durch folgende Sch\"atzer charakterisiert werden:
\begin{itemize} \begin{itemize}
\item Histogramm der counts $n_i$. \item Histogram of the counts $n_i$.
\item Mittlere Anzahl von Ereignissen: $\mu_N = \langle n \rangle$. \item Average number of events: $\mu_N = \langle n \rangle$.
\item Varianz der Anzahl: $\sigma_n^2 = \langle (n - \langle n \rangle)^2 \rangle$. \item Variance of the counts: $\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 \determ{Fano Faktor} (The variance divided by the average): $F = \frac{\sigma_n^2}{\mu_n}$.
\end{itemize} \end{itemize}
Insbesondere ist die mittlere Rate der Ereignisse $r$ (Spikes pro And in particular the average firing rate $r$ (spike count per time interval
Zeit, \determ{Feuerrate}) gemessen in Hertz \sindex[term]{Feuerrate!mittlere Rate} , \determ{Feuerrate}) that is given in Hertz \sindex[term]{Feuerrate!mittlere Rate}
\begin{equation} \begin{equation}
\label{firingrate} \label{firingrate}
r = \frac{\langle n \rangle}{W} \; . r = \frac{\langle n \rangle}{W} \; .
@ -200,315 +200,333 @@ Zeit, \determ{Feuerrate}) gemessen in Hertz \sindex[term]{Feuerrate!mittlere Rat
% \end{figure} % \end{figure}
\begin{exercise}{counthist.m}{} \begin{exercise}{counthist.m}{}
Schreibe eine Funktion \code{counthist()}, die aus mehreren trials Implement a function \code{counthist()} that calculates and plots
von Spiketrains die Verteilung der Anzahl der Spikes in Fenstern the distribution of spike counts observed in a certain time
einer der Funktion \"ubergegebenen Breite bestimmt, das Histogramm window. The function should take two input arguments: (i) a
plottet und zur\"uckgibt. Jeder trial der Spiketrains ist ein Vektor \codeterm{cell-array} of vectors containing the spike times in
mit den Spikezeiten gegeben in Sekunden als Element in einem seconds observed in a number of trials and (ii) the duration of the
\codeterm{cell-array}. time window that is used to evaluate the counts.
\end{exercise} \end{exercise}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Homogener Poisson Prozess} \section{Homogeneous Poisson process}
F\"ur kontinuierliche Me{\ss}gr\"o{\ss}en ist die Normalverteilung
u.a. wegen dem Zentralen Grenzwertsatz die Standardverteilung. Eine The Gaussian distribution is, due to the central limit theorem, the
\"ahnliche Rolle spielt bei Punktprozessen der \determ{Poisson standard for continuous measures. The equivalent in the realm of point
Prozess}. processes is the \enterm{Poisson distribution}.
Beim \determ[Poisson Prozess!homogener]{homogenen Poisson Prozess} In a \enterm[Poisson process!homogeneous]{homogeneous Poisson process}
treten Ereignisse mit einer festen Rate $\lambda=\text{const.}$ auf the events occur at a fixed rate $\lambda=\text{const.}$ and are
und sind unabh\"angig von der Zeit $t$ und unabh\"angig von den independent of both the time $t$ and occurrence of previous events
Zeitpunkten fr\"uherer Ereignisse (\figref{hompoissonfig}). Die (\figref{hompoissonfig}). The probability of observing an even within a
Wahrscheinlichkeit zu irgendeiner Zeit ein Ereigniss in einem kleinen small time window of width $\Delta t$ is given by
Zeitfenster der Breite $\Delta t$ zu bekommen ist
\begin{equation} \begin{equation}
\label{hompoissonprob} \label{hompoissonprob}
P = \lambda \cdot \Delta t \; . P = \lambda \cdot \Delta t \; .
\end{equation} \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}{} \begin{exercise}{poissonspikes.m}{}
Schreibe eine Funktion \code{poissonspikes()}, die die Spikezeiten Implement a function \code{poissonspikes()} that uses a homogeneous
eines homogenen Poisson-Prozesses mit gegebener Rate in Hertz f\"ur Poisson process to generate events at a given rate for a certain
eine Anzahl von trials gegebener maximaler L\"ange in Sekunden in duration and a number of trials. The rate should be given in Hertz
einem \codeterm{cell-array} zur\"uckgibt. Benutze \eqnref{hompoissonprob} and the duration of the trials is given in seconds. The function
um die Spikezeiten zu bestimmen. 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} \end{exercise}
\begin{figure}[t] \begin{figure}[t]
\includegraphics[width=1\textwidth]{poissonraster100hz} \includegraphics[width=1\textwidth]{poissonraster100hz}
\titlecaption{\label{hompoissonfig}Rasterplot von Spikes eines homogenen \titlecaption{\label{hompoissonfig}Rasterplot of spikes of a
Poisson Prozesses mit $\lambda=100$\,Hz.}{} homogeneous Poisson process with a rate $\lambda=100$\,Hz.}{}
\end{figure} \end{figure}
\begin{figure}[t] \begin{figure}[t]
\includegraphics[width=0.45\textwidth]{poissonisihexp20hz}\hfill \includegraphics[width=0.45\textwidth]{poissonisihexp20hz}\hfill
\includegraphics[width=0.45\textwidth]{poissonisihexp100hz} \includegraphics[width=0.45\textwidth]{poissonisihexp100hz}
\titlecaption{\label{hompoissonisihfig}Interspikeintervallverteilungen \titlecaption{\label{hompoissonisihfig}Distribution of interspike intervals of two Poisson processes.}{}
zweier Poissonprozesse.}{}
\end{figure} \end{figure}
Der homogene Poissonprozess hat folgende Eigenschaften: The homogeneous Poisson process has the following properties:
\begin{itemize} \begin{itemize}
\item Die Intervalle $T$ sind exponentiell verteilt (\figref{hompoissonisihfig}): \item Intervals $T$ are exponentially distributed (\figref{hompoissonisihfig}):
\begin{equation} \begin{equation}
\label{poissonintervals} \label{poissonintervals}
p(T) = \lambda e^{-\lambda T} \; . p(T) = \lambda e^{-\lambda T} \; .
\end{equation} \end{equation}
\item Das mittlere Intervall ist $\mu_{ISI} = \frac{1}{\lambda}$ . \item The average interval is $\mu_{ISI} = \frac{1}{\lambda}$ .
\item Die Varianz der Intervalle ist $\sigma_{ISI}^2 = \frac{1}{\lambda^2}$ . \item The variance of the intervals is $\sigma_{ISI}^2 = \frac{1}{\lambda^2}$ .
\item Der Variationskoeffizient ist also immer $CV_{ISI} = 1$ . \item Thus, the coefficient of variation is always $CV_{ISI} = 1$ .
\item Die \determ[serielle Korrelationen]{seriellen Korrelationen} \item The serial correlation is $\rho_k =0$ for $k>0$, since the
$\rho_k =0$ f\"ur $k>0$, da das Auftreten der Ereignisse occurrence of an event is independent of all previous events. Such a
unabh\"angig von der Vorgeschichte ist. Ein solcher Prozess wird process is also called a \enterm{renewal process}.
auch \determ{Erneuerungsprozess} genannt (\enterm{renewal process}). \item The number of events $k$ within a temporal window of duration
\item Die Anzahl der Ereignisse $k$ innerhalb eines Fensters der $W$ is Poisson distributed:
L\"ange W ist \determ[Poisson-Verteilung]{Poissonverteilt}:
\[ P(k) = \frac{(\lambda W)^ke^{\lambda W}}{k!} \] \[ P(k) = \frac{(\lambda W)^ke^{\lambda W}}{k!} \]
(\figref{hompoissoncountfig}) (\figref{hompoissoncountfig})
\item Der \determ{Fano Faktor} ist immer $F=1$ . \item The Fano Faktor is always $F=1$ .
\end{itemize} \end{itemize}
\begin{exercise}{hompoissonspikes.m}{} \begin{exercise}{hompoissonspikes.m}{}
Schreibe eine Funktion \code{hompoissonspikes()}, die die Spikezeiten Implement a function \code{hompoissonspikes()} that uses a
eines homogenen Poisson-Prozesses mit gegebener Rate in Hertz f\"ur homogeneous Poisson process to generate spike events at a given rate
eine Anzahl von trials gegebener maximaler L\"ange in Sekunden in for a certain duration and a number of trials. The rate should be
einem \codeterm{cell-array} zur\"uckgibt. Benutze die exponentiell-verteilten given in Hertz and the duration of the trials is given in
Interspikeintervalle \eqnref{poissonintervals}, um die Spikezeiten zu erzeugen. 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} \end{exercise}
\begin{figure}[t] \begin{figure}[t]
\includegraphics[width=0.48\textwidth]{poissoncounthistdist100hz10ms}\hfill \includegraphics[width=0.48\textwidth]{poissoncounthistdist100hz10ms}\hfill
\includegraphics[width=0.48\textwidth]{poissoncounthistdist100hz100ms} \includegraphics[width=0.48\textwidth]{poissoncounthistdist100hz100ms}
\titlecaption{\label{hompoissoncountfig}Z\"ahlstatistik von Poisson Spikes.}{} \titlecaption{\label{hompoissoncountfig}Count statistics of Poisson
spiketrains.}{}
\end{figure} \end{figure}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Zeitabh\"angige Feuerraten} \section{Time-dependent firing rate}
Bisher haben wir station\"are Spiketrains betrachtet, deren Statistik So far we have discussed stationary spiketrains. The statistical properties
sich innerhalb der Analysezeit nicht ver\"andert (station\"are of these did not change within the observation time (stationary point
Punktprozesse). Meistens jedoch \"andert sich die Statistik der processes). Most commonly, however, this is not the case. A sensory
Spiketrains eines Neurons mit der Zeit. Z.B. kann ein sensorisches neuron, for example, might respond to a stimulus by modulating its
Neuron auf einen Reiz hin mit einer erh\"ohten Feuerrate antworten firing rate (non-stationary point process).
(nichtstation\"arer Punktprozess).
How the firing rate $r(t)$ changes over time is the most important
Wie die mittlere Anzahl der Spikes sich mit der Zeit ver\"andert, die measure, when analyzing non-stationary spike trains. The unit of the
\determ{Feuerrate} $r(t)$, ist die wichtigste Gr\"o{\ss}e bei firing rate is Hertz, i.e. the number of action potentials per
nicht-station\"aren Spiketrains. Die Einheit der Feuerrate ist Hertz, second. There are different ways to estimate the firing rate and three
also Anzahl Aktionspotentiale pro Sekunde. Es gibt verschiedene of these methods are illustrated in \figref{psthfig}. All of
Methoden diese zu bestimmen. Drei solcher Methoden sind in Abbildung these have their own justifications and pros- and cons. In the
\ref{psthfig} dargestellt. Alle Methoden haben ihre Berechtigung und following we will discuss these methods more
ihre Vor- und Nachteile. Im folgenden werden die drei Methoden aus closely.
Abbildung \ref{psthfig} n\"aher erl\"autert.
\begin{figure}[tp] \begin{figure}[tp]
\includegraphics[width=\columnwidth]{firingrates} \includegraphics[width=\columnwidth]{firingrates}
\titlecaption{Bestimmung der zeitabh\"angigen \titlecaption{Estimating the time-dependent firing
Feuerrate.}{\textbf{A)} Rasterplot eines Spiketrains. \textbf{B)} rate.}{\textbf{A)} Rasterplot depicting the spiking activity of a
Feurerrate aus der instantanen Feuerrate bestimmt. \textbf{C)} neuron. \textbf{B)} Firing rate calculated from the
klassisches PSTH mit der Binning Methode. \textbf{D)} Feuerrate \emph{instantaneous rate}. \textbf{C)} classical PSTH with the
durch Faltung mit einem Gauss Kern bestimmt.}\label{psthfig} \emph{binning} method. \textbf{D)} Firing rate estimated by
\emph{convolution} of the activity with a Gaussian
kernel.}\label{psthfig}
\end{figure} \end{figure}
\subsection{Instantane Feuerrate} \subsection{Instantaneous firing rate}
\begin{figure}[tp] \begin{figure}[tp]
\includegraphics[width=\columnwidth]{isimethod} \includegraphics[width=\columnwidth]{isimethod}
\titlecaption{Instantane Feuerrate.}{Skizze eines Spiketrains \titlecaption{Instantaneous firing rate.}{Sketch of the recorded
(oben). Die Pfeile zwischen aufeinanderfolgenden spiketrain (top). Arrows illustrate the interspike intervals and
Aktionspotentialen mit den Zahlen in Millisekunden illustrieren number give the intervals in milliseconds. The inverse of the
die Interspikeintervalle. Der Kehrwert des Interspikeintervalle interspike interval is the \emph{instantaneous firing
ergibt die instantane Feuerrate.}\label{instrate} rate}.}\label{instratefig}
\end{figure} \end{figure}
Ein sehr einfacher Weg, die zeitabh\"angige Feuerrate zu bestimmen ist A very simple method for estimating the time-dependent firing rate is
die sogenannte \determ[Feuerrate!instantane]{instantane Feuerrate} the \enterm[firing rate!instantaneous]{instantaneous firing rate}. The
(\enterm[firing rate!instantaneous]{instantaneous firing rate}). Dabei firing rate can be directly estimated as the inverse of the time
wird die Feuerrate aus dem Kehrwert der Interspikeintervalle, der Zeit between successive spikes, the interspike-interval
zwischen zwei aufeinander folgenden Aktionspotentialen (\figref{instratefig}).
(\figref{instrate} A), bestimmt. Die abgesch\"atzte Feuerrate
(\figref{instrate} B) ist g\"ultig f\"ur das gesammte \begin{equation}
Interspikeintervall. Diese Methode hat den Vorteil, dass sie sehr \label{instantaneousrateeqn}
einfach zu berechnen ist und keine Annahme \"uber eine relevante r_i = \frac{1}{T_i} .
Zeitskala (der Kodierung oder des Auslesemechanismus der \end{equation}
postsynaptischen Zelle) macht. $r(t)$ ist allerdings keine
kontinuierliche Funktion, die Spr\"unge in der Feuerrate k\"onnen The instantaneous rate $r_i$ is valid for the whole interspike
f\"ur manche Analysen nachteilig sein. Au{\ss}erdem wird die Feuerrate interval. The method has the advantage of being extremely easy to
nie gleich Null, auch wenn lange keine Aktionspotentiale generiert compute and that it does not make any assumptions about the relevant
wurden. 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}{} \begin{exercise}{instantaneousRate.m}{}
Implementiere die Absch\"atzung der Feuerrate auf Basis der Implement a function that computes the instantaneous firing
instantanen Feuerrate. Plotte die Feuerrate als Funktion der Zeit. rate. Plot the firing rate as a function of time.
\note{TODO: example data!!!}
\end{exercise} \end{exercise}
\subsection{Peri-Stimulus-Zeit-Histogramm} \subsection{Peri-stimulus-time-histogram}
W\"ahrend die Instantane Rate den Kehrwert der Zeit von einem bis zum While the \emph{instantaneous firing rate} uses the interspike
n\"achsten Aktionspotential misst, sch\"atzt das sogenannte interval, the \enterm{peri stimulus time histogram} (PSTH) uses the
\determ{Peri-Stimulus-Zeit-Histogramm} (\enterm{peri stimulus time spike count within observation windows of the duration $W$. It
histogram}, \determ[PSTH|see{Peri-Stimulus-Zeit-Histogramm}]{PSTH}) estimates the probability of observing a spike within that observation
die Wahrscheinlichkeit ab, zu einem Zeitpunkt Aktionspotentiale time. It tries to estimat the average rate in the limit of small
anzutreffen. Es wird versucht die mittlere Rate \eqnref{firingrate} im obersvation times \eqnref{psthrate}:
Grenzwert kleiner Beobachtungszeiten abzusch\"atzen:
\begin{equation} \begin{equation}
\label{psthrate} \label{psthrate}
r(t) = \lim_{W \to 0} \frac{\langle n \rangle}{W} \; , r(t) = \lim_{W \to 0} \frac{\langle n \rangle}{W} \; ,
\end{equation} \end{equation}
wobei die Anzahl $n$ der Aktionspotentiale, die im Zeitintervall where $\langle n \rangle$ is the across trial average number of action
$(t,t+W)$ aufgetreten sind, \"uber trials gemittelt wird. Eine solche potentials observed within the interval $(t, t+W)$. Such description
Rate enspricht der zeitabh\"angigen Rate $\lambda(t)$ des inhomogenen matches the time-dependent firing rate $\lambda(t)$ of an
Poisson-Prozesses. inhomogeneous Poisson process.
Das PSTH \eqnref{psthrate} kann entweder \"uber die Binning-Methode The firing probability can be estimated using the \emph{binning
oder durch Verfaltung mit einem Kern bestimmt werden. Beiden Methoden method} or by using \emph{kernel density estimations}. Both methods
gemeinsam ist die Notwendigkeit der Wahl einer zus\"atzlichen Zeitskala, make an assumption about the relevant observation time-scale ($W$ in
die der Beobachtungszeit $W$ in \eqnref{psthrate} entspricht. \eqnref{psthrate}).
\subsubsection{Binning-Methode} \subsubsection{Binning-method}
\begin{figure}[tp] \begin{figure}[tp]
\includegraphics[width=\columnwidth]{binmethod} \includegraphics[width=\columnwidth]{binmethod}
\titlecaption{Bestimmung des PSTH mit der Binning Methode.}{Der \titlecaption{Estimating the PSTH using the binning method.}{ The
gleiche Spiketrain wie in \figref{instrate}. Die grauen Linien same spiketrain as shown in \figref{instratefig} is used
markieren die Grenzen der Bins und die Zahlen geben die Anzahl der Spikes here. Vertical gray lines indicate the borders between adjacent
in jedem Bin an (oben). Die Feuerrate ergibt sich aus dem bins in which the number of action potentials is counted (red
mit der Binbreite normierten Zeithistogramm (unten).}\label{binpsth} numbers). The firing rate is then the histogram normalized to the
binwidth.}\label{binpsthfig}
\end{figure} \end{figure}
Bei der Binning-Methode wird die Zeitachse in gleichm\"aßige The binning method separates the time axis into regular bins of the
Abschnitte (Bins) eingeteilt und die Anzahl Aktionspotentiale, die in bin width $W$ and counts for each bin the number of observed action
die jeweiligen Bins fallen, gez\"ahlt (\figref{binpsth} A). Um diese potentials (\figref{binpsthfig} top). The resulting histogram is then
Z\"ahlungen in die Feuerrate umzurechnen muss noch mit der Binweite normalized with the bin width $W$ to yield the firing rate shown in
normiert werden. Das ist \"aquivalent zur Absch\"atzung einer the bottom trace of figure \ref{binpsthfig}. The above sketched
Wahrscheinlichkeitsdichte. Es kann auch die \code{hist()} Funktion zur process is equivalent to estimating the probability density. It is
Bestimmung des PSTHs verwendet werden. \sindex[term]{Feuerrate!Binningmethode} possible to estimate the PSTH using the \code{hist()} method
\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 The estimated firing rate is valid for the total duration of each
Wahl der Binweite anh\"angt. $r(t)$ ist also keine stetige bin. This leads to the step-like plot shown in
Funktion. Die Binweite bestimmt die zeitliche Aufl\"osung der \figref{binpsthfig}. $r(t)$ is thus not a contiunous function in
Absch\"atzung. \"Anderungen in der Feuerrate, die innerhalb eines Bins time. The binwidth defines the temporal resolution of the firing rate
vorkommen k\"onnen nicht aufgl\"ost werden. Mit der Wahl der Binweite estimation Changes that happen within a bin cannot be resolved. Thus
wird somit eine Annahme \"uber die relevante Zeitskala des Spiketrains chosing a bin width implies an assumption about the relevant
gemacht. time-scale.
\pagebreak[4] \pagebreak[4]
\begin{exercise}{binnedRate.m}{} \begin{exercise}{binnedRate.m}{}
Implementiere die Absch\"atzung der Feuerrate mit der ``binning'' Implement a function that estimates the firing rate using the
Methode. Plotte das PSTH. ``binning method''. The method should take the spike-times as an
input argument and returns the firing rate. Plot the PSTH.
\end{exercise} \end{exercise}
\subsubsection{Faltungsmethode} \subsubsection{Convolution method --- Kernel density estimation}
\begin{figure}[tp] \begin{figure}[tp]
\includegraphics[width=\columnwidth]{convmethod} \includegraphics[width=\columnwidth]{convmethod}
\titlecaption{Bestimmung des PSTH mit der Faltungsmethode.}{Der \titlecaption{Estimating the firing rate using the convolution
gleiche Spiketrain wie in \figref{instrate}. Bei der Verfaltung method.}{The same spiketrain as in \figref{instratefig}. The
des Spiketrains mit einem Faltungskern wird jeder Spike durch den convolution of the spiketrain with a convolution kernel basically
Faltungskern ersetzt (oben). Bei korrekter Normierung des replaces each spike event with the kernel (top). A Gaussian kernel
Kerns ergibt sich die Feuerrate direkt aus der \"Uberlagerung der is used here, but other kernels are also possible. If the kernel
Kerne.}\label{convrate} is properly normalized the firing rate results directly form the
superposition of the kernels.}\label{convratefig}
\end{figure} \end{figure}
Bei der Faltungsmethode werden die harten Kanten der Bins der With the convolution method we avoid the sharp edges of the binning
Binning-Methode vermieden. Der Spiketrain wird mit einem Kern method. The spiketrain is convolved with a \enterm{convolution
verfaltet, d.h. jedes Aktionspotential wird durch den Kern ersetzt. kernel}. Technically speaking we need to first create a binary
Zur Berechnung wird die Aktionspotentialfolge zun\"achst representation of the spike train. This binary representation is a
``bin\"ar'' dargestellt. Dabei wird ein Spiketrain als 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:
(Zeit-)Vektor dargestellt, in welchem die Zeitpunkte der
Aktionspotentiale als 1 notiert werden. Alle anderen Elemente des \[r(t) = \int_{-\infty}^{\infty} \omega(\tau) \, \rho(t-\tau) \, {\rm d}\tau \; , \]
Vektors sind 0. Anschlie{\ss}end wir dieser bin\"are Spiketrain mit where $\omega(\tau)$ represents the kernel and $\rho(t)$ the binary
einem Gau{\ss}-Kern bestimmter Breite verfaltet: representation of the response. The process of convolution can be
\[r(t) = \int_{-\infty}^{\infty} \omega(\tau) \, \rho(t-\tau) \, {\rm d}\tau \; , \] imagined as replacing each event of the spiketrain with the kernel
wobei $\omega(\tau)$ der Filterkern und $\rho(t)$ die bin\"are Antwort (figure \ref{convratefig} top). The superposition of the replaced
ist. Bildlich geprochen wird jede 1 in $\rho(t)$ durch den Filterkern kernels is then the firing rate (if the kerel is correctly normalized
ersetzt (Abbildung \ref{convrate} A). Wenn der Kern richtig normiert to an integral of one, figure \ref{convreffig}
wurde (Integral gleich Eins), ergibt sich die Feuerrate direkt aus der bottom). \sindex[term]{Feuerrate!Faltungsmethode}
\"Uberlagerung der Kerne (Abb. \ref{convrate} B). \sindex[term]{Feuerrate!Faltungsmethode}
In contrast to the other methods the convolution methods leads to a
Die Faltungsmethode f\"uhrt, anders als die anderen Methoden, zu einer continuous function which is often desirable (in particular when
stetigen Funktion was insbesondere f\"ur spektrale Analysen von applying methods in the frequency domain). The choice of the kernel
Vorteil sein kann. Die Wahl der Kernbreite bestimmt, \"ahnlich zur width defines, similar to the bin width of the binning method, the
Binweite, die zeitliche Aufl\"osung von $r(t)$. Die Breite des Kerns temporal resolution of the method and thus makes assumptions about the
macht also auch wieder eine Annahme \"uber die relevante Zeitskala des relevate time-scale.
Spiketrains.
\pagebreak[4] \pagebreak[4]
\begin{exercise}{convolutionRate.m}{} \begin{exercise}{convolutionRate.m}{}
Verwende die Faltungsmethode um die Feuerrate zu bestimmen. Plotte Implement the function that estimates the firing rate using the
das Ergebnis. 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} \end{exercise}
\section{Spike-triggered Average} \section{Spike-triggered Average}
Die graphischer Darstellung der Feuerrate allein reicht nicht aus um
den Zusammenhang zwischen neuronaler Antwort und einem Stimulus zu The graphical representation of the neuronal activity alone is not
analysieren. Eine Methode um mehr \"uber diesen Zusammenhang zu sufficient tot investigate the relation between the neuronal response
erfahren, ist der \enterm{spike-triggered average} and a stimulus. One method to do this is the \enterm[STA|see
(\enterm[STA|see{spike-triggered average}]{STA}). Der STA spike-triggered average]{spike-triggered average}. The STA
\begin{equation} \begin{equation}
STA(\tau) = \langle s(t - \tau) \rangle = \frac{1}{N} \sum_{i=1}^{N} s(t_i - \tau) STA(\tau) = \langle s(t - \tau) \rangle = \frac{1}{N} \sum_{i=1}^{N} s(t_i - \tau)
\end{equation} \end{equation}
der $N$ Aktionspotentiale zu den Zeiten $t_i$ in Anwort auf den of $N$ action potentials observed at the times $t_i$ in response to
Stimulus $s(t)$ ist der mittlere Stimulus, der zu einem the stimulus $s(t)$ is the average stimulus that led to a spike in the
Aktionspotential in der neuronalen Antwort f\"uhrt. neuron. The STA can be easily extracted by cutting snippets out of
$s(t)$ that surround the times of the spikes. The resulting stimulus
Der STA l\"a{\ss}t sich relativ einfach berechnen, indem aus dem snippets are then averaged (\figref{stafig}).
Stimulus f\"ur jeden beobachteten Spike ein entsprechender Abschnitt
ausgeschnitten wird und diese dann gemittelt werde (\figref{stafig}).
\begin{figure}[t] \begin{figure}[t]
\includegraphics[width=\columnwidth]{sta} \includegraphics[width=\columnwidth]{sta}
\titlecaption{Spike-triggered Average eines P-Typ Elektrorezeptors \titlecaption{Spike-triggered average of a P-type electroreceptors
und Stimulusrekonstruktion.}{Der STA (links): der Rezeptor and the stimulus reconstruction.}{The neuron was driven by a
wurde mit einem ``white-noise'' Stimulus getrieben. Zeitpunkt 0 ``white-noise'' stimulus (blue curve, right panel). The STA (left)
ist der Zeitpunkt des beobachteten Aktionspotentials. Die Kurve is the average stimulus that surrounds the times of the recorded
ergibt sich aus dem gemittelten Stimulus je 50\,ms vor und nach action potentials (40\,ms before and 20\,ms after the
einem Aktionspotential. Stimulusrekonstruktion mittels spike). Using the STA as a convolution kernel for convolving the
STA (rechts). Die Zellantwort wird mit dem STA gefaltet um eine spiketrain we can reconstruct the stimulus from the neuronal
Rekonstruktion des Stimulus zu erhalten.}\label{stafig} 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} \end{figure}
Aus dem STA k\"onnen verschiedene Informationen \"uber den From the STA we can extract several pieces of information about the
Zusammenhang zwischen Stimulus und neuronaler Antwort gewonnen relation of stimulus and response. The width of the STA represents the
werden. Die Breite des STA repr\"asentiert die zeitliche Pr\"azision, temporal precision with which the neuron encodes the stimulus waveform
mit der Stimulus und Antwort zusammenh\"angen und wie weit das Neuron and in which temporal window the neuron integrates the (sensory)
zeitlich integriert. Die Amplitude des STA (gegeben in der gleichen input. The amplitude of the STA tells something about the sensitivity
Einheit wie der Stimulus) deutet auf die Empfindlichkeit des Neurons of the neuron. The STA is given in the same units as the stimulus and
bez\"uglich des Stimulus hin. Eine hohe Amplitude besagt, dass es a small amplitude indicates that the neuron needs only a small
einen starken Stimulus ben\"otigt, um ein Aktionspotential stimulus amplitude to create a spike, a large amplitude on the
hervorzurufen. Aus dem zeitlichen Versatz des STA kann die Zeit contrary suggests the opposite. The temporal delay between the STA and
abgelesen werden, die das System braucht, um auf den Stimulus zu the time of the spike is a consequence of the time the system (neuron)
antworten. needs to process the stimulus.
Der STA kann auch dazu benutzt werden, aus den Antworten der Zelle den We can further use the STA to do a \enterm{reverse reconstruction} and
Stimulus zu rekonstruieren (\figref{stafig} B). Bei der estimate the stimulus from the neuronal response (\figref{stafig},
\determ[invertierte Rekonstruktion]{invertierten Rekonstruktion} wird right). For this, the spiketrain is convolved with the STA as a
die Zellantwort mit dem STA verfaltet. kernel.
\begin{exercise}{spikeTriggeredAverage.m}{} \begin{exercise}{spikeTriggeredAverage.m}{}
Implementiere eine Funktion, die den STA ermittelt. Verwende dazu Implement a function that calculates the STA. Use the dataset
den Datensatz \file{sta\_data.mat}. Die Funktion sollte folgende \file{sta\_data.mat}. The function expects the spike train, the
R\"uckgabewerte haben: stimulus and the temporal resolution of the recording as input
arguments and should return the following information:
\vspace{-1ex} \vspace{-1ex}
\begin{itemize} \begin{itemize}
\setlength{\itemsep}{0ex} \setlength{\itemsep}{0ex}
\item den Spike-Triggered-Average. \item the spike-triggered average.
\item die Standardabweichung der individuellen STAs. \item the standard deviation of the STA across the individual snippets.
\item die Anzahl Aktionspotentiale, die zur Berechnung des STA verwendet wurden. \item The number of action potentials used to estimate the STA.
\vspace{-2ex}
\end{itemize} \end{itemize}
\end{exercise} \end{exercise}
\begin{exercise}{reconstructStimulus.m}{} \begin{exercise}{reconstructStimulus.m}{}
Rekonstruiere den Stimulus mithilfe des STA und der Spike Do the reverse reconstruction using the STA and the spike times. The
Zeiten. Die Funktion soll Vektor als R\"uckgabewert haben, der function should return the estimated stimulus in a vector that has
genauso gro{\ss} ist wie der Originalstimulus aus der Datei the same size as the original stimulus contained in file
\file{sta\_data.mat}. \file{sta\_data.mat}.
\end{exercise} \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 = [] times = []
v = vreset v = vreset
noise = np.sqrt(2.0*D)*np.random.randn(len(input))/np.sqrt(dt) 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 v += (input[k]+noise[k])*dt/tau
if v >= vthresh : if v >= vthresh :
v = vreset v = vreset
@ -55,7 +55,7 @@ rng = np.random.RandomState(54637281)
time = np.arange(0.0, duration, dt) time = np.arange(0.0, duration, dt)
x = np.zeros(time.shape)+rate x = np.zeros(time.shape)+rate
n = rng.randn(len(time))*drate*tau/np.sqrt(dt)+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[k] = x[k-1] + (n[k]-x[k-1])*dt/tau
x[x<0.0] = 0.0 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 = [] times = []
v = vreset v = vreset
noise = np.sqrt(2.0*D)*np.random.randn(len(input))/np.sqrt(dt) 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 v += (input[k]+noise[k])*dt/tau
if v >= vthresh : if v >= vthresh :
v = vreset v = vreset
@ -41,7 +41,7 @@ def pifspikes(input, trials, dt, D=0.1) :
def isis( spikes ) : def isis( spikes ) :
isi = [] isi = []
for k in xrange(len(spikes)) : for k in range(len(spikes)) :
isi.extend(np.diff(spikes[k])) isi.extend(np.diff(spikes[k]))
return np.array( isi ) return np.array( isi )
@ -84,7 +84,7 @@ rng = np.random.RandomState(54637281)
time = np.arange(0.0, duration, dt) time = np.arange(0.0, duration, dt)
x = np.zeros(time.shape)+rate x = np.zeros(time.shape)+rate
n = rng.randn(len(time))*drate*tau/np.sqrt(dt)+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[k] = x[k-1] + (n[k]-x[k-1])*dt/tau
x[x<0.0] = 0.0 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 = [] times = []
v = vreset v = vreset
noise = np.sqrt(2.0*D)*np.random.randn(len(input))/np.sqrt(dt) 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 v += (input[k]+noise[k])*dt/tau
if v >= vthresh : if v >= vthresh :
v = vreset v = vreset
@ -41,7 +41,7 @@ def pifspikes(input, trials, dt, D=0.1) :
def isis( spikes ) : def isis( spikes ) :
isi = [] isi = []
for k in xrange(len(spikes)) : for k in range(len(spikes)) :
isi.extend(np.diff(spikes[k])) isi.extend(np.diff(spikes[k]))
return np.array( isi ) return np.array( isi )
@ -95,7 +95,7 @@ rng = np.random.RandomState(54637281)
time = np.arange(0.0, duration, dt) time = np.arange(0.0, duration, dt)
x = np.zeros(time.shape)+rate x = np.zeros(time.shape)+rate
n = rng.randn(len(time))*drate*tau/np.sqrt(dt)+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[k] = x[k-1] + (n[k]-x[k-1])*dt/tau
x[x<0.0] = 0.0 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): def plot_sta(times, stim, dt, t_min=-0.1, t_max=.1):
count = 0 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) time = np.arange(t_min, t_max, dt)
if len(stim.shape) > 1 and stim.shape[1] > 1: if len(stim.shape) > 1 and stim.shape[1] > 1:
stim = stim[:,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+bitand(10, 8)+
\part Execute and explain: \verb+bitor(10, 8)+ \part Execute and explain: \verb+bitor(10, 8)+
\end{parts} \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} \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. \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{parts}
\end{questions} \end{questions}
@ -82,7 +84,7 @@ matrices that match in certain criteria. This process is called
the following commands and explain. the following commands and explain.
\begin{parts} \begin{parts}
\part \verb+x < 5+ \part \verb+x < 5+
\part \verb+x( x < 5) )+ \part \verb+x( (x < 5) )+
\part \verb+x( (y <= 2) )+ \part \verb+x( (y <= 2) )+
\part \verb+x( (x > 2) | (y < 8) )+ \part \verb+x( (x > 2) | (y < 8) )+
\part \verb+x( (x == 0) & (y == 0) )+ \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. \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). \part Count the number of elements in each class using Boolean expressions (\verb+sum+ can be used to count the matches).
\end{parts} \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{questions}
\end{document} \end{document}

47
programming/exercises/control_flow.tex
View File

@ -14,7 +14,7 @@
%%%%% text size %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%% text size %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\usepackage[left=20mm,right=20mm,top=25mm,bottom=25mm]{geometry} \usepackage[left=20mm,right=20mm,top=25mm,bottom=25mm]{geometry}
\pagestyle{headandfoot} \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: \firstpagefooter{Dr. Jan Grewe}{Phone: 29 74588}{Email:
jan.grewe@uni-tuebingen.de} jan.grewe@uni-tuebingen.de}
\runningfooter{}{\thepage}{} \runningfooter{}{\thepage}{}
@ -33,39 +33,44 @@
\begin{center} \begin{center}
\textbf{\Large Introduction to scientific computing}\\[1ex] \textbf{\Large Introduction to scientific computing}\\[1ex]
{\large Jan Grewe, Jan Benda}\\[-3ex] {\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} \end{center}
The exercises are meant for self-monitoring, revision of the lecture The exercises are meant for self-monitoring and revision of the
topic. You should try to solve them on your own. Your solution should lecture. You should try to solve them on your own. Your solution
be submitted as a single script (m-file) in the Ilias system. Each should be submitted as a single script (m-file) in the Ilias
task should be solved in its own ``cell''. Each cell must be system. Each task should be solved in its own ``cell''. Each cell must
executable on its own. The file should be named according to the following pattern: be executable on its own. The file should be named according to the
``variables\_datatypes\_\{lastname\}.m'' benannt werden following pattern: ``variables\_datatypes\_\{lastname\}.m''
(e.g. variables\_datentypes\_mueller.m). (e.g. variables\_datentypes\_mueller.m).
\begin{questions} \begin{questions}
\question Implement \code{for} loops in which the \emph{running variable}: \question Implement \code{for} loops in which the \emph{running variable}:
\begin{parts} \begin{parts}
\part ... assumes values from 0 to 10. Print the value of the running variable for each iteration step. \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 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. \part ... assumes values from 0 to 1 in steps of 0.1. Print out the value for each iteration.
\end{parts} \end{parts}
\question Indexing in vectros: \question Indexing in vectors:
\begin{parts} \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. \part Use a \code{for} loop to print each element of \code{x}. Use the running variable for the indexing.
\begin{solution} \begin{solution}
\code{for i = 1:length(x); disp(x(i)); end} \code{for i = 1:length(x); disp(x(i)); end}
\end{solution} \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} \begin{solution}
\code{for i : x; disp(i); end;} \code{for i : x; disp(i); end;}
\end{solution} \end{solution}
\end{parts} \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} \begin{parts}
\part Use a loop to calculate the arithmetic mean. The mean is defined as: \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 $. $\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 :-). \part Search the MATLAB help for functions that do the calculations for you :-).
\end{parts} \end{parts}
\question Implement a\code{while} loop \question Implement a \code{while} loop
\begin{parts} \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}. \part ... iterates endlessly. You can interrupt the execution with the command \code{Strg + C}.
\end{parts} \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 \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} \begin{parts}
\part Count the number of required iterations. \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 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 schows for each of the 1000 trials, how often you had to draw. \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? \part Play around with the threshold. What happens?
\end{parts} \end{parts}
@ -100,19 +105,19 @@ executable on its own. The file should be named according to the following patte
\end{parts} \end{parts}
\question Test the random number generator! To do so count the \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 [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 loop to draw 1000 random numbers with \verb+rand()+ (see help). What
would be the expectation? What is the result? 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} \begin{parts}
\part Create a variable \part Create a variable
\verb+filename = '2015-10-12_100Hz_1.25V.dat'+. Obviously, the \verb+filename = '2015-10-12_100Hz_1.25V.dat'+. Obviously, the
underscore was used as a delimiter. 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 \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 \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. \part Print out the individual parts.
\end{parts} \end{parts}

12
programming/exercises/scripts_functions.tex
View File

@ -16,7 +16,7 @@
%%%%% text size %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%% text size %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\usepackage[left=20mm,right=20mm,top=25mm,bottom=25mm]{geometry} \usepackage[left=20mm,right=20mm,top=25mm,bottom=25mm]{geometry}
\pagestyle{headandfoot} \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: \firstpagefooter{Prof. Jan Benda}{Phone: 29 74 573}{Email:
jan.benda@uni-tuebingen.de} jan.benda@uni-tuebingen.de}
\runningfooter{}{\thepage}{} \runningfooter{}{\thepage}{}
@ -46,7 +46,6 @@
aboveskip=10pt aboveskip=10pt
} }
\newcommand{\code}[1]{\texttt{#1}} \newcommand{\code}[1]{\texttt{#1}}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
@ -60,7 +59,7 @@
\end{center} \end{center}
The exercises are meant for self-monitoring and revision of the 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 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 but each question needs an individual file. Combine the files into a
single zip archive and submit it via ILIAS. Name the archive according single zip archive and submit it via ILIAS. Name the archive according
@ -108,7 +107,7 @@ to the pattern: ``scripts\_functions\_\{surname\}.zip''.
\part{} \part{}
Improve the function that it takes the duration of the time axis, 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 calculation should use a temporal stepsize that is 0.01 of the
frequency. frequency.
\begin{solution} \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 \question Random Walk. In a 1-D random walk an \emph{agent} walks
randomly either in the one ($+1$) or the other ($-1$) randomly either in the one ($+1$) or the other ($-1$)
direction. With each simulation step one direction is chosen and the 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} \begin{parts}
\part{} \part{}
@ -172,7 +172,7 @@ to the pattern: ``scripts\_functions\_\{surname\}.zip''.
Now we want to know how the probability of $p_{+1}$ (the Now we want to know how the probability of $p_{+1}$ (the
probability to walk into the $+1$ direction) impacts the random 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 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). stopping the simulations as before).
\begin{solution} \begin{solution}
\lstinputlisting{randomwalkscriptc.m} \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[german]{babel}
\usepackage{natbib} \usepackage{natbib}

24
programming/lecture/programming.tex
View File

@ -40,14 +40,12 @@ variable.
\begin{figure} \begin{figure}
\centering \centering
\begin{subfigure}{.5\textwidth} \begin{subfigure}{.5\textwidth}
\includegraphics[width=0.8\textwidth]{variable} \includegraphics[width=0.8\textwidth]{variable}\label{variable:a}
\label{variable:a}
\end{subfigure}% \end{subfigure}%
\begin{subfigure}{.5\textwidth} \begin{subfigure}{.5\textwidth}
\includegraphics[width=.8\textwidth]{variableB} \includegraphics[width=.8\textwidth]{variableB}\label{variable:b}
\label{variable:b}
\end{subfigure} \end{subfigure}
\titlecaption{Variables} point to a memory \titlecaption{Variables}{ point to a memory
address. They further are described by their name and address. They further are described by their name and
data type. The variable's value is stored as a pattern of binary data type. The variable's value is stored as a pattern of binary
values (0 or 1). When reading the variable this pattern is 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. matrices.
To request the length of a vector we used the function 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, information about the size of a matrix. As mentioned above,
\code{length()} would return the length of the largest dimension. The \code{length()} would return the length of the largest dimension. The
function \code{size()} however, returns the length in each dimension 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{x} represent the measured data at the times in
\varcode{t}. \varcode{t}.
\item Use logical indexing to select those values that have been \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{itemize}
\end{exercise} \end{exercise}
@ -1095,7 +1093,7 @@ ans =
Generally, a program is executed line by line from top to Generally, a program is executed line by line from top to
bottom. Sometimes this is not the desired behavior, or the other way bottom. Sometimes this is not the desired behavior, or the other way
round, it is needed to skip certain parts or execute others 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 statements that allow to manipulate the control flow. There are two
major classes of such statements: major classes of such statements:
@ -1120,7 +1118,7 @@ x =
120 120
\end{lstlisting} \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 thing that changes is the increasing factor. The repetition of such
very similar lines of code is bad programming style. This is not only 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: 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. repetitions. It is easy to forget a single change.
\item Readability: repetitive code is terrible to read and to \item Readability: repetitive code is terrible to read and to
understand. (I) one tends to skip repetitions (its the same, 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. code leads to long and hard-to-parse programs.
\end{enumerate} \end{enumerate}
All imperative programming languages offer a solution: the loop. It is 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} \subsubsection{The \code{for} --- loop}
The most common type of loop is the \codeterm{for-loop}. It The most common type of loop is the \codeterm{for-loop}. It
consists of a \codeterm[Loop!head]{head} and the 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 body is executed. In \matlab{} the head begins with the keyword
\code{for} which is followed by the \codeterm{running variable}. In \code{for} which is followed by the \codeterm{running variable}. In
\matlab{} a for-loop always operates on vectors. With each \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 % execute this code if the expression yields true
end end
\end{lstlisting} \end{lstlisting}
\begin{exercise}{factorialWhileLoop.m}{} \begin{exercise}{factorialWhileLoop.m}{}
Implement the factorial of a number \varcode{n} using a \code{while} Implement the factorial of a number \varcode{n} using a \code{while}
-- loop. -- loop.
@ -1201,7 +1199,7 @@ end
\begin{exercise}{neverendingWhile.m}{} \begin{exercise}{neverendingWhile.m}{}
Implement a \code{while}--loop that is never-ending. Hint: the body Implement a \code{while}--loop that is never-ending. Hint: the body
is executed as long as the Boolean expression in the head is 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} \end{exercise}

8
programmingstyle/code/calculateSines.m
View File

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

157
programmingstyle/lecture/programmingstyle.tex
View File

@ -6,20 +6,19 @@
%\selectlanguage{ngerman} %\selectlanguage{ngerman}
Cultivating a good code style not a matter of good taste but is 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 facilitates reproducibility of scientific results. Programs should be
written and structured in a way that supports outsiders as well the written and structured in a way that supports outsiders as well the
author himself --- a few weeks or months after it was written --- to author himself --- a few weeks or months after it was written --- to
understand the programs' rationale. Clean code pays off for the understand the programs' rationale. Clean code pays off for the
original author as well as others that are supposed to use the code. original author as well as others that are supposed to use the code.
Clean code addresses several issues: Clean code addresses several issues:
\begin{enumerate} \begin{enumerate}
\item The programs' structure. \item The programs' structure.
\item Naming of scripts and functions. \item Naming of scripts and functions.
\item Naming of variables and constants. \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 Use of comments and inline documentation.
\item Delegation of repeated code to functions and dedicated \item Delegation of repeated code to functions and dedicated
subroutines. subroutines.
@ -29,13 +28,12 @@ Clean code addresses several issues:
While introducing scripts and functions we suggested a typical program While introducing scripts and functions we suggested a typical program
layout (box\,\ref{whenscriptsbox}). The idea is to create a single layout (box\,\ref{whenscriptsbox}). The idea is to create a single
entry point by having one script that controls the rest of program by entry point by having one script that controls the rest of the program
managing data and results and calling functions that work on the data by calling functions that work on the data and managing the
and produce the results. Applying this structure makes it easy to results. Applying this structure makes it easy to understand the flow
understand the flow of the program but two questions remain: (i) How of the program but two questions remain: (i) How to organize the files
to organize the files on the file system and (ii) how to name them on the file system and (ii) how to name them that the controlling
that the controlling script is easily identified among the other script is easily identified among the other \codeterm{m-files}.
\codeterm{m-files}.
Upon installation ``MATLAB'' creates a folder called \emph{MATLAB} in Upon installation ``MATLAB'' creates a folder called \emph{MATLAB} in
the user space (Windows: My files, Linux: Documents, MacOS: the user space (Windows: My files, Linux: Documents, MacOS:
@ -110,24 +108,11 @@ Box~\ref{matlabpathbox}).
\end{lstlisting} \end{lstlisting}
\end{ibox} \end{ibox}
\section{Naming scripts and functions} \section{Naming things}
\matlab{} will search the search path (Box \ref{matlabpathbox}) The dictum of good code style is: ``Program code must be readable.''
exclusively by name. It is case-sensitive this implies that the files Expressive names are extraordinarily important in this respect. Even
\file{test\_function.m} and \file{Test\_function.m} are two different if it is tricky to find expressive names that are not overly long,
things. It is self-evident that choosing such names is nonsensical naming should be taken seriously.
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}
\matlab{} has a few rules about names: Names must not start with a \matlab{} has a few rules about names: Names must not start with a
number, they must not contain blanks or other special characters like 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 There are other common patterns such as the \emph{camelCase} in which
the first character of compound words is capitalized. Other the first character of compound words is capitalized. Other
conventions use the underscore to separate the individual words conventions use the underscore to separate the individual words (
\emph{snake\_case}. A function that counts the number of action \emph{snake\_case}). A function that counts the number of action
potentials could be named \file{spikeCount.m} or potentials could be named \file{spikeCount.m} or
\file{spike\_count.m}. \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 \subsection{Naming variables and constants}
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.
While the names of scripts and functions describe the purpose, names While the names of scripts and functions describe the purpose, names
of variables describe the stored content. A variable storing the 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 \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}). (\varcode{average\_spike\_counts}).
The control variables used in the head of a \code{for} loop are often 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 does not enforce indentation, indentation is very powerful for
defining coherent blocks. The \matlab{} editor supports this by an defining coherent blocks. The \matlab{} editor supports this by an
auto-indentation mechanism. A selected section of the code and be 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 Interspersing empty lines is very helpful to separate regions in the
code that belong together. Too many empty lines, however lead to 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. by the screen and thus takes overview.
The following two listings show basically the same implementation of a The following two listings show basically the same implementation of a
random walk once in a rather chaotic version (listing random walk\footnote{A random walk is a simple simulation of Brownian
\ref{chaoticcode}) then in cleaner way (listing \ref{cleancode}) 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.}] \begin{lstlisting}[label=chaoticcode, caption={Chaotic implementation of the random-walk.}]
num_runs = 10; max_steps = 1000; num_runs = 10; max_steps = 1000;
@ -245,25 +248,24 @@ end
\section{Using comments} \section{Using comments}
It is common to provide extra information about the meaning of program It is common to provide extra information about the meaning of program
code by adding comments to it. In \matlab{} comments are indicated by code by adding comments. In \matlab{} comments are indicated by the
the percent character \code{\%}. Anything that is written in the percent character \code{\%}. Anything that follows the percent
respective line following the percent is ignored and considered a character in a line is ignored and considered a comment. When used
comment. When used sparsely comments can immensely important for sparsely comments can be immensely helpful. Comments
understanding. Comments are short sentences that describe the meaning are short sentences that describe the meaning of the (following) lines
of the (following) lines in the program code. During the initial in the program code. During the initial implementation of a function
implementation of a function they can be used to guide the development they can be used to guide the development but have the tendency to
but have the tendency to blow up the code and decrease readability. By blow up the code and decrease readability. By choosing expressive
choosing expressive variable and function names, most lines should be variable and function names, most lines should be self-explanatory.
self-explanatory.
For example stating the obvious does not really help and should be
For example stating the obvious does not really help:\\ avoided:\\ \varcode{ x = x + 2; \% add two to x}\\
\varcode{ x = x + 2; \% add two to x}\\
\begin{important}[Using comments] \begin{important}[Using comments]
\begin{itemize} \begin{itemize}
\item Comments describe the rationale of the respective code block. \item Comments describe the rationale of the respective code block.
\item Comments are good and helpful --- they have to be true, however! \item Comments are good and helpful --- they must be true, however!
\item A wrong comment is worse than a non-existent comment! \item A wrong comment is worse than a non-existent one!
\item Comments must be maintained just as the code. Otherwise they \item Comments must be maintained just as the code. Otherwise they
may become wrong and worse than meaningless! may become wrong and worse than meaningless!
\end{itemize} \end{itemize}
@ -303,7 +305,7 @@ well documented function.
Comments and empty lines are used to organize code into logical blocks 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 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 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 Not delegating the tasks leads to very long \codeterm{m-files} which
can be confusing. Sometimes such a code is called ``spaghetti 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 the same name as the function itself. Delegating tasks to functions
thus leads to a large set of \codeterm{m-files} which increases thus leads to a large set of \codeterm{m-files} which increases
complexity and may lead to confusion. If the delegated functionality complexity and may lead to confusion. If the delegated functionality
is used in multiple instances, it is advisable to do so. On the other is used in multiple instances, it is still advisable to do so. On the
hand, when the delegated functionality is only used within the context other hand, when the delegated functionality is only used within the
of another function \matlab{} allows to define context of another function \matlab{} allows to define
\codeterm[function!local]{local functions} and \codeterm[function!local]{local functions} and
\codeterm[function!nested]{nested functions} within the same \codeterm[function!nested]{nested functions} within the same
file. Listing \ref{localfunctions} shows an example of a local file. Listing \ref{localfunctions} shows an example of a local
@ -336,17 +338,19 @@ function definition.
\pagebreak[3] \lstinputlisting[label=localfunctions, caption={Example \pagebreak[3] \lstinputlisting[label=localfunctions, caption={Example
for local functions.}]{calculateSines.m} for local functions.}]{calculateSines.m}
Local function live in the same \codeterm{m-file} as the main function \emph{Local function} live in the same \codeterm{m-file} as the main
and are only available in this context. Each local function has its function and are only available in this context. Each local function
own \codeterm{scope}, that is, the local function can not access (read has its own \codeterm{scope}, that is, the local function can not
or write) variables of the calling function. 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 \emph{Nested functions} are different in this respect. They are
functions}. These are defined within the body of the parent function defined within the body of the parent function (between the keywords
(between the keywords \code{function} and \code{end}) and have full \code{function} and \code{end}) and have full access to all variables
access to all variables defined in the parent function. Working (in defined in the parent function. Working (in particular changing) the
particular changing) the parent's variables is handy on the one side, parent's variables is handy on the one side, but is also risky. One
but is also risky. One should take care when defining nested functions. should take care when defining nested functions.
\section{Specifics when using scripts} \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 become available in the global \codeterm{Workspace}. There is the risk
of name conflicts, that is, a called sub-script redefines or uses the of name conflicts, that is, a called sub-script redefines or uses the
same variable name and may \emph{silently} change its content. 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 expect a completely different content. Bugs that are based on such
mistakes are hard to find since the program itself looks perfectly mistakes are hard to find since the program itself looks perfectly
fine. fine.
@ -379,12 +383,15 @@ provides important information to track and fix the bug.
\item Scripts should work independently of existing variables in the \item Scripts should work independently of existing variables in the
global workspace. global workspace.
\item It is advisable to start a script with deleting variables \item Often it is advisable to start a script with deleting
(\code{clear}) from the workspace and most of the times it is also variables (\code{clear}) from the workspace and most of the times
good to close all open figures (\code{close all}). 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 \item Clean up the workspace at the end of a script. Delete
(\code{clear}) all variables that are no longer needed. (\code{clear}) all variables that are no longer needed.
\item Consider to write functions instead of scripts.
\end{itemize} \end{itemize}
\end{important} \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 $f_{cost}(\{(x_i, y_i)\}|\{y^{est}_i\})$ is a so called
\enterm{objective function} or \enterm{cost function}. We aim to adapt \enterm{objective function} or \enterm{cost function}. We aim to adapt
the model parameters to minimize the error (mean square error) and 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 we will show that the minimization of the mean square error is
equivalent to maximizing the likelihood that the observations equivalent to maximizing the likelihood that the observations
originate from the model (assuming a normal distribution of the data 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$ (Box~\ref{partialderivativebox}) with respect to the parameters $m$
and $b$ of the linear equation. There is no need to calculate it and $b$ of the linear equation. There is no need to calculate it
analytically but it can be estimated from the partial derivatives 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 small steps $\Delta m$ und $\Delta b$. For example the partial
derivative with respect to $m$: 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 Problems that involve nonlinear computations on parameters, e.g. the
rate $\lambda$ in the exponential function $f(x;\lambda) = 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 in such functions numerical methods such as the gradient descent have
to be applied. to be applied.
@ -396,7 +396,7 @@ objective functions while more specialized functions are specifically
designed for optimizations in the least square error sense designed for optimizations in the least square error sense
\matlabfun{lsqcurvefit()}. \matlabfun{lsqcurvefit()}.
\newpage %\newpage
\begin{important}[Beware of secondary minima!] \begin{important}[Beware of secondary minima!]
Finding the absolute minimum is not always as easy as in the case of 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 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} \lstset{inputpath=regression/code}
\include{regression/lecture/regression} \include{regression/lecture/regression}
\setboolean{showexercisesolutions}{true} \setboolean{showexercisesolutions}{false}
\graphicspath{{likelihood/lecture/}{likelihood/lecture/figures/}} \graphicspath{{likelihood/lecture/}{likelihood/lecture/figures/}}
\lstset{inputpath=likelihood/code} \lstset{inputpath=likelihood/code}
\include{likelihood/lecture/likelihood} \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), arrowprops=dict(arrowstyle="->", relpos=(0.5,0.0),
connectionstyle="angle3,angleA=-60,angleB=80") ) 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['right'].set_visible(False)
ax.spines['top'].set_visible(False) ax.spines['top'].set_visible(False)
ax.spines['left'].set_visible(False) ax.spines['left'].set_visible(False)
@ -92,7 +92,7 @@ ax.annotate('median',
arrowprops=dict(arrowstyle="->", relpos=(0.8,0.0), arrowprops=dict(arrowstyle="->", relpos=(0.8,0.0),
connectionstyle="angle3,angleA=-60,angleB=20") ) 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['right'].set_visible(False)
ax.spines['top'].set_visible(False) ax.spines['top'].set_visible(False)
ax.xaxis.set_ticks_position('bottom') ax.xaxis.set_ticks_position('bottom')