diff --git a/bootstrap/lecture/bootstrapsem.py b/bootstrap/lecture/bootstrapsem.py index 347f40f..66154ce 100644 --- a/bootstrap/lecture/bootstrapsem.py +++ b/bootstrap/lecture/bootstrapsem.py @@ -15,13 +15,13 @@ x = rng.randn(nsamples) # bootstrap the mean: mus = [] -for i in xrange(nresamples) : +for i in range(nresamples) : mus.append(np.mean(x[rng.randint(0, nsamples, nsamples)])) hmus, _ = np.histogram(mus, bins, density=True) # many SRS: musrs = [] -for i in xrange(nresamples) : +for i in range(nresamples) : musrs.append(np.mean(rng.randn(nsamples))) hmusrs, _ = np.histogram(musrs, bins, density=True) diff --git a/bootstrap/lecture/permutecorrelation.py b/bootstrap/lecture/permutecorrelation.py index ab1a583..b6bc1bd 100644 --- a/bootstrap/lecture/permutecorrelation.py +++ b/bootstrap/lecture/permutecorrelation.py @@ -19,7 +19,7 @@ rd = np.corrcoef(x, y)[0, 1] # permutation: nperm = 1000 rs = [] -for i in xrange(nperm) : +for i in range(nperm) : xr=rng.permutation(x) yr=rng.permutation(y) rs.append( np.corrcoef(xr, yr)[0, 1] ) diff --git a/chapter.mk b/chapter.mk index 8a0e56e..d9aa865 100644 --- a/chapter.mk +++ b/chapter.mk @@ -8,7 +8,8 @@ PYPDFFILES=$(PYFILES:.py=.pdf) pythonplots : $(PYPDFFILES) $(PYPDFFILES) : %.pdf: %.py - python $< + echo $$(which python) + python3 $< cleanpythonplots : rm -f $(PYPDFFILES) diff --git a/likelihood/lecture/mlemean.py b/likelihood/lecture/mlemean.py index 892e663..5fd4932 100644 --- a/likelihood/lecture/mlemean.py +++ b/likelihood/lecture/mlemean.py @@ -41,7 +41,7 @@ for mu in mus : ax.text(mu-0.1, 0.04, '?', zorder=1, ha='right') else : ax.text(mu+0.1, 0.04, '?', zorder=1) -for k in xrange(len(mus)) : +for k in range(len(mus)) : ax.plot(x, g[:,k], zorder=5) ax.scatter(xd, 0.05*rng.rand(len(xd))+0.2, s=30, zorder=10) diff --git a/plotting/code/simple_plot.m b/plotting/code/simple_plot.m new file mode 100644 index 0000000..71bb549 --- /dev/null +++ b/plotting/code/simple_plot.m @@ -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') diff --git a/plotting/lecture/plotting.tex b/plotting/lecture/plotting.tex index 07e3ae4..f631189 100644 --- a/plotting/lecture/plotting.tex +++ b/plotting/lecture/plotting.tex @@ -11,123 +11,6 @@ results. severe.}{\url{www.xkcd.com}}\label{xkcdplotting} \end{figure} -\section{What makes a good plot?} -Plot should help/enable the interested reader to get a grasp of the -data and to understand the performed analysis and to critically assess -the presented results. The most important rule is the correct and -complete annotation of the plots. This starts with axis labels and -units and and extends to legends. Incomplete annotation can have -terrible consequences (\figref{xkcdplotting}). - -The principle of \emph{ink minimization} may be used a a guiding -principle for appealing plots. It requires that the relation of amount -of ink spent on the data and that spent on other parts of the plot -should be strongly in favor of the data. Ornamental of otherwise -unnecessary gimicks should not be used in scientific contexts. An -exception can be made if the particular figure was designed for -didactic purposes and sometimes for presentations. - -\begin{important}[Correct labeling of plots] - A data plot must be sufficiently labeled: - \begin{itemize} - \item Every axis must have a label and the correct unit, if it has - one.\\ (e.g. \code[xlabel()]{xlabel('Speed [m/s]'}). - \item When more than one line is plotted, they have to be labeled - using the figure legend, or similar \matlabfun{legend()}. - \item If using subplots that show similar information on the axes, - they should be scaled to show the same ranges to ease comparison - between plots. (e.g. \code[xlim()]{xlim([0 100])}.\\ If one - chooses to ignore this rule one should explicitly state this in - the figure caption and/or the descriptions in the text. - \item Labels must be large enough to be readable. In particular, - when using the figure in a presentation use large enough fonts. - \end{itemize} -\end{important} - - -\section{Things that should be avoided.} - -When plotting scientific data we should take great care to avoid -suggestive or misleading presentations. Unnecessary additions and -fancy graphical effects make a plot frivolous and also violate the -\emph{ink minimization principle}. Illustrations in comic style -(\figref{comicexamplefig}) are not suited for scientific data in most -instances. For presentations or didactic purposes, however, using a -comic style may be helpful to indicate that the figure is a mere -sketch and the exact position of the data points is of no importance. - -\begin{figure}[t] - \includegraphics[width=0.7\columnwidth]{outlier}\vspace{-3ex} - \titlecaption{Comic-like illustration.}{Obviously not suited to - present scientific data. In didactic or illustrative contexts they - can be helpful to focus on the important - aspects.}\label{comicexamplefig} -\end{figure} - - -The following figures show examples of misleading or suggestive -presentations of data. Several of the effects have been exaggerated to -make the point. A little more subtlety these methods are employed to -nudge the viewers experience into the desired direction. You can find -more examples on \url{https://en.wikipedia.org/wiki/Misleading_graph}. - -\begin{figure}[p] - \includegraphics[width=0.35\textwidth]{misleading_pie} - \hspace{0.05\textwidth} - \includegraphics[width=0.35\textwidth]{sample_pie} - \titlecaption{Perspective distortion influences the perceived - size.}{By changing the perspective of the 3-D illustration the - highlighted segment \textbf{C} gains more weight than it should - have. In the left graph segments \textbf{A} and \textbf{C} appear - very similar. The 2-D plot on the right-hand side shows that this - is an - illusion. \url{https://en.wikipedia.org/wiki/Misleading_graph}}\label{misleadingpiefig} -\end{figure} - -\begin{figure}[p] - \includegraphics[width=0.9\textwidth]{plot_scaling.pdf} - \titlecaption{Choosing the figure format and scaling of the axes - influences the perceived strength of a correlation.}{All subplots - show the same data. By choosing a certain figure size we can - pronounce or reduce the perceived strength of the correlation - in the data. Technically all three plots are correct. - }\label{misleadingscalingfig} -\end{figure} - -\begin{figure}[p] - \begin{minipage}[t]{0.3\textwidth} - \includegraphics[width=0.8\textwidth]{improperly_scaled_graph} - \end{minipage} - \begin{minipage}[t]{0.3\textwidth} - \includegraphics[width=0.8\textwidth]{comparison_properly_improperly_graph} - \end{minipage} - \begin{minipage}[t]{0.3\textwidth} - \includegraphics[width=0.7\textwidth]{properly_scaled_graph} - \end{minipage} - \titlecaption{Scaling of markers and symbols.} {In these graphs - symbols have been used to illustrate the measurements made in two - categories. The measured value for category \textbf{B} is actually - three times the measured value for category \textbf{A}. In the - left graph the symbol for category \textbf{B} has been scaled to - triple height while maintaining the proportions. This appears just - fair and correct but leads to the effect that the covered surface - is not increased to the 3-fold but the 9-fold (center plot). The - plot on the right shows how it could have been done correctly. - \url{https://en.wikipedia.org/wiki/Misleading_graph}}\label{misleadingsymbolsfig} -\end{figure} - -By using perspective effects in 3-D plot the perceived size can be -distorted into the desired direction. While the plot is correct in a -strict sense it is rather suggestive -(\figref{misleadingpiefig}). Similarly the choice of figure size and -proportions can lead to different interpretations of the -data. Stretching the y-extent of a graph leads to a stronger -impression of the correlation in the data. Compressing this axis will -lead to a much weaker perceived correlation -(\figref{misleadingscalingfig}). When using symbols to illustrate a -quantity we have to take care not to overrate of difference due to -symbol scaling (\figref{misleadingsymbolsfig}). - \section{The \matlab{} plotting system} Plotting data in \matlab{} is rather straight forward for simple line @@ -191,19 +74,20 @@ the data was changed or the same kind of plot has to be created for a number of datasets. \begin{important}[Why manual editing should be avoided.] - On first glance the manual editing of a figure using common tools - like Corel draw, Illustrator, etc.\,appears some much more - convenient and less complex. This, however, is not entirely - true. What if the figure has to be re-drawn or updated? Then the - editing work starts all over again. Rather, there is a great risk - associated with this approach. Axes are shifted, fonts have not been - embedded into the final document, annotations have been copy pasted - between figures and are not valid. All of these mistakes can be - found in publications and then require an erratum, which is not + On first glance the manual editing of a figure using tools such as + Corel draw, Illustrator, etc.\,appears much more convenient and less + complex than coding everything into the analysis scripts. This, + however, is not entirely true. What if the figure has to be re-drawn + or updated? Then the editing work starts all over again. Rather, + there is a great risk associated with the manual editing + approach. Axes may be shifted, fonts have not been embedded into the + final document, annotations have been copy pasted between figures + and are not valid. All of these mistakes can be found in + publications and then require an erratum, which is not desirable. Even if it appears more cumbersome in the beginning one should always try to create publication-ready figures directly from - the data analysis tool using scripts or functions to properly - layout the plot. + the data analysis tool using scripts or functions to properly layout + the plot. \end{important} \subsection{Simple plotting} @@ -259,6 +143,12 @@ additional options consult the help. \end{tabular} \end{table} +The following listing shows a simple line plot with axis labeling and a title + +\lstinputlisting[caption={A simple plot showing a sinewave.}, + label=simpleplotlisting]{simple_plot.m} + + \subsection{Changing properties of a line plot} The properties of line plots can be changed by passing more arguments @@ -419,14 +309,14 @@ output format (box\,\ref{graphicsformatbox}). \end{tabular} \end{minipage} - It is often meaningful to store of data plots generated by \matlab{} - using a vector graphics format. When in doubt they can usually be + It is advisable to store of data plots generated by \matlab{} + using a vector graphics format. In doubt they can usually be easily converted to a bitmap format. The way from a bitmap to a vector graphic is not possible without a loss in quality. Storing a - plot that contains a very large set of graphical elements (e.g.\,a + plot that contains very large sets of graphical elements (e.g.\,a raster-plot showing thousands of action potentials) may, on the other hand, lead to very large files that can be hard to - handle. Saving such a plot using a bitmap format may be more + handle. Saving such plots using a bitmap format may be more efficient. \end{ibox} @@ -584,24 +474,24 @@ its properties. See the \matlab{} help for more information. \begin{figure}[ht] \includegraphics[width=0.9\linewidth]{errorbars} - \titlecaption{Adding error bars to a line plot}{\textbf{A} + \titlecaption{Indicating the estimation error in plots.}{\textbf{A} symmetrical error around the mean (e.g.\ using the standard deviation). \textbf{B} Errorbars of an asymmetrical distribution of the data (note: the average value is now the median and the errors are the lower and upper quartiles). \textbf{C} A shaded area is used to illustrate the spread of the data. See - listing\,\ref{errorbarlisting}}\label{errorbarplot} + listing\,\ref{errorbarlisting} for A and C and listing\,\ref{errorbarlisting2} }\label{errorbarplot} \end{figure} -\lstinputlisting[caption={Illustrating estimation errors. Script that - creates \figref{errorbarplot}.}, +\lstinputlisting[caption={Illustrating estimation errors using error bars. Script that + creates \figref{errorbarplot}. A, B}, label=errorbarlisting, firstline=13, lastline=29, basicstyle=\ttfamily\scriptsize]{errorbarplot.m} \subsubsection{Fill} For a few years now it has become fancy to illustrate the error not using errorbars but by drawing a shaded area around the mean. Beside -their fancyness there is also a real argument in favor of using error +the fancyness there is also a real argument in favor of using error areas instead of errorbars: In case you have a lot of data points with respective errorbars such that they would merge in the figure it is cleaner and probably easier to read and handle if one uses an error @@ -613,8 +503,8 @@ with the vertex points of the polygon. For each x-value we now have two y-values (average minus error and average plus error). Further, we want the vertices to be connected in a defined order. One can achieve this by going back and forth on the x-axis; we append a reversed -version of the x-values to the original x-values using the \code{cat} -and inversion is done using the \code{fliplr} command (line 3 in +version of the x-values to the original x-values using \code{cat} and +\code{fliplr} for concatenation and inversion, respectively (line 3 in listing \ref{errorbarlisting2}; Depending on the layout of your data you may need concatenate along a different dimension of the data and use \code{flipud} instead). The y-coordinates of the polygon vertices @@ -625,27 +515,26 @@ property defines the transparency (or rather the opaqueness) of the area. The provided alpha value is a number between 0 and 1 with zero leading to invisibility and a value of one to complete opaqueness. Finally, we use the normal plot command to draw a line -connecting the average values. +connecting the average values (line 12). -\lstinputlisting[caption={Illustrating estimation errors. Script that - creates \figref{errorbarplot}.}, label=errorbarlisting2, +\lstinputlisting[caption={Illustrating estimation errors using a shaded area. Script that + creates \figref{errorbarplot} C.}, label=errorbarlisting2, firstline=30, basicstyle=\ttfamily\scriptsize]{errorbarplot.m} \subsection{Annotations, text} -Sometimes want to highlight certain parts of a plot or simply add an -annotation that does not fit or belong to the legend. In these cases -we can use the \code[text()]{text()} or -\code[annotation()]{annotation()} function to add this information to -the plot. While \varcode{text} simply prints out the given text string -at the defined position (for example line in +The \code[text()]{text()} or \code[annotation()]{annotation()} are +used for highlighting certain parts of a plot or simply adding an +annotation that does not fit or does not belong into the legend. +While \varcode{text} simply prints out the given text string at the +defined position (for example line in listing\,\ref{regularsubplotlisting}) the \varcode{annotation} function allows to add some more advanced highlights like arrows, lines, ellipses, or rectangles. Figure\,\ref{annotationsplot} shows some examples, the respective code can be found in listing\,\ref{annotationsplotlisting}. For more options consult the -documentation. +\matlab{} help. \begin{figure}[ht] \includegraphics[width=0.5\linewidth]{annotations} @@ -714,20 +603,133 @@ Lissajous figure. The basic steps are: movie.}, label=animationlisting, firstline=3, lastline=33, basicstyle=\ttfamily\scriptsize]{movie_example.m} +\section{What makes a good plot?} +Plot should help/enable the interested reader to get a grasp of the +data and to understand the performed analysis and to critically assess +the presented results. The most important rule is the correct and +complete annotation of the plots. This starts with axis labels and +units and and extends to legends. Incomplete annotation can have +terrible consequences (\figref{xkcdplotting}). + +The principle of \emph{ink minimization} may be used a a guiding +principle for appealing plots. It requires that the relation of amount +of ink spent on the data and that spent on other parts of the plot +should be strongly in favor of the data. Ornamental of otherwise +unnecessary gimicks should not be used in scientific contexts. An +exception can be made if the particular figure was designed for +didactic purposes and sometimes for presentations. + +\begin{important}[Correct labeling of plots] + A data plot must be sufficiently labeled: + \begin{itemize} + \item Every axis must have a label and the correct unit, if it has + one.\\ (e.g. \code[xlabel()]{xlabel('Speed [m/s]'}). + \item When more than one line is plotted, they have to be labeled + using the figure legend, or similar \matlabfun{legend()}. + \item If using subplots that show similar information on the axes, + they should be scaled to show the same ranges to ease comparison + between plots. (e.g. \code[xlim()]{xlim([0 100])}.\\ If one + chooses to ignore this rule one should explicitly state this in + the figure caption and/or the descriptions in the text. + \item Labels must be large enough to be readable. In particular, + when using the figure in a presentation use large enough fonts. + \end{itemize} +\end{important} + + +\section{Things that should be avoided.} + +When plotting scientific data we should take great care to avoid +suggestive or misleading presentations. Unnecessary additions and +fancy graphical effects make a plot frivolous and also violate the +\emph{ink minimization principle}. Illustrations in comic style +(\figref{comicexamplefig}) are not suited for scientific data in most +instances. For presentations or didactic purposes, however, using a +comic style may be helpful to indicate that the figure is a mere +sketch and the exact position of the data points is of no importance. + +\begin{figure}[t] + \includegraphics[width=0.7\columnwidth]{outlier}\vspace{-3ex} + \titlecaption{Comic-like illustration.}{Obviously not suited to + present scientific data. In didactic or illustrative contexts they + can be helpful to focus on the important + aspects.}\label{comicexamplefig} +\end{figure} + + +The following figures show examples of misleading or suggestive +presentations of data. Several of the effects have been exaggerated to +make the point. A little more subtlety these methods are employed to +nudge the viewers experience into the desired direction. You can find +more examples on \url{https://en.wikipedia.org/wiki/Misleading_graph}. + +\begin{figure}[p] + \includegraphics[width=0.35\textwidth]{misleading_pie} + \hspace{0.05\textwidth} + \includegraphics[width=0.35\textwidth]{sample_pie} + \titlecaption{Perspective distortion influences the perceived + size.}{By changing the perspective of the 3-D illustration the + highlighted segment \textbf{C} gains more weight than it should + have. In the left graph segments \textbf{A} and \textbf{C} appear + very similar. The 2-D plot on the right-hand side shows that this + is an + illusion. \url{https://en.wikipedia.org/wiki/Misleading_graph}}\label{misleadingpiefig} +\end{figure} +\begin{figure}[p] + \includegraphics[width=0.9\textwidth]{plot_scaling.pdf} + \titlecaption{Choosing the figure format and scaling of the axes + influences the perceived strength of a correlation.}{All subplots + show the same data. By choosing a certain figure size we can + pronounce or reduce the perceived strength of the correlation + in the data. Technically all three plots are correct. + }\label{misleadingscalingfig} +\end{figure} +\begin{figure}[p] + \begin{minipage}[t]{0.3\textwidth} + \includegraphics[width=0.8\textwidth]{improperly_scaled_graph} + \end{minipage} + \begin{minipage}[t]{0.3\textwidth} + \includegraphics[width=0.8\textwidth]{comparison_properly_improperly_graph} + \end{minipage} + \begin{minipage}[t]{0.3\textwidth} + \includegraphics[width=0.7\textwidth]{properly_scaled_graph} + \end{minipage} + \titlecaption{Scaling of markers and symbols.} {In these graphs + symbols have been used to illustrate the measurements made in two + categories. The measured value for category \textbf{B} is actually + three times the measured value for category \textbf{A}. In the + left graph the symbol for category \textbf{B} has been scaled to + triple height while maintaining the proportions. This appears just + fair and correct but leads to the effect that the covered surface + is not increased to the 3-fold but the 9-fold (center plot). The + plot on the right shows how it could have been done correctly. + \url{https://en.wikipedia.org/wiki/Misleading_graph}}\label{misleadingsymbolsfig} +\end{figure} +By using perspective effects in 3-D plot the perceived size can be +distorted into the desired direction. While the plot is correct in a +strict sense it is rather suggestive +(\figref{misleadingpiefig}). Similarly the choice of figure size and +proportions can lead to different interpretations of the +data. Stretching the y-extent of a graph leads to a stronger +impression of the correlation in the data. Compressing this axis will +lead to a much weaker perceived correlation +(\figref{misleadingscalingfig}). When using symbols to illustrate a +quantity we have to take care not to overrate of difference due to +symbol scaling (\figref{misleadingsymbolsfig}). \section{Summary} A good plot of scientific data displays the data completely and seriously without too many distractions. Misleading or suggestive plots as may result from perspective presentations, inappropriate -scaling of axes of symbols should be avoided. +scaling of axes and symbols should be avoided. \noindent When combining several line plots within the same figure one should consider adapting color \textbf{and} line style (solid, dashed, dotted. etc.) to make the distinguishable even in black-and-white -prints. Combinations of red and green are no good choice since they +prints. Combinations of red and green are not a good choice since they cannot be distinguished by people with red-green blindness. \vspace{2ex} @@ -737,6 +739,8 @@ Key ingredients for a good data plot: \item Complete labeling. \item Plotted lines and curves must be distinguishable. \item No suggestive or misleading presentation. -\item The right balance of line width, font size and size of the figure. +\item The right balance of line width, font size and size of the + figure, this may depend on the purpose, for presentations slightly + thicker lines help. \item Error bars wherever they are appropriate. \end{itemize} diff --git a/pointprocesses/lecture/firingrates.py b/pointprocesses/lecture/firingrates.py index 6b98732..d6d1666 100644 --- a/pointprocesses/lecture/firingrates.py +++ b/pointprocesses/lecture/firingrates.py @@ -86,7 +86,7 @@ def plot_isi_rate(spike_times, max_t=30, dt=1e-4): def get_binned_rate(times, bin_width=0.05, max_t=30., dt=1e-4): time = np.arange(0., max_t, dt) bins = np.arange(0., max_t, bin_width) - bin_indices = bins / dt + bin_indices = np.asarray(bins / dt, np.int) hist, _ = sp.histogram(times, bins) rate = np.zeros(time.shape) diff --git a/pointprocesses/lecture/isihexamples.py b/pointprocesses/lecture/isihexamples.py index b03d2a9..7d8d30e 100644 --- a/pointprocesses/lecture/isihexamples.py +++ b/pointprocesses/lecture/isihexamples.py @@ -31,7 +31,7 @@ def pifspikes(input, trials, dt, D=0.1) : times = [] v = vreset noise = np.sqrt(2.0*D)*np.random.randn(len(input))/np.sqrt(dt) - for k in xrange(len(noise)) : + for k in range(len(noise)) : v += (input[k]+noise[k])*dt/tau if v >= vthresh : v = vreset @@ -41,7 +41,7 @@ def pifspikes(input, trials, dt, D=0.1) : def isis( spikes ) : isi = [] - for k in xrange(len(spikes)) : + for k in range(len(spikes)) : isi.extend(np.diff(spikes[k])) return isi @@ -76,7 +76,7 @@ rng = np.random.RandomState(54637281) time = np.arange(0.0, duration, dt) x = np.zeros(time.shape)+rate n = rng.randn(len(time))*drate*tau/np.sqrt(dt)+rate -for k in xrange(1,len(x)) : +for k in range(1,len(x)) : x[k] = x[k-1] + (n[k]-x[k-1])*dt/tau x[x<0.0] = 0.0 diff --git a/pointprocesses/lecture/pointprocesses.tex b/pointprocesses/lecture/pointprocesses.tex index 4b0850f..1a78dfe 100644 --- a/pointprocesses/lecture/pointprocesses.tex +++ b/pointprocesses/lecture/pointprocesses.tex @@ -1,186 +1,186 @@ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -\chapter{Analyse von Spiketrains} +\chapter{Spiketrain analysis} -\selectlanguage{ngerman} +\selectlanguage{english} -\determ[Aktionspotential]{Aktionspotentiale} (\enterm{spikes}) sind die Tr\"ager der -Information in Nervensystemen. Dabei ist in erster Linie nur der -Zeitpunkt des Auftretens eines Aktionspotentials von Bedeutung. Die -genaue Form des Aktionspotentials spielt keine oder nur eine -untergeordnete Rolle. +\enterm[Actionspotentials]{Actionspotentials} (\enterm{spikes}) are +the carriers of information in the nervous system. Thereby it is +mainly the time at which the spikes are generated that is of +importance. The waveform of the action potential is largely +stereotyped and does not carry information. -Nach etwas Vorverarbeitung haben elektrophysiologische Messungen -deshalb Listen von Spikezeitpunkten als Ergebniss --- sogenannte -\enterm{spiketrains}. Diese Messungen k\"onnen wiederholt werden und -es ergeben sich mehrere \enterm{trials} von Spiketrains -(\figref{rasterexamplesfig}). +The result of the processing of electrophysiological recordings are +series of spike times, which are then termed \enterm{spiketrains}. If +measurements are repeated we yield several \enterm{trials} of +spiketrains (\figref{rasterexamplesfig}). -Spiketrains sind Zeitpunkte von Ereignissen --- den Aktionspotentialen ---- und deren Analyse f\"allt daher in das Gebiet der Statistik von -sogenannten \determ[Punktprozess]{Punktprozessen}. +Spiketrains are times of events, the action potentials. The analysis +of these leads into the realm of the so called \enterm[point + process]{point processes}. \begin{figure}[ht] \includegraphics[width=1\textwidth]{rasterexamples} - \titlecaption{\label{rasterexamplesfig}Raster-Plot.}{Raster-Plot von - jeweils 10 Realisierungen eines station\"arenen Punktprozesses - (homogener Poisson Prozess mit Rate $\lambda=20$\;Hz, links) und - eines nicht-station\"aren Punktprozesses (perfect - integrate-and-fire Neuron getrieben mit Ohrnstein-Uhlenbeck - Rauschen mit Zeitkonstante $\tau=100$\,ms, rechts). Jeder - vertikale Strich markiert den Zeitpunkt eines Ereignisses. - Jede Zeile zeigt die Ereignisse eines trials.} + \titlecaption{\label{rasterexamplesfig}Raster-plot.}{Raster-plot of + ten realizations of a stationary point process (homogeneous point + process with a rate $\lambda=20$\;Hz, left) and an inhomogeneous + point process (perfect integrate-and-fire neuron dirven by + Ohrnstein-Uhlenbeck noise with a time-constant $\tau=100$\,ms, + right). Each vertical dash illustrates the time at which the + action potential was observed. Each line represents the event of + each trial.} \end{figure} %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -\section{Punktprozesse} +\section{Point processes} -Ein zeitlicher Punktprozess (\enterm{point process}) ist ein -stochastischer Prozess, der eine Abfolge von Ereignissen zu den Zeiten -$\{t_i\}$, $t_i \in \reZ$, generiert. +A temporal \enterm{point process} is a stochastic process that +generates a sequence of events at times $\{t_i\}$, $t_i \in +\reZ$. -\begin{ibox}{Beispiele von Punktprozessen} -Jeder Punktprozess wird durch einen sich in der Zeit kontinuierlich -entwickelnden Prozess generiert. Wann immer dieser Prozess eine -Schwelle \"uberschreitet wird ein Ereigniss des Punktprozesses -erzeugt. Zum Beispiel: -\begin{itemize} -\item Aktionspotentiale/Herzschlag: wird durch die Dynamik des - Membranpotentials eines Neurons/Herzzelle erzeugt. -\item Erdbeben: wird durch die Dynamik des Druckes zwischen - tektonischen Platten auf beiden Seiten einer geologischen Verwerfung - erzeugt. -\item Zeitpunkt eines Grillen/Frosch/Vogelgesangs: wird durch die - Dynamik des Nervensystems und des Muskelapparates erzeugt. -\end{itemize} +\begin{ibox}{Examples of point processes} + Every point process is generated by a temporally continuously + developing process. An event is generated whenever this process + reaches a certain threshold. For example: + \begin{itemize} + \item Action potentials/heart beat: created by the dynamics of the + neuron/sinoatrial node + \item Earthquake: defined by the dynamics of the pressure between + tectonical plates. + \item Evoked communication calls in crickets/frogs/birds: shaped by + the dynamics of nervous system and the muscle appartus. + \end{itemize} \end{ibox} \begin{figure}[t] \texpicture{pointprocessscetch} - \titlecaption{\label{pointprocessscetchfig} Statistik von - Punktprozessen.}{Ein Punktprozess ist eine Abfolge von - Zeitpunkten $t_i$ die auch durch die Intervalle $T_i=t_{i+1}-t_i$ - oder die Anzahl der Ereignisse $n_i$ beschrieben werden kann. } + \titlecaption{\label{pointprocessscetchfig} Statistics of point + processes.}{A point process is a sequence of instances in time + $t_i$ that can be characterized through the inter-event-intervals + $T_i=t_{i+1}-t_i$ and the number of events $n_i$. } \end{figure} -F\"ur die Neurowissenschaften ist die Statistik der Punktprozesse -besonders wichtig, da die Zeitpunkte der Aktionspotentiale als -zeitlicher Punktprozess betrachtet werden k\"onnen und entscheidend -f\"ur die Informations\"ubertragung sind. - -Bei Punktprozessen k\"onnen wir die Zeitpunkte $t_i$ ihres Auftretens, -die Intervalle zwischen diesen Zeitpunkten $T_i=t_{i+1}-t_i$, sowie -die Anzahl der Ereignisse $n_i$ bis zu einer bestimmten Zeit betrachten -(\figref{pointprocessscetchfig}). +In the neurosciences, the statistics of point processes is of +importance since the timing of the neuronal events (the action +potentials) is crucial for information transmission and can be treated +as such a process. -Zwei Punktprozesse mit verschiedenen Eigenschaften sind in -\figref{rasterexamplesfig} als Rasterplot dargestellt, bei dem die -Zeitpunkte der Ereignisse durch senkrechte Striche markiert werden. +Point processes can be described using the intervals between +successive events $T_i=t_{i+1}-t_i$ and the number of observed events +within a certain time window $n_i$ (\figref{pointprocessscetchfig}). +The events originating from a point process can be illustrated in form +of a scatter- or raster plot in which each vertical line indicates the +time of an event. The event from two different point processes are +shown in \figref{rasterexamplesfig}. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -\section{Intervallstatistik} +\section{Intervalstatistics} -Die Intervalle $T_i=t_{i+1}-t_i$ zwischen aufeinanderfolgenden -Ereignissen sind reelle, positive Zahlen. Bei Aktionspotentialen -heisen die Intervalle auch \determ{Interspikeintervalle} -(\enterm{interspike intervals}). Deren Statistik kann mit den -\"ublichen Gr\"o{\ss}en beschrieben werden. +The intervals $T_i=t_{i+1}-t_i$ between successive events are real +positive numbers. In the context of action potentials they are +referred to as \enterm{interspike intervals}. The statistics of these +are described using the common measures. \begin{figure}[t] \includegraphics[width=0.96\textwidth]{isihexamples}\vspace{-2ex} - \titlecaption{\label{isihexamplesfig}Interspikeintervall Histogramme}{der in - \figref{rasterexamplesfig} gezeigten Spikes.} + \titlecaption{\label{isihexamplesfig}Interspike interval + histogram}{of the spikes depicted in \figref{rasterexamplesfig}.} \end{figure} \begin{exercise}{isis.m}{} - Schreibe eine Funktion \code{isis()}, die aus mehreren trials von Spiketrains die - Interspikeintervalle bestimmt und diese in einem Vektor - zur\"uckgibt. Jeder trial der Spiketrains ist ein Vektor mit den - Spikezeiten gegeben in Sekunden als Element in einem \codeterm{cell-array}. + Implement a function \code{isis()} that calculates the interspike + intervals from several spike trains. The function should return a + single vector of intervals. The action potentials recorded in the + individual trials are stored as vectors of spike times within a + \codeterm{cell-array}. Spike times are given in seconds. \end{exercise} -\subsection{Intervallstatistik erster Ordnung} +\subsection{First order interval statistics} \begin{itemize} -\item Wahrscheinlichkeitsdichte $p(T)$ der Intervalle $T$ - (\figref{isihexamplesfig}). Normiert auf $\int_0^{\infty} p(T) \; dT +\item Probability density $p(T)$ of the intervals $T$ + (\figref{isihexamplesfig}). Normalized to $\int_0^{\infty} p(T) \; dT = 1$. -\item Mittleres Intervall: $\mu_{ISI} = \langle T \rangle = +\item Average interval: $\mu_{ISI} = \langle T \rangle = \frac{1}{n}\sum\limits_{i=1}^n T_i$. -\item Standardabweichung der Intervalle: $\sigma_{ISI} = \sqrt{\langle (T - \langle T - \rangle)^2 \rangle}$\vspace{1ex} -\item \determ{Variationskoeffizient} (\enterm{coefficient of variation}): $CV_{ISI} = +\item Standard deviation of the interspike intervals: $\sigma_{ISI} = \sqrt{\langle (T - \langle T + \rangle)^2 \rangle}$\vspace{1ex} +\item \enterm{Coefficient of variation}: $CV_{ISI} = \frac{\sigma_{ISI}}{\mu_{ISI}}$. -\item \determ{Diffusionskoeffizient} (\enterm{diffusion coefficient}): $D_{ISI} = +\item \enterm{Diffusion coefficient}: $D_{ISI} = \frac{\sigma_{ISI}^2}{2\mu_{ISI}^3}$. \end{itemize} \begin{exercise}{isihist.m}{} - Schreibe eine Funktion \code{isiHist()}, die einen Vektor mit Interspikeintervallen - entgegennimmt und daraus ein normiertes Histogramm der Interspikeintervalle - berechnet. + Implement a function \code{isiHist()} that calculates the normalized + interspike interval histogram. The function should take two input + arguments; (i) a vector of interspike intervals and (ii) the width + of the bins used for the histogram. It further returns the + probability density as well as the centers of the bins. \end{exercise} \begin{exercise}{plotisihist.m}{} - Schreibe eine Funktion, die die Histogrammdaten der Funktion - \code{isiHist()} entgegennimmt, um das Histogramm zu plotten. Im - Plot sollen die Interspikeintervalle in Millisekunden aufgetragen - werden. Das Histogramm soll zus\"atzlich mit Mittelwert, - Standardabweichung und Variationskoeffizient der - Interspikeintervalle annotiert werden. + Implement a function that takes the return values of + \code{isiHist()} as input arguments and then plots the data. The + plot should show the histogram with the x-axis scaled to + milliseconds and should be annotated with the average ISI, the + standard deviation and the coefficient of variation. \end{exercise} -\subsection{Korrelationen der Intervalle} -In \enterm{return maps} werden die um das \enterm{lag} $k$ verz\"ogerten -Intervalle $T_{i+k}$ gegen die Intervalle $T_i$ geplottet. Dies macht -m\"ogliche Abh\"angigkeiten von aufeinanderfolgenden Intervallen -sichtbar. +\subsection{Interval correlations} +So called \enterm{return maps} are used to illustrate +interdependencies between successive interspike intervals. The return +map plots the delayed interval $T_{i+k}$ against the interval +$T_i$. The parameter $k$ is called the \enterm{lag} $k$. Stationary +and non-stationary return maps are distinctly different +\figref{returnmapfig}. \begin{figure}[t] \includegraphics[width=1\textwidth]{returnmapexamples} \includegraphics[width=1\textwidth]{serialcorrexamples} - \titlecaption{\label{returnmapfig}Interspikeintervall return maps und - serielle Korrelationen}{zwischen aufeinander folgenden Intervallen - im Abstand des Lags $k$.} + \titlecaption{\label{returnmapfig}Interspike interval analyses of a + stationary and a non-stationary pointprocess.}{Upper plots show the + return maps and the lower panels depict the serial correlation of + successive intervals separated by the lag $k$.} \end{figure} -Solche Ab\"angigkeiten werden durch die \determ{serielle - Korrelationen} (\enterm{serial correlations}) der Intervalle -quantifiziert. Das ist der \determ{Korrelationskoeffizient} zwischen -aufeinander folgenden Intervallen getrennt durch lag $k$: +Such dependencies can be further quantified using the \enterm{serial + correlations} \figref{returnmapfig}. The serial correlation is the +correlation coefficient of the intervals $T_i$ and the intervals +delayed by the lag $T_{i+k}$: \[ \rho_k = \frac{\langle (T_{i+k} - \langle T \rangle)(T_i - \langle T \rangle) \rangle}{\langle (T_i - \langle T \rangle)^2\rangle} = \frac{{\rm cov}(T_{i+k}, T_i)}{{\rm var}(T_i)} -= {\rm corr}(T_{i+k}, T_i) \] -\"Ublicherweise wird die Korrelation $\rho_k$ gegen den Lag $k$ -aufgetragen (\figref{returnmapfig}). $\rho_0=1$ (Korrelation jedes -Intervalls mit sich selber). += {\rm corr}(T_{i+k}, T_i) \] The resulting correlation coefficient +$\rho_k$ is usually plotted against the lag $k$ +\figref{returnmapfig}. $\rho_0=1$ is the correlation of each interval +with itself and is always 1. \begin{exercise}{isiserialcorr.m}{} - Schreibe eine Funktion \code{isiserialcorr()}, die einen Vektor mit Interspikeintervallen - entgegennimmt und daraus die seriellen Korrelationen berechnet und plottet. - \pagebreak[4] + Implement a function \code{isiserialcorr()} that takes a vector of + interspike intervals as input argument and calculates the serial + correlation. The function should further plot the serial + correlation. \pagebreak[4] \end{exercise} %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -\section{Z\"ahlstatistik} +\section{Count statistics} % \begin{figure}[t] % \includegraphics[width=0.48\textwidth]{poissoncounthist100hz10ms}\hfill % \includegraphics[width=0.48\textwidth]{poissoncounthist100hz100ms} % \titlecaption{\label{countstatsfig}Count Statistik.}{} % \end{figure} - -Die Anzahl der Ereignisse $n_i$ in Zeifenstern $i$ der -L\"ange $W$ ergeben ganzzahlige, positive Zufallsvariablen, die meist -durch folgende Sch\"atzer charakterisiert werden: +The number of events $n_i$ (counts) in a time window $i$ of the duration $W$ +yields positive integer random numbers that are commonly quantified +using the following measures: \begin{itemize} -\item Histogramm der counts $n_i$. -\item Mittlere Anzahl von Ereignissen: $\mu_N = \langle n \rangle$. -\item Varianz der Anzahl: $\sigma_n^2 = \langle (n - \langle n \rangle)^2 \rangle$. -\item \determ{Fano Faktor} (Varianz geteilt durch Mittelwert): $F = \frac{\sigma_n^2}{\mu_n}$. +\item Histogram of the counts $n_i$. +\item Average number of events: $\mu_N = \langle n \rangle$. +\item Variance of the counts: $\sigma_n^2 = \langle (n - \langle n \rangle)^2 \rangle$. +\item \determ{Fano Faktor} (The variance divided by the average): $F = \frac{\sigma_n^2}{\mu_n}$. \end{itemize} -Insbesondere ist die mittlere Rate der Ereignisse $r$ (Spikes pro -Zeit, \determ{Feuerrate}) gemessen in Hertz \sindex[term]{Feuerrate!mittlere Rate} +And in particular the average firing rate $r$ (spike count per time interval +, \determ{Feuerrate}) that is given in Hertz \sindex[term]{Feuerrate!mittlere Rate} \begin{equation} \label{firingrate} r = \frac{\langle n \rangle}{W} \; . @@ -200,315 +200,333 @@ Zeit, \determ{Feuerrate}) gemessen in Hertz \sindex[term]{Feuerrate!mittlere Rat % \end{figure} \begin{exercise}{counthist.m}{} - Schreibe eine Funktion \code{counthist()}, die aus mehreren trials - von Spiketrains die Verteilung der Anzahl der Spikes in Fenstern - einer der Funktion \"ubergegebenen Breite bestimmt, das Histogramm - plottet und zur\"uckgibt. Jeder trial der Spiketrains ist ein Vektor - mit den Spikezeiten gegeben in Sekunden als Element in einem - \codeterm{cell-array}. + Implement a function \code{counthist()} that calculates and plots + the distribution of spike counts observed in a certain time + window. The function should take two input arguments: (i) a + \codeterm{cell-array} of vectors containing the spike times in + seconds observed in a number of trials and (ii) the duration of the + time window that is used to evaluate the counts. \end{exercise} %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -\section{Homogener Poisson Prozess} -F\"ur kontinuierliche Me{\ss}gr\"o{\ss}en ist die Normalverteilung -u.a. wegen dem Zentralen Grenzwertsatz die Standardverteilung. Eine -\"ahnliche Rolle spielt bei Punktprozessen der \determ{Poisson - Prozess}. - -Beim \determ[Poisson Prozess!homogener]{homogenen Poisson Prozess} -treten Ereignisse mit einer festen Rate $\lambda=\text{const.}$ auf -und sind unabh\"angig von der Zeit $t$ und unabh\"angig von den -Zeitpunkten fr\"uherer Ereignisse (\figref{hompoissonfig}). Die -Wahrscheinlichkeit zu irgendeiner Zeit ein Ereigniss in einem kleinen -Zeitfenster der Breite $\Delta t$ zu bekommen ist +\section{Homogeneous Poisson process} + +The Gaussian distribution is, due to the central limit theorem, the +standard for continuous measures. The equivalent in the realm of point +processes is the \enterm{Poisson distribution}. + +In a \enterm[Poisson process!homogeneous]{homogeneous Poisson process} +the events occur at a fixed rate $\lambda=\text{const.}$ and are +independent of both the time $t$ and occurrence of previous events +(\figref{hompoissonfig}). The probability of observing an even within a +small time window of width $\Delta t$ is given by \begin{equation} \label{hompoissonprob} - P = \lambda \cdot \Delta t \; . + P = \lambda \cdot \Delta t \; . \end{equation} -Beim \determ[Poisson Prozess!inhomogener]{inhomogenen Poisson Prozess} -h\"angt die Rate $\lambda$ von der Zeit ab: $\lambda = \lambda(t)$. + +In an \enterm[Poisson process!inhomogeneous]{inhomogeneous Poisson + process}, however, the rate $\lambda$ depends on the time: $\lambda = +\lambda(t)$. \begin{exercise}{poissonspikes.m}{} - Schreibe eine Funktion \code{poissonspikes()}, die die Spikezeiten - eines homogenen Poisson-Prozesses mit gegebener Rate in Hertz f\"ur - eine Anzahl von trials gegebener maximaler L\"ange in Sekunden in - einem \codeterm{cell-array} zur\"uckgibt. Benutze \eqnref{hompoissonprob} - um die Spikezeiten zu bestimmen. + Implement a function \code{poissonspikes()} that uses a homogeneous + Poisson process to generate events at a given rate for a certain + duration and a number of trials. The rate should be given in Hertz + and the duration of the trials is given in seconds. The function + should return the event times in a cell-array. Each entry in this + array represents the events observed in one trial. Apply + \eqnref{hompoissonprob} to generate the event times. \end{exercise} \begin{figure}[t] \includegraphics[width=1\textwidth]{poissonraster100hz} - \titlecaption{\label{hompoissonfig}Rasterplot von Spikes eines homogenen - Poisson Prozesses mit $\lambda=100$\,Hz.}{} + \titlecaption{\label{hompoissonfig}Rasterplot of spikes of a + homogeneous Poisson process with a rate $\lambda=100$\,Hz.}{} \end{figure} \begin{figure}[t] \includegraphics[width=0.45\textwidth]{poissonisihexp20hz}\hfill \includegraphics[width=0.45\textwidth]{poissonisihexp100hz} - \titlecaption{\label{hompoissonisihfig}Interspikeintervallverteilungen - zweier Poissonprozesse.}{} + \titlecaption{\label{hompoissonisihfig}Distribution of interspike intervals of two Poisson processes.}{} \end{figure} -Der homogene Poissonprozess hat folgende Eigenschaften: +The homogeneous Poisson process has the following properties: \begin{itemize} -\item Die Intervalle $T$ sind exponentiell verteilt (\figref{hompoissonisihfig}): +\item Intervals $T$ are exponentially distributed (\figref{hompoissonisihfig}): \begin{equation} \label{poissonintervals} p(T) = \lambda e^{-\lambda T} \; . \end{equation} -\item Das mittlere Intervall ist $\mu_{ISI} = \frac{1}{\lambda}$ . -\item Die Varianz der Intervalle ist $\sigma_{ISI}^2 = \frac{1}{\lambda^2}$ . -\item Der Variationskoeffizient ist also immer $CV_{ISI} = 1$ . -\item Die \determ[serielle Korrelationen]{seriellen Korrelationen} - $\rho_k =0$ f\"ur $k>0$, da das Auftreten der Ereignisse - unabh\"angig von der Vorgeschichte ist. Ein solcher Prozess wird - auch \determ{Erneuerungsprozess} genannt (\enterm{renewal process}). -\item Die Anzahl der Ereignisse $k$ innerhalb eines Fensters der - L\"ange W ist \determ[Poisson-Verteilung]{Poissonverteilt}: +\item The average interval is $\mu_{ISI} = \frac{1}{\lambda}$ . +\item The variance of the intervals is $\sigma_{ISI}^2 = \frac{1}{\lambda^2}$ . +\item Thus, the coefficient of variation is always $CV_{ISI} = 1$ . +\item The serial correlation is $\rho_k =0$ for $k>0$, since the + occurrence of an event is independent of all previous events. Such a + process is also called a \enterm{renewal process}. +\item The number of events $k$ within a temporal window of duration + $W$ is Poisson distributed: \[ P(k) = \frac{(\lambda W)^ke^{\lambda W}}{k!} \] - (\figref{hompoissoncountfig}) -\item Der \determ{Fano Faktor} ist immer $F=1$ . +(\figref{hompoissoncountfig}) +\item The Fano Faktor is always $F=1$ . \end{itemize} \begin{exercise}{hompoissonspikes.m}{} - Schreibe eine Funktion \code{hompoissonspikes()}, die die Spikezeiten - eines homogenen Poisson-Prozesses mit gegebener Rate in Hertz f\"ur - eine Anzahl von trials gegebener maximaler L\"ange in Sekunden in - einem \codeterm{cell-array} zur\"uckgibt. Benutze die exponentiell-verteilten - Interspikeintervalle \eqnref{poissonintervals}, um die Spikezeiten zu erzeugen. + Implement a function \code{hompoissonspikes()} that uses a + homogeneous Poisson process to generate spike events at a given rate + for a certain duration and a number of trials. The rate should be + given in Hertz and the duration of the trials is given in + seconds. The function should return the event times in a + cell-array. Each entry in this array represents the events observed + in one trial. Apply \eqnref{poissonintervals} to generate the event + times. \end{exercise} \begin{figure}[t] \includegraphics[width=0.48\textwidth]{poissoncounthistdist100hz10ms}\hfill \includegraphics[width=0.48\textwidth]{poissoncounthistdist100hz100ms} - \titlecaption{\label{hompoissoncountfig}Z\"ahlstatistik von Poisson Spikes.}{} + \titlecaption{\label{hompoissoncountfig}Count statistics of Poisson + spiketrains.}{} \end{figure} %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -\section{Zeitabh\"angige Feuerraten} - -Bisher haben wir station\"are Spiketrains betrachtet, deren Statistik -sich innerhalb der Analysezeit nicht ver\"andert (station\"are -Punktprozesse). Meistens jedoch \"andert sich die Statistik der -Spiketrains eines Neurons mit der Zeit. Z.B. kann ein sensorisches -Neuron auf einen Reiz hin mit einer erh\"ohten Feuerrate antworten -(nichtstation\"arer Punktprozess). - -Wie die mittlere Anzahl der Spikes sich mit der Zeit ver\"andert, die -\determ{Feuerrate} $r(t)$, ist die wichtigste Gr\"o{\ss}e bei -nicht-station\"aren Spiketrains. Die Einheit der Feuerrate ist Hertz, -also Anzahl Aktionspotentiale pro Sekunde. Es gibt verschiedene -Methoden diese zu bestimmen. Drei solcher Methoden sind in Abbildung -\ref{psthfig} dargestellt. Alle Methoden haben ihre Berechtigung und -ihre Vor- und Nachteile. Im folgenden werden die drei Methoden aus -Abbildung \ref{psthfig} n\"aher erl\"autert. +\section{Time-dependent firing rate} + +So far we have discussed stationary spiketrains. The statistical properties +of these did not change within the observation time (stationary point +processes). Most commonly, however, this is not the case. A sensory +neuron, for example, might respond to a stimulus by modulating its +firing rate (non-stationary point process). + +How the firing rate $r(t)$ changes over time is the most important +measure, when analyzing non-stationary spike trains. The unit of the +firing rate is Hertz, i.e. the number of action potentials per +second. There are different ways to estimate the firing rate and three +of these methods are illustrated in \figref{psthfig}. All of +these have their own justifications and pros- and cons. In the +following we will discuss these methods more +closely. \begin{figure}[tp] \includegraphics[width=\columnwidth]{firingrates} - \titlecaption{Bestimmung der zeitabh\"angigen - Feuerrate.}{\textbf{A)} Rasterplot eines Spiketrains. \textbf{B)} - Feurerrate aus der instantanen Feuerrate bestimmt. \textbf{C)} - klassisches PSTH mit der Binning Methode. \textbf{D)} Feuerrate - durch Faltung mit einem Gauss Kern bestimmt.}\label{psthfig} + \titlecaption{Estimating the time-dependent firing + rate.}{\textbf{A)} Rasterplot depicting the spiking activity of a + neuron. \textbf{B)} Firing rate calculated from the + \emph{instantaneous rate}. \textbf{C)} classical PSTH with the + \emph{binning} method. \textbf{D)} Firing rate estimated by + \emph{convolution} of the activity with a Gaussian + kernel.}\label{psthfig} \end{figure} -\subsection{Instantane Feuerrate} +\subsection{Instantaneous firing rate} \begin{figure}[tp] \includegraphics[width=\columnwidth]{isimethod} - \titlecaption{Instantane Feuerrate.}{Skizze eines Spiketrains - (oben). Die Pfeile zwischen aufeinanderfolgenden - Aktionspotentialen mit den Zahlen in Millisekunden illustrieren - die Interspikeintervalle. Der Kehrwert des Interspikeintervalle - ergibt die instantane Feuerrate.}\label{instrate} + \titlecaption{Instantaneous firing rate.}{Sketch of the recorded + spiketrain (top). Arrows illustrate the interspike intervals and + number give the intervals in milliseconds. The inverse of the + interspike interval is the \emph{instantaneous firing + rate}.}\label{instratefig} \end{figure} -Ein sehr einfacher Weg, die zeitabh\"angige Feuerrate zu bestimmen ist -die sogenannte \determ[Feuerrate!instantane]{instantane Feuerrate} -(\enterm[firing rate!instantaneous]{instantaneous firing rate}). Dabei -wird die Feuerrate aus dem Kehrwert der Interspikeintervalle, der Zeit -zwischen zwei aufeinander folgenden Aktionspotentialen -(\figref{instrate} A), bestimmt. Die abgesch\"atzte Feuerrate -(\figref{instrate} B) ist g\"ultig f\"ur das gesammte -Interspikeintervall. Diese Methode hat den Vorteil, dass sie sehr -einfach zu berechnen ist und keine Annahme \"uber eine relevante -Zeitskala (der Kodierung oder des Auslesemechanismus der -postsynaptischen Zelle) macht. $r(t)$ ist allerdings keine -kontinuierliche Funktion, die Spr\"unge in der Feuerrate k\"onnen -f\"ur manche Analysen nachteilig sein. Au{\ss}erdem wird die Feuerrate -nie gleich Null, auch wenn lange keine Aktionspotentiale generiert -wurden. +A very simple method for estimating the time-dependent firing rate is +the \enterm[firing rate!instantaneous]{instantaneous firing rate}. The +firing rate can be directly estimated as the inverse of the time +between successive spikes, the interspike-interval +(\figref{instratefig}). + +\begin{equation} + \label{instantaneousrateeqn} + r_i = \frac{1}{T_i} . +\end{equation} + +The instantaneous rate $r_i$ is valid for the whole interspike +interval. The method has the advantage of being extremely easy to +compute and that it does not make any assumptions about the relevant +timescale (of the encoding in the neuron or the decoding of a +postsynaptic neuron). The resulting $r(t)$, however, is no continuous +function, the firing rate jumps from one level to the next. Since the +interspike interval between successive spikes is never infinitely +long, the firing rate never reaches zero despite that the neuron may +not fire an action potential for a long time. \begin{exercise}{instantaneousRate.m}{} - Implementiere die Absch\"atzung der Feuerrate auf Basis der - instantanen Feuerrate. Plotte die Feuerrate als Funktion der Zeit. + Implement a function that computes the instantaneous firing + rate. Plot the firing rate as a function of time. + \note{TODO: example data!!!} \end{exercise} -\subsection{Peri-Stimulus-Zeit-Histogramm} -W\"ahrend die Instantane Rate den Kehrwert der Zeit von einem bis zum -n\"achsten Aktionspotential misst, sch\"atzt das sogenannte -\determ{Peri-Stimulus-Zeit-Histogramm} (\enterm{peri stimulus time - histogram}, \determ[PSTH|see{Peri-Stimulus-Zeit-Histogramm}]{PSTH}) -die Wahrscheinlichkeit ab, zu einem Zeitpunkt Aktionspotentiale -anzutreffen. Es wird versucht die mittlere Rate \eqnref{firingrate} im -Grenzwert kleiner Beobachtungszeiten abzusch\"atzen: +\subsection{Peri-stimulus-time-histogram} +While the \emph{instantaneous firing rate} uses the interspike +interval, the \enterm{peri stimulus time histogram} (PSTH) uses the +spike count within observation windows of the duration $W$. It +estimates the probability of observing a spike within that observation +time. It tries to estimat the average rate in the limit of small +obersvation times \eqnref{psthrate}: \begin{equation} \label{psthrate} r(t) = \lim_{W \to 0} \frac{\langle n \rangle}{W} \; , \end{equation} -wobei die Anzahl $n$ der Aktionspotentiale, die im Zeitintervall -$(t,t+W)$ aufgetreten sind, \"uber trials gemittelt wird. Eine solche -Rate enspricht der zeitabh\"angigen Rate $\lambda(t)$ des inhomogenen -Poisson-Prozesses. +where $\langle n \rangle$ is the across trial average number of action +potentials observed within the interval $(t, t+W)$. Such description +matches the time-dependent firing rate $\lambda(t)$ of an +inhomogeneous Poisson process. -Das PSTH \eqnref{psthrate} kann entweder \"uber die Binning-Methode -oder durch Verfaltung mit einem Kern bestimmt werden. Beiden Methoden -gemeinsam ist die Notwendigkeit der Wahl einer zus\"atzlichen Zeitskala, -die der Beobachtungszeit $W$ in \eqnref{psthrate} entspricht. +The firing probability can be estimated using the \emph{binning + method} or by using \emph{kernel density estimations}. Both methods +make an assumption about the relevant observation time-scale ($W$ in +\eqnref{psthrate}). -\subsubsection{Binning-Methode} +\subsubsection{Binning-method} \begin{figure}[tp] \includegraphics[width=\columnwidth]{binmethod} - \titlecaption{Bestimmung des PSTH mit der Binning Methode.}{Der - gleiche Spiketrain wie in \figref{instrate}. Die grauen Linien - markieren die Grenzen der Bins und die Zahlen geben die Anzahl der Spikes - in jedem Bin an (oben). Die Feuerrate ergibt sich aus dem - mit der Binbreite normierten Zeithistogramm (unten).}\label{binpsth} + \titlecaption{Estimating the PSTH using the binning method.}{ The + same spiketrain as shown in \figref{instratefig} is used + here. Vertical gray lines indicate the borders between adjacent + bins in which the number of action potentials is counted (red + numbers). The firing rate is then the histogram normalized to the + binwidth.}\label{binpsthfig} \end{figure} -Bei der Binning-Methode wird die Zeitachse in gleichm\"aßige -Abschnitte (Bins) eingeteilt und die Anzahl Aktionspotentiale, die in -die jeweiligen Bins fallen, gez\"ahlt (\figref{binpsth} A). Um diese -Z\"ahlungen in die Feuerrate umzurechnen muss noch mit der Binweite -normiert werden. Das ist \"aquivalent zur Absch\"atzung einer -Wahrscheinlichkeitsdichte. Es kann auch die \code{hist()} Funktion zur -Bestimmung des PSTHs verwendet werden. \sindex[term]{Feuerrate!Binningmethode} - -Die bestimmte Feuerrate gilt f\"ur das gesamte Bin (\figref{binpsth} -B). Das so berechnete PSTH hat wiederum eine stufige Form, die von der -Wahl der Binweite anh\"angt. $r(t)$ ist also keine stetige -Funktion. Die Binweite bestimmt die zeitliche Aufl\"osung der -Absch\"atzung. \"Anderungen in der Feuerrate, die innerhalb eines Bins -vorkommen k\"onnen nicht aufgl\"ost werden. Mit der Wahl der Binweite -wird somit eine Annahme \"uber die relevante Zeitskala des Spiketrains -gemacht. +The binning method separates the time axis into regular bins of the +bin width $W$ and counts for each bin the number of observed action +potentials (\figref{binpsthfig} top). The resulting histogram is then +normalized with the bin width $W$ to yield the firing rate shown in +the bottom trace of figure \ref{binpsthfig}. The above sketched +process is equivalent to estimating the probability density. It is +possible to estimate the PSTH using the \code{hist()} method +\sindex[term]{Feuerrate!Binningmethode} + +The estimated firing rate is valid for the total duration of each +bin. This leads to the step-like plot shown in +\figref{binpsthfig}. $r(t)$ is thus not a contiunous function in +time. The binwidth defines the temporal resolution of the firing rate +estimation Changes that happen within a bin cannot be resolved. Thus +chosing a bin width implies an assumption about the relevant +time-scale. \pagebreak[4] \begin{exercise}{binnedRate.m}{} - Implementiere die Absch\"atzung der Feuerrate mit der ``binning'' - Methode. Plotte das PSTH. + Implement a function that estimates the firing rate using the + ``binning method''. The method should take the spike-times as an + input argument and returns the firing rate. Plot the PSTH. \end{exercise} -\subsubsection{Faltungsmethode} +\subsubsection{Convolution method --- Kernel density estimation} \begin{figure}[tp] \includegraphics[width=\columnwidth]{convmethod} - \titlecaption{Bestimmung des PSTH mit der Faltungsmethode.}{Der - gleiche Spiketrain wie in \figref{instrate}. Bei der Verfaltung - des Spiketrains mit einem Faltungskern wird jeder Spike durch den - Faltungskern ersetzt (oben). Bei korrekter Normierung des - Kerns ergibt sich die Feuerrate direkt aus der \"Uberlagerung der - Kerne.}\label{convrate} + \titlecaption{Estimating the firing rate using the convolution + method.}{The same spiketrain as in \figref{instratefig}. The + convolution of the spiketrain with a convolution kernel basically + replaces each spike event with the kernel (top). A Gaussian kernel + is used here, but other kernels are also possible. If the kernel + is properly normalized the firing rate results directly form the + superposition of the kernels.}\label{convratefig} \end{figure} -Bei der Faltungsmethode werden die harten Kanten der Bins der -Binning-Methode vermieden. Der Spiketrain wird mit einem Kern -verfaltet, d.h. jedes Aktionspotential wird durch den Kern ersetzt. -Zur Berechnung wird die Aktionspotentialfolge zun\"achst -``bin\"ar'' dargestellt. Dabei wird ein Spiketrain als -(Zeit-)Vektor dargestellt, in welchem die Zeitpunkte der -Aktionspotentiale als 1 notiert werden. Alle anderen Elemente des -Vektors sind 0. Anschlie{\ss}end wir dieser bin\"are Spiketrain mit -einem Gau{\ss}-Kern bestimmter Breite verfaltet: - \[r(t) = \int_{-\infty}^{\infty} \omega(\tau) \, \rho(t-\tau) \, {\rm d}\tau \; , \] -wobei $\omega(\tau)$ der Filterkern und $\rho(t)$ die bin\"are Antwort -ist. Bildlich geprochen wird jede 1 in $\rho(t)$ durch den Filterkern -ersetzt (Abbildung \ref{convrate} A). Wenn der Kern richtig normiert -wurde (Integral gleich Eins), ergibt sich die Feuerrate direkt aus der -\"Uberlagerung der Kerne (Abb. \ref{convrate} B). \sindex[term]{Feuerrate!Faltungsmethode} - -Die Faltungsmethode f\"uhrt, anders als die anderen Methoden, zu einer -stetigen Funktion was insbesondere f\"ur spektrale Analysen von -Vorteil sein kann. Die Wahl der Kernbreite bestimmt, \"ahnlich zur -Binweite, die zeitliche Aufl\"osung von $r(t)$. Die Breite des Kerns -macht also auch wieder eine Annahme \"uber die relevante Zeitskala des -Spiketrains. +With the convolution method we avoid the sharp edges of the binning +method. The spiketrain is convolved with a \enterm{convolution + kernel}. Technically speaking we need to first create a binary +representation of the spike train. This binary representation is a +series of zeros and ones in which the ones denote the spike. Then this binary vector is convolved with a kernel of a certain width: + +\[r(t) = \int_{-\infty}^{\infty} \omega(\tau) \, \rho(t-\tau) \, {\rm d}\tau \; , \] +where $\omega(\tau)$ represents the kernel and $\rho(t)$ the binary +representation of the response. The process of convolution can be +imagined as replacing each event of the spiketrain with the kernel +(figure \ref{convratefig} top). The superposition of the replaced +kernels is then the firing rate (if the kerel is correctly normalized +to an integral of one, figure \ref{convreffig} +bottom). \sindex[term]{Feuerrate!Faltungsmethode} + +In contrast to the other methods the convolution methods leads to a +continuous function which is often desirable (in particular when +applying methods in the frequency domain). The choice of the kernel +width defines, similar to the bin width of the binning method, the +temporal resolution of the method and thus makes assumptions about the +relevate time-scale. \pagebreak[4] \begin{exercise}{convolutionRate.m}{} - Verwende die Faltungsmethode um die Feuerrate zu bestimmen. Plotte - das Ergebnis. + Implement the function that estimates the firing rate using the + convolution method. The method takes the spiketrain, the temporal + resolution of the recording (as the stepsize $dt$, in seconds) and + the width of the kernel (the standard deviation $\sigma$ of the + Gaussian kernel, in seconds) as input arguments. It returns the + firing rate. Plot the result. \end{exercise} \section{Spike-triggered Average} -Die graphischer Darstellung der Feuerrate allein reicht nicht aus um -den Zusammenhang zwischen neuronaler Antwort und einem Stimulus zu -analysieren. Eine Methode um mehr \"uber diesen Zusammenhang zu -erfahren, ist der \enterm{spike-triggered average} -(\enterm[STA|see{spike-triggered average}]{STA}). Der STA + +The graphical representation of the neuronal activity alone is not +sufficient tot investigate the relation between the neuronal response +and a stimulus. One method to do this is the \enterm[STA|see + spike-triggered average]{spike-triggered average}. The STA \begin{equation} STA(\tau) = \langle s(t - \tau) \rangle = \frac{1}{N} \sum_{i=1}^{N} s(t_i - \tau) \end{equation} -der $N$ Aktionspotentiale zu den Zeiten $t_i$ in Anwort auf den -Stimulus $s(t)$ ist der mittlere Stimulus, der zu einem -Aktionspotential in der neuronalen Antwort f\"uhrt. - -Der STA l\"a{\ss}t sich relativ einfach berechnen, indem aus dem -Stimulus f\"ur jeden beobachteten Spike ein entsprechender Abschnitt -ausgeschnitten wird und diese dann gemittelt werde (\figref{stafig}). +of $N$ action potentials observed at the times $t_i$ in response to +the stimulus $s(t)$ is the average stimulus that led to a spike in the +neuron. The STA can be easily extracted by cutting snippets out of +$s(t)$ that surround the times of the spikes. The resulting stimulus +snippets are then averaged (\figref{stafig}). \begin{figure}[t] \includegraphics[width=\columnwidth]{sta} - \titlecaption{Spike-triggered Average eines P-Typ Elektrorezeptors - und Stimulusrekonstruktion.}{Der STA (links): der Rezeptor - wurde mit einem ``white-noise'' Stimulus getrieben. Zeitpunkt 0 - ist der Zeitpunkt des beobachteten Aktionspotentials. Die Kurve - ergibt sich aus dem gemittelten Stimulus je 50\,ms vor und nach - einem Aktionspotential. Stimulusrekonstruktion mittels - STA (rechts). Die Zellantwort wird mit dem STA gefaltet um eine - Rekonstruktion des Stimulus zu erhalten.}\label{stafig} + \titlecaption{Spike-triggered average of a P-type electroreceptors + and the stimulus reconstruction.}{The neuron was driven by a + ``white-noise'' stimulus (blue curve, right panel). The STA (left) + is the average stimulus that surrounds the times of the recorded + action potentials (40\,ms before and 20\,ms after the + spike). Using the STA as a convolution kernel for convolving the + spiketrain we can reconstruct the stimulus from the neuronal + response. In this way we can get an impression of the stimulus + features that are encoded in the neuronal response (right + panel).}\label{stafig} \end{figure} -Aus dem STA k\"onnen verschiedene Informationen \"uber den -Zusammenhang zwischen Stimulus und neuronaler Antwort gewonnen -werden. Die Breite des STA repr\"asentiert die zeitliche Pr\"azision, -mit der Stimulus und Antwort zusammenh\"angen und wie weit das Neuron -zeitlich integriert. Die Amplitude des STA (gegeben in der gleichen -Einheit wie der Stimulus) deutet auf die Empfindlichkeit des Neurons -bez\"uglich des Stimulus hin. Eine hohe Amplitude besagt, dass es -einen starken Stimulus ben\"otigt, um ein Aktionspotential -hervorzurufen. Aus dem zeitlichen Versatz des STA kann die Zeit -abgelesen werden, die das System braucht, um auf den Stimulus zu -antworten. - -Der STA kann auch dazu benutzt werden, aus den Antworten der Zelle den -Stimulus zu rekonstruieren (\figref{stafig} B). Bei der -\determ[invertierte Rekonstruktion]{invertierten Rekonstruktion} wird -die Zellantwort mit dem STA verfaltet. +From the STA we can extract several pieces of information about the +relation of stimulus and response. The width of the STA represents the +temporal precision with which the neuron encodes the stimulus waveform +and in which temporal window the neuron integrates the (sensory) +input. The amplitude of the STA tells something about the sensitivity +of the neuron. The STA is given in the same units as the stimulus and +a small amplitude indicates that the neuron needs only a small +stimulus amplitude to create a spike, a large amplitude on the +contrary suggests the opposite. The temporal delay between the STA and +the time of the spike is a consequence of the time the system (neuron) +needs to process the stimulus. + +We can further use the STA to do a \enterm{reverse reconstruction} and +estimate the stimulus from the neuronal response (\figref{stafig}, +right). For this, the spiketrain is convolved with the STA as a +kernel. \begin{exercise}{spikeTriggeredAverage.m}{} - Implementiere eine Funktion, die den STA ermittelt. Verwende dazu - den Datensatz \file{sta\_data.mat}. Die Funktion sollte folgende - R\"uckgabewerte haben: + Implement a function that calculates the STA. Use the dataset + \file{sta\_data.mat}. The function expects the spike train, the + stimulus and the temporal resolution of the recording as input + arguments and should return the following information: \vspace{-1ex} \begin{itemize} \setlength{\itemsep}{0ex} - \item den Spike-Triggered-Average. - \item die Standardabweichung der individuellen STAs. - \item die Anzahl Aktionspotentiale, die zur Berechnung des STA verwendet wurden. - \vspace{-2ex} + \item the spike-triggered average. + \item the standard deviation of the STA across the individual snippets. + \item The number of action potentials used to estimate the STA. \end{itemize} + \end{exercise} \begin{exercise}{reconstructStimulus.m}{} - Rekonstruiere den Stimulus mithilfe des STA und der Spike - Zeiten. Die Funktion soll Vektor als R\"uckgabewert haben, der - genauso gro{\ss} ist wie der Originalstimulus aus der Datei + Do the reverse reconstruction using the STA and the spike times. The + function should return the estimated stimulus in a vector that has + the same size as the original stimulus contained in file \file{sta\_data.mat}. \end{exercise} diff --git a/pointprocesses/lecture/pointprocesses_de.tex b/pointprocesses/lecture/pointprocesses_de.tex new file mode 100644 index 0000000..4b0850f --- /dev/null +++ b/pointprocesses/lecture/pointprocesses_de.tex @@ -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} diff --git a/pointprocesses/lecture/rasterexamples.py b/pointprocesses/lecture/rasterexamples.py index 9c66897..042bda2 100644 --- a/pointprocesses/lecture/rasterexamples.py +++ b/pointprocesses/lecture/rasterexamples.py @@ -31,7 +31,7 @@ def pifspikes(input, trials, dt, D=0.1) : times = [] v = vreset noise = np.sqrt(2.0*D)*np.random.randn(len(input))/np.sqrt(dt) - for k in xrange(len(noise)) : + for k in range(len(noise)) : v += (input[k]+noise[k])*dt/tau if v >= vthresh : v = vreset @@ -55,7 +55,7 @@ rng = np.random.RandomState(54637281) time = np.arange(0.0, duration, dt) x = np.zeros(time.shape)+rate n = rng.randn(len(time))*drate*tau/np.sqrt(dt)+rate -for k in xrange(1,len(x)) : +for k in range(1,len(x)) : x[k] = x[k-1] + (n[k]-x[k-1])*dt/tau x[x<0.0] = 0.0 diff --git a/pointprocesses/lecture/returnmapexamples.py b/pointprocesses/lecture/returnmapexamples.py index 96803ff..829cf6d 100644 --- a/pointprocesses/lecture/returnmapexamples.py +++ b/pointprocesses/lecture/returnmapexamples.py @@ -31,7 +31,7 @@ def pifspikes(input, trials, dt, D=0.1) : times = [] v = vreset noise = np.sqrt(2.0*D)*np.random.randn(len(input))/np.sqrt(dt) - for k in xrange(len(noise)) : + for k in range(len(noise)) : v += (input[k]+noise[k])*dt/tau if v >= vthresh : v = vreset @@ -41,7 +41,7 @@ def pifspikes(input, trials, dt, D=0.1) : def isis( spikes ) : isi = [] - for k in xrange(len(spikes)) : + for k in range(len(spikes)) : isi.extend(np.diff(spikes[k])) return np.array( isi ) @@ -84,7 +84,7 @@ rng = np.random.RandomState(54637281) time = np.arange(0.0, duration, dt) x = np.zeros(time.shape)+rate n = rng.randn(len(time))*drate*tau/np.sqrt(dt)+rate -for k in xrange(1,len(x)) : +for k in range(1,len(x)) : x[k] = x[k-1] + (n[k]-x[k-1])*dt/tau x[x<0.0] = 0.0 diff --git a/pointprocesses/lecture/serialcorrexamples.py b/pointprocesses/lecture/serialcorrexamples.py index e7fff9c..a2dc853 100644 --- a/pointprocesses/lecture/serialcorrexamples.py +++ b/pointprocesses/lecture/serialcorrexamples.py @@ -31,7 +31,7 @@ def pifspikes(input, trials, dt, D=0.1) : times = [] v = vreset noise = np.sqrt(2.0*D)*np.random.randn(len(input))/np.sqrt(dt) - for k in xrange(len(noise)) : + for k in range(len(noise)) : v += (input[k]+noise[k])*dt/tau if v >= vthresh : v = vreset @@ -41,7 +41,7 @@ def pifspikes(input, trials, dt, D=0.1) : def isis( spikes ) : isi = [] - for k in xrange(len(spikes)) : + for k in range(len(spikes)) : isi.extend(np.diff(spikes[k])) return np.array( isi ) @@ -95,7 +95,7 @@ rng = np.random.RandomState(54637281) time = np.arange(0.0, duration, dt) x = np.zeros(time.shape)+rate n = rng.randn(len(time))*drate*tau/np.sqrt(dt)+rate -for k in xrange(1,len(x)) : +for k in range(1,len(x)) : x[k] = x[k-1] + (n[k]-x[k-1])*dt/tau x[x<0.0] = 0.0 diff --git a/pointprocesses/lecture/sta.py b/pointprocesses/lecture/sta.py index c5d702b..127c5ae 100644 --- a/pointprocesses/lecture/sta.py +++ b/pointprocesses/lecture/sta.py @@ -6,7 +6,7 @@ from IPython import embed def plot_sta(times, stim, dt, t_min=-0.1, t_max=.1): count = 0 - sta = np.zeros((abs(t_min) + abs(t_max))/dt) + sta = np.zeros(int((abs(t_min) + abs(t_max))/dt)) time = np.arange(t_min, t_max, dt) if len(stim.shape) > 1 and stim.shape[1] > 1: stim = stim[:,1] diff --git a/programming/exercises/boolean_logical_indexing.tex b/programming/exercises/boolean_logical_indexing.tex index d91b3d8..315efe7 100644 --- a/programming/exercises/boolean_logical_indexing.tex +++ b/programming/exercises/boolean_logical_indexing.tex @@ -62,9 +62,11 @@ following pattern: ``variables\_datatypes\_\{lastname\}.m'' \part Execute and explain: \verb+bitand(10, 8)+ \part Execute and explain: \verb+bitor(10, 8)+ \end{parts} -\item Implement the following Boolean expressions. Test using randomly selected integer values for \verb+x+ and \verb+y+. +\item Implement the following Boolean expressions. Test your + implementations using random integer + numbers for \verb+x+ and \verb+y+ (\verb+randi+). \begin{parts} - \part The result should be \verb+true+ if \verb+x+ greater than \verb+y+ and the sum of \verb+x+ and \verb+y+ is not less than 100. + \part The result should be \verb+true+ if \verb+x+ is greater than \verb+y+ and the sum of \verb+x+ and \verb+y+ is not less than 100. \part The result should be \verb+true+ if \verb+x+ and \verb+y+ are not equal zero or \verb+x+ and \verb+y+ are equal. \end{parts} \end{questions} @@ -82,7 +84,7 @@ matrices that match in certain criteria. This process is called the following commands and explain. \begin{parts} \part \verb+x < 5+ - \part \verb+x( x < 5) )+ + \part \verb+x( (x < 5) )+ \part \verb+x( (y <= 2) )+ \part \verb+x( (x > 2) | (y < 8) )+ \part \verb+x( (x == 0) & (y == 0) )+ @@ -94,6 +96,19 @@ matrices that match in certain criteria. This process is called \verb+x >= 33 and x < 66+ to 1 and all \verb+x >= 66+ to 2. \part Count the number of elements in each class using Boolean expressions (\verb+sum+ can be used to count the matches). \end{parts} + + + \question Plotting a periodic signal. + \begin{parts} + \part Load the file ``signal.mat'' into the workspace (use + \verb+load('signal.mat')+ or use the UI for it). It contains two + variables \verb+signal+ and \verb+time+. The signal has a period of + 0.2\,s + \part What is the size of the data? + \part What is the temporal resolution of the time axis? + \part Plot the full data. Make sure that the axes are properly labeled (script chapter 3, or the matlab documentation). + \part Use logical indexing to select the periods individually and plot them into the same plot (first period starts at time 0.0\,s, the next at 0.2\,s and so on). + \end{parts} \end{questions} \end{document} diff --git a/programming/exercises/control_flow.tex b/programming/exercises/control_flow.tex index 725a81e..84d48da 100644 --- a/programming/exercises/control_flow.tex +++ b/programming/exercises/control_flow.tex @@ -14,7 +14,7 @@ %%%%% text size %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \usepackage[left=20mm,right=20mm,top=25mm,bottom=25mm]{geometry} \pagestyle{headandfoot} -\header{{\bfseries\large Exercise 4}}{{\bfseries\large Control Flow}}{{\bfseries\large 30. October, 2017}} +\header{{\bfseries\large Exercise 5}}{{\bfseries\large Control Flow}}{{\bfseries\large 30. October, 2018}} \firstpagefooter{Dr. Jan Grewe}{Phone: 29 74588}{Email: jan.grewe@uni-tuebingen.de} \runningfooter{}{\thepage}{} @@ -33,39 +33,44 @@ \begin{center} \textbf{\Large Introduction to scientific computing}\\[1ex] {\large Jan Grewe, Jan Benda}\\[-3ex] - Abteilung Neuroethologie \hfill --- \hfill Institut f\"ur Neurobiologie \hfill --- \hfill \includegraphics[width=0.28\textwidth]{UT_WBMW_Black_RGB} \\ + Neuroethologie \hfill --- \hfill Institute for Neurobiology \hfill --- \hfill \includegraphics[width=0.28\textwidth]{UT_WBMW_Black_RGB} \\ \end{center} -The exercises are meant for self-monitoring, revision of the lecture -topic. You should try to solve them on your own. Your solution should -be submitted as a single script (m-file) in the Ilias system. Each -task should be solved in its own ``cell''. Each cell must be -executable on its own. The file should be named according to the following pattern: -``variables\_datatypes\_\{lastname\}.m'' benannt werden +The exercises are meant for self-monitoring and revision of the +lecture. You should try to solve them on your own. Your solution +should be submitted as a single script (m-file) in the Ilias +system. Each task should be solved in its own ``cell''. Each cell must +be executable on its own. The file should be named according to the +following pattern: ``variables\_datatypes\_\{lastname\}.m'' (e.g. variables\_datentypes\_mueller.m). \begin{questions} \question Implement \code{for} loops in which the \emph{running variable}: \begin{parts} \part ... assumes values from 0 to 10. Print the value of the running variable for each iteration step. + \begin{solution} + for i = 1:10 + disp(i); + end; + \end{solution} \part ... assumes values from 10 to 0. Print out the value for each iteration. \part ... assumes values from 0 to 1 in steps of 0.1. Print out the value for each iteration. \end{parts} - \question Indexing in vectros: + \question Indexing in vectors: \begin{parts} - \part Define a vector \code{x} that contains values from 1:100. + \part Define a vector \code{x} that contains values ranging from 101 to 200. \part Use a \code{for} loop to print each element of \code{x}. Use the running variable for the indexing. \begin{solution} \code{for i = 1:length(x); disp(x(i)); end} \end{solution} - \part Do the same without using the running variable for indexing. + \part Do the same without use the running variable directly (not for indexing). \begin{solution} \code{for i : x; disp(i); end;} \end{solution} \end{parts} - \question Create a vector \verb+x+ that contains 50 random numbers in the range 0 - 10. + \question Create a vector \verb+x+ that contains 50 random numbers in the value range 0 - 10. \begin{parts} \part Use a loop to calculate the arithmetic mean. The mean is defined as: $\overline{x}=\frac{1}{n}\sum\limits_{i=0}^{n}x_i $. @@ -74,20 +79,20 @@ executable on its own. The file should be named according to the following patte \part Search the MATLAB help for functions that do the calculations for you :-). \end{parts} - \question Implement a\code{while} loop + \question Implement a \code{while} loop \begin{parts} - \part ... that iterates 100 times. Print out the current iteration number. + \part ... that iterates 100 times. Print out the current iteration number. \part ... iterates endlessly. You can interrupt the execution with the command \code{Strg + C}. \end{parts} - \question Use an endless \code{while} loop to individually print out the elements of a vector that contains 10 elements. + \question Use an endless \code{while} loop to individually print out the elements of a vector that contains 10 elements. Stop the execution after all elements have been printed. \question Use the endless \verb+while+ loop to draw random numbers - (\verb+randn+) until you it is larger than 1.33. + (\verb+randn+) until you hit one that is larger than 1.33. \begin{parts} \part Count the number of required iterations. - \part Use a \code{for} loop to run the previous test 1000 times. Remeber all counts and calculate the mean of it. - \part Create a plot that schows for each of the 1000 trials, how often you had to draw. + \part Use a \code{for} loop to run the previous test 1000 times. Remember all counts and calculate the mean of it. + \part Create a plot that shows for each of the 1000 trials, how often you had to draw. \part Play around with the threshold. What happens? \end{parts} @@ -100,19 +105,19 @@ executable on its own. The file should be named according to the following patte \end{parts} \question Test the random number generator! To do so count the - number of elements that fall in the classes defined by the edges + number of elements that fall into the classes defined by the edges [0.0, 0.2, 0.4, 0.6, 0.8, 1.0]. Store the results in a vector. Use a loop to draw 1000 random numbers with \verb+rand()+ (see help). What would be the expectation? What is the result? - \question String parsing: It is quite common to use the names of datasets to encode the conditions under which they were recorded. To find out the required information we have to \emph{parse} he name. + \question String parsing: It is quite common to use the names of datasets to encode the conditions under which they were recorded. To find out the required information we have to \emph{parse} the name. \begin{parts} \part Create a variable \verb+filename = '2015-10-12_100Hz_1.25V.dat'+. Obviously, the underscore was used as a delimiter. \part Use a \verb+for+ loop to find the positions of the underscores and store these in a vector. \part Use a second loop to iterate over the previously filled `position' vector and use this information to - cut filename appropriately. + cut \verb+filename+ appropriately. \part Print out the individual parts. \end{parts} diff --git a/programming/exercises/scripts_functions.tex b/programming/exercises/scripts_functions.tex index 652a88c..2b3282a 100644 --- a/programming/exercises/scripts_functions.tex +++ b/programming/exercises/scripts_functions.tex @@ -16,7 +16,7 @@ %%%%% text size %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \usepackage[left=20mm,right=20mm,top=25mm,bottom=25mm]{geometry} \pagestyle{headandfoot} -\header{{\bfseries\large \"Ubung 5}}{{\bfseries\large Scripts and functions}}{{\bfseries\large 07. November, 2017}} +\header{{\bfseries\large Exercise 6}}{{\bfseries\large Scripts and functions}}{{\bfseries\large 06. November, 2018}} \firstpagefooter{Prof. Jan Benda}{Phone: 29 74 573}{Email: jan.benda@uni-tuebingen.de} \runningfooter{}{\thepage}{} @@ -46,7 +46,6 @@ aboveskip=10pt } - \newcommand{\code}[1]{\texttt{#1}} %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% @@ -60,7 +59,7 @@ \end{center} The exercises are meant for self-monitoring and revision of the -lecture topic. You should try to solve them on your own. In contrast +lecture. You should try to solve them on your own. In contrast to previous exercises, the solutions can not be saved in a single file but each question needs an individual file. Combine the files into a single zip archive and submit it via ILIAS. Name the archive according @@ -108,7 +107,7 @@ to the pattern: ``scripts\_functions\_\{surname\}.zip''. \part{} Improve the function that it takes the duration of the time axis, - the amplitude and the frequency as input arguments. The + the amplitude, and the frequency as input arguments. The calculation should use a temporal stepsize that is 0.01 of the frequency. \begin{solution} @@ -144,7 +143,8 @@ to the pattern: ``scripts\_functions\_\{surname\}.zip''. \question Random Walk. In a 1-D random walk an \emph{agent} walks randomly either in the one ($+1$) or the other ($-1$) direction. With each simulation step one direction is chosen and the - position is updated accordingly. + position is updated accordingly. \textbf{There are some differences + to the solution discussed in the lecture!} \begin{parts} \part{} @@ -172,7 +172,7 @@ to the pattern: ``scripts\_functions\_\{surname\}.zip''. Now we want to know how the probability of $p_{+1}$ (the probability to walk into the $+1$ direction) impacts the random walk. Vary $p_{+1}$ in the range $0.5 \le p_{+1} < 0.8$. Do 10 - random walks for four probabilities (apply the same thresholds for + random walks for the four probabilities (apply the same thresholds for stopping the simulations as before). \begin{solution} \lstinputlisting{randomwalkscriptc.m} diff --git a/programming/exercises/variables_types.tex b/programming/exercises/variables_types.tex index b2b7dca..fb6c092 100644 --- a/programming/exercises/variables_types.tex +++ b/programming/exercises/variables_types.tex @@ -1,4 +1,4 @@ -\documentclass[12pt,a4paper,pdftex]{exam} +\documentclass[12pt,a4paper,pdftex, answers]{exam} \usepackage[german]{babel} \usepackage{natbib} diff --git a/programming/lecture/programming.tex b/programming/lecture/programming.tex index 65a907d..91c816a 100644 --- a/programming/lecture/programming.tex +++ b/programming/lecture/programming.tex @@ -40,14 +40,12 @@ variable. \begin{figure} \centering \begin{subfigure}{.5\textwidth} - \includegraphics[width=0.8\textwidth]{variable} - \label{variable:a} + \includegraphics[width=0.8\textwidth]{variable}\label{variable:a} \end{subfigure}% \begin{subfigure}{.5\textwidth} - \includegraphics[width=.8\textwidth]{variableB} - \label{variable:b} + \includegraphics[width=.8\textwidth]{variableB}\label{variable:b} \end{subfigure} - \titlecaption{Variables} point to a memory + \titlecaption{Variables}{ point to a memory address. They further are described by their name and data type. The variable's value is stored as a pattern of binary values (0 or 1). When reading the variable this pattern is @@ -630,7 +628,7 @@ matrix). The function \code{cat()} allows to concatenate n-dimensional matrices. To request the length of a vector we used the function -\code{length()}. This function is \tetbf{not} suited to request +\code{length()}. This function is \textbf{not} suited to request information about the size of a matrix. As mentioned above, \code{length()} would return the length of the largest dimension. The function \code{size()} however, returns the length in each dimension @@ -1031,7 +1029,7 @@ segment of data of a certain time span (the stimulus was on, \varcode{x} represent the measured data at the times in \varcode{t}. \item Use logical indexing to select those values that have been - recorded in the time span form 5--6\,s. + recorded in the time span from 5--6\,s. \end{itemize} \end{exercise} @@ -1095,7 +1093,7 @@ ans = Generally, a program is executed line by line from top to bottom. Sometimes this is not the desired behavior, or the other way round, it is needed to skip certain parts or execute others -repeatedly. High-level programming languages like \matlab{} offer +repeatedly. High-level programming languages such as \matlab{} offer statements that allow to manipulate the control flow. There are two major classes of such statements: @@ -1120,7 +1118,7 @@ x = 120 \end{lstlisting} -This kind of program is fine but it is rather repetitive. The only +This kind of program solves the taks but it is rather repetitive. The only thing that changes is the increasing factor. The repetition of such very similar lines of code is bad programming style. This is not only a matter of taste but there are severe drawbacks to this style: @@ -1136,7 +1134,7 @@ a matter of taste but there are severe drawbacks to this style: repetitions. It is easy to forget a single change. \item Readability: repetitive code is terrible to read and to understand. (I) one tends to skip repetitions (its the same, - anyways) and misses the essential change. (II), the duplication of + anyway) and misses the essential change. (II), the duplication of code leads to long and hard-to-parse programs. \end{enumerate} All imperative programming languages offer a solution: the loop. It is @@ -1146,7 +1144,7 @@ used whenever the same commands have to be repeated. \subsubsection{The \code{for} --- loop} The most common type of loop is the \codeterm{for-loop}. It consists of a \codeterm[Loop!head]{head} and the -\codeterm[Loop!body]{body}. The head defines how often the code of the +\codeterm[Loop!body]{body}. The head defines how often the code in the body is executed. In \matlab{} the head begins with the keyword \code{for} which is followed by the \codeterm{running variable}. In \matlab{} a for-loop always operates on vectors. With each @@ -1191,7 +1189,7 @@ while x == true % head with a Boolean expression % execute this code if the expression yields true end \end{lstlisting} - + \begin{exercise}{factorialWhileLoop.m}{} Implement the factorial of a number \varcode{n} using a \code{while} -- loop. @@ -1201,7 +1199,7 @@ end \begin{exercise}{neverendingWhile.m}{} Implement a \code{while}--loop that is never-ending. Hint: the body is executed as long as the Boolean expression in the head is - true. You can escape the loop by pressing \keycode{Ctrl+C}. + \code{true}. You can escape the loop by pressing \keycode{Ctrl+C}. \end{exercise} diff --git a/programmingstyle/code/calculateSines.m b/programmingstyle/code/calculateSines.m index 8a2139a..b25388e 100644 --- a/programmingstyle/code/calculateSines.m +++ b/programmingstyle/code/calculateSines.m @@ -1,5 +1,7 @@ function sines = calculateSines(x, amplitudes, frequencies) - % Function calculates sinewaves with all combinations of + % sines = calculateSines(x, amplitudes, frequencies) + % + % Function calculates sinewaves with all combinations of % given amplitudes and frequencies. % Arguments: x, a vector of radiants for which the sine should be % computed. @@ -8,11 +10,11 @@ function sines = calculateSines(x, amplitudes, frequencies) % % Returns: a 3-D Matrix of sinewaves, 2nd dimension represents % the amplitudes, 3rd the frequencies. - + sines = zeros(length(x), length(amplitudes), length(frequencies)); for i = 1:length(amplitudes) - sines(:,i,:) = sinesWithFrequencies(x, amplitudes(i), frequencies); + sines(:,i,:) = sinesWithFrequencies(x, amplitudes(i), frequencies); end end diff --git a/programmingstyle/lecture/programmingstyle.tex b/programmingstyle/lecture/programmingstyle.tex index 9a048e7..a4e29ac 100644 --- a/programmingstyle/lecture/programmingstyle.tex +++ b/programmingstyle/lecture/programmingstyle.tex @@ -6,20 +6,19 @@ %\selectlanguage{ngerman} Cultivating a good code style not a matter of good taste but is -a key ingredient for understandability, maintainability and in the end +a key ingredient for understandability, maintainability and, in the end, facilitates reproducibility of scientific results. Programs should be written and structured in a way that supports outsiders as well the author himself --- a few weeks or months after it was written --- to understand the programs' rationale. Clean code pays off for the original author as well as others that are supposed to use the code. - Clean code addresses several issues: \begin{enumerate} \item The programs' structure. \item Naming of scripts and functions. \item Naming of variables and constants. -\item Application of indentation empty lines to define blocks. +\item Application of indentation and empty lines to define blocks. \item Use of comments and inline documentation. \item Delegation of repeated code to functions and dedicated subroutines. @@ -29,13 +28,12 @@ Clean code addresses several issues: While introducing scripts and functions we suggested a typical program layout (box\,\ref{whenscriptsbox}). The idea is to create a single -entry point by having one script that controls the rest of program by -managing data and results and calling functions that work on the data -and produce the results. Applying this structure makes it easy to -understand the flow of the program but two questions remain: (i) How -to organize the files on the file system and (ii) how to name them -that the controlling script is easily identified among the other -\codeterm{m-files}. +entry point by having one script that controls the rest of the program +by calling functions that work on the data and managing the +results. Applying this structure makes it easy to understand the flow +of the program but two questions remain: (i) How to organize the files +on the file system and (ii) how to name them that the controlling +script is easily identified among the other \codeterm{m-files}. Upon installation ``MATLAB'' creates a folder called \emph{MATLAB} in the user space (Windows: My files, Linux: Documents, MacOS: @@ -110,24 +108,11 @@ Box~\ref{matlabpathbox}). \end{lstlisting} \end{ibox} -\section{Naming scripts and functions} -\matlab{} will search the search path (Box \ref{matlabpathbox}) -exclusively by name. It is case-sensitive this implies that the files -\file{test\_function.m} and \file{Test\_function.m} are two different -things. It is self-evident that choosing such names is nonsensical -because the name contains no cue about the difference between the two -and it further tells close to nothing about the purpose. Finding good -names is not trivial sometimes it is harder than the programming -itself. Expressive names, however, pay off! Expressive means that the -name provides information about the purpose. - -\begin{important}[Naming scripts and functions] - Function and script names should be expressive in the sense that the - name provides information about the function's purpose - (\file{estimate\_firingrate.m} tells much more than - \file{exercise1.m}). Choosing a good name replaces large parts of - the documentation. -\end{important} +\section{Naming things} +The dictum of good code style is: ``Program code must be readable.'' +Expressive names are extraordinarily important in this respect. Even +if it is tricky to find expressive names that are not overly long, +naming should be taken seriously. \matlab{} has a few rules about names: Names must not start with a number, they must not contain blanks or other special characters like @@ -144,26 +129,41 @@ patterns: There are other common patterns such as the \emph{camelCase} in which the first character of compound words is capitalized. Other -conventions use the underscore to separate the individual words -\emph{snake\_case}. A function that counts the number of action +conventions use the underscore to separate the individual words ( +\emph{snake\_case}). A function that counts the number of action potentials could be named \file{spikeCount.m} or \file{spike\_count.m}. +The same naming rules apply for scripts and functions as well as +variables and constants. -\section{Naming variables and constants} +\subsection{Naming scripts and functions} +\matlab{} will search the search path (Box \ref{matlabpathbox}) +exclusively by name. This search is case-sensitive which implies that +the files \file{test\_function.m} and \file{Test\_function.m} are two +different things. It is self-evident that choosing such names is +nonsensical because the tiny difference in the name contains no cue +about the difference between the two versions and the function names +themselves tell close to nothing about the purpose. Finding good names +is not trivial. Sometimes it is harder than the programming +itself. Choosing \emph{expressive names} that provide information about a +function's purpose, however, pays off! + +\begin{important}[Naming scripts and functions] + Names of functions and scripts should be expressive in the sense + that the name provides information about the function's purpose. + (\file{estimate\_firingrate.m} tells much more than + \file{exercise1.m}). Choosing a good name replaces large parts of + the documentation. +\end{important} -\matlab{} applies the same rules for naming variables and constants as -for the naming of scripts and functions. The dictum of good -code style is: ``Program code must be readable.'' Expressive -names are extraordinarily important in this respect. Even if it is -tricky to find expressive names that are not overly long, naming -should be taken seriously. +\subsection{Naming variables and constants} While the names of scripts and functions describe the purpose, names of variables describe the stored content. A variable storing the -average number of actions potentials could be called +average number of actions potentials could be called\\ \varcode{average\_spike\_count}. If this variable is meant to store -multiple spike counts the plural form would be appropriate +multiple spike counts the plural form would be appropriate\\ (\varcode{average\_spike\_counts}). The control variables used in the head of a \code{for} loop are often @@ -190,7 +190,7 @@ to comprehend. Even though the \matlab{} language (as many others) does not enforce indentation, indentation is very powerful for defining coherent blocks. The \matlab{} editor supports this by an auto-indentation mechanism. A selected section of the code and be -automatically indented by pressing the \keycode{Ctrl-I} combination. +automatically indented by pressing \keycode{Ctrl-I}. Interspersing empty lines is very helpful to separate regions in the code that belong together. Too many empty lines, however lead to @@ -198,8 +198,11 @@ hard-to-read code because it might require more space than a granted by the screen and thus takes overview. The following two listings show basically the same implementation of a -random walk once in a rather chaotic version (listing -\ref{chaoticcode}) then in cleaner way (listing \ref{cleancode}) +random walk\footnote{A random walk is a simple simulation of Brownian + motion. In each simulation step an agent takes a step into a + randomly chosen direction.} once in a rather chaotic version +(listing \ref{chaoticcode}) then in cleaner way (listing +\ref{cleancode}) \begin{lstlisting}[label=chaoticcode, caption={Chaotic implementation of the random-walk.}] num_runs = 10; max_steps = 1000; @@ -245,25 +248,24 @@ end \section{Using comments} It is common to provide extra information about the meaning of program -code by adding comments to it. In \matlab{} comments are indicated by -the percent character \code{\%}. Anything that is written in the -respective line following the percent is ignored and considered a -comment. When used sparsely comments can immensely important for -understanding. Comments are short sentences that describe the meaning -of the (following) lines in the program code. During the initial -implementation of a function they can be used to guide the development -but have the tendency to blow up the code and decrease readability. By -choosing expressive variable and function names, most lines should be -self-explanatory. - -For example stating the obvious does not really help:\\ -\varcode{ x = x + 2; \% add two to x}\\ +code by adding comments. In \matlab{} comments are indicated by the +percent character \code{\%}. Anything that follows the percent +character in a line is ignored and considered a comment. When used +sparsely comments can be immensely helpful. Comments +are short sentences that describe the meaning of the (following) lines +in the program code. During the initial implementation of a function +they can be used to guide the development but have the tendency to +blow up the code and decrease readability. By choosing expressive +variable and function names, most lines should be self-explanatory. + +For example stating the obvious does not really help and should be +avoided:\\ \varcode{ x = x + 2; \% add two to x}\\ \begin{important}[Using comments] \begin{itemize} \item Comments describe the rationale of the respective code block. - \item Comments are good and helpful --- they have to be true, however! - \item A wrong comment is worse than a non-existent comment! + \item Comments are good and helpful --- they must be true, however! + \item A wrong comment is worse than a non-existent one! \item Comments must be maintained just as the code. Otherwise they may become wrong and worse than meaningless! \end{itemize} @@ -303,7 +305,7 @@ well documented function. Comments and empty lines are used to organize code into logical blocks and to briefly explain what they do. Whenever one feels tempted to do this, one could also consider to delegate the respective task to a -function. In most cases this is preferable. +function. In most cases this is preferable. Not delegating the tasks leads to very long \codeterm{m-files} which can be confusing. Sometimes such a code is called ``spaghetti @@ -325,9 +327,9 @@ Generally, functions live in their own \codeterm{m-files} that have the same name as the function itself. Delegating tasks to functions thus leads to a large set of \codeterm{m-files} which increases complexity and may lead to confusion. If the delegated functionality -is used in multiple instances, it is advisable to do so. On the other -hand, when the delegated functionality is only used within the context -of another function \matlab{} allows to define +is used in multiple instances, it is still advisable to do so. On the +other hand, when the delegated functionality is only used within the +context of another function \matlab{} allows to define \codeterm[function!local]{local functions} and \codeterm[function!nested]{nested functions} within the same file. Listing \ref{localfunctions} shows an example of a local @@ -336,17 +338,19 @@ function definition. \pagebreak[3] \lstinputlisting[label=localfunctions, caption={Example for local functions.}]{calculateSines.m} -Local function live in the same \codeterm{m-file} as the main function -and are only available in this context. Each local function has its -own \codeterm{scope}, that is, the local function can not access (read -or write) variables of the calling function. +\emph{Local function} live in the same \codeterm{m-file} as the main +function and are only available in this context. Each local function +has its own \codeterm{scope}, that is, the local function can not +access (read or write) variables of the calling function. Interaction +with the local function requires to pass all required arguments and to +take care of the return values of the function. -This is different in so called \codeterm[function!nested]{nested - functions}. These are defined within the body of the parent function -(between the keywords \code{function} and \code{end}) and have full -access to all variables defined in the parent function. Working (in -particular changing) the parent's variables is handy on the one side, -but is also risky. One should take care when defining nested functions. +\emph{Nested functions} are different in this respect. They are +defined within the body of the parent function (between the keywords +\code{function} and \code{end}) and have full access to all variables +defined in the parent function. Working (in particular changing) the +parent's variables is handy on the one side, but is also risky. One +should take care when defining nested functions. \section{Specifics when using scripts} @@ -355,7 +359,7 @@ A similar problem as with nested function arises when using scripts become available in the global \codeterm{Workspace}. There is the risk of name conflicts, that is, a called sub-script redefines or uses the same variable name and may \emph{silently} change its content. The -user will not be notified by this change and the calling script may +user will not be notified about this change and the calling script may expect a completely different content. Bugs that are based on such mistakes are hard to find since the program itself looks perfectly fine. @@ -379,12 +383,15 @@ provides important information to track and fix the bug. \item Scripts should work independently of existing variables in the global workspace. - \item It is advisable to start a script with deleting variables - (\code{clear}) from the workspace and most of the times it is also - good to close all open figures (\code{close all}). + \item Often it is advisable to start a script with deleting + variables (\code{clear}) from the workspace and most of the times + it is also good to close all open figures (\code{close all}). Be + careful if a the respective script has been called by another one. \item Clean up the workspace at the end of a script. Delete (\code{clear}) all variables that are no longer needed. + + \item Consider to write functions instead of scripts. \end{itemize} \end{important} diff --git a/regression/lecture/regression.tex b/regression/lecture/regression.tex index d94374c..c60b926 100644 --- a/regression/lecture/regression.tex +++ b/regression/lecture/regression.tex @@ -86,7 +86,7 @@ large deviations. $f_{cost}(\{(x_i, y_i)\}|\{y^{est}_i\})$ is a so called \enterm{objective function} or \enterm{cost function}. We aim to adapt the model parameters to minimize the error (mean square error) and -thus the \emph{objective function}. In Chapter~\ref{maximumlikelihood} +thus the \emph{objective function}. In Chapter~\ref{maximumlikelihoodchapter} we will show that the minimization of the mean square error is equivalent to maximizing the likelihood that the observations originate from the model (assuming a normal distribution of the data @@ -270,7 +270,7 @@ The gradient is given by partial derivatives (Box~\ref{partialderivativebox}) with respect to the parameters $m$ and $b$ of the linear equation. There is no need to calculate it analytically but it can be estimated from the partial derivatives -using the difference quotient (Box~\ref{differentialquotient}) for +using the difference quotient (Box~\ref{differentialquotientbox}) for small steps $\Delta m$ und $\Delta b$. For example the partial derivative with respect to $m$: @@ -383,7 +383,7 @@ introduce how this can be done without using the gradient descent Problems that involve nonlinear computations on parameters, e.g. the rate $\lambda$ in the exponential function $f(x;\lambda) = -\exp(\lambda x)$, do not have an analytical solution. To find minima +e^{\lambda x}$, do not have an analytical solution. To find minima in such functions numerical methods such as the gradient descent have to be applied. @@ -396,7 +396,7 @@ objective functions while more specialized functions are specifically designed for optimizations in the least square error sense \matlabfun{lsqcurvefit()}. -\newpage +%\newpage \begin{important}[Beware of secondary minima!] Finding the absolute minimum is not always as easy as in the case of the linear equation. Often, the error surface has secondary or local diff --git a/scientificcomputing-script.tex b/scientificcomputing-script.tex index 22bbed8..03d8c8e 100644 --- a/scientificcomputing-script.tex +++ b/scientificcomputing-script.tex @@ -67,7 +67,7 @@ \lstset{inputpath=regression/code} \include{regression/lecture/regression} -\setboolean{showexercisesolutions}{true} +\setboolean{showexercisesolutions}{false} \graphicspath{{likelihood/lecture/}{likelihood/lecture/figures/}} \lstset{inputpath=likelihood/code} \include{likelihood/lecture/likelihood} diff --git a/statistics/lecture/displayunivariatedata.py b/statistics/lecture/displayunivariatedata.py index 546d4ff..55130a8 100644 --- a/statistics/lecture/displayunivariatedata.py +++ b/statistics/lecture/displayunivariatedata.py @@ -53,7 +53,7 @@ ax.annotate('mean plus\nstd. dev.', arrowprops=dict(arrowstyle="->", relpos=(0.5,0.0), connectionstyle="angle3,angleA=-60,angleB=80") ) -ax = fig.add_axes([xpos, ypos, width, height], axis_bgcolor='none') +ax = fig.add_axes([xpos, ypos, width, height]) ax.spines['right'].set_visible(False) ax.spines['top'].set_visible(False) ax.spines['left'].set_visible(False) @@ -92,7 +92,7 @@ ax.annotate('median', arrowprops=dict(arrowstyle="->", relpos=(0.8,0.0), connectionstyle="angle3,angleA=-60,angleB=20") ) -ax = fig.add_axes([xpos+width+0.03, ypos, 0.98-(xpos+width+0.03), height], axis_bgcolor='none') +ax = fig.add_axes([xpos+width+0.03, ypos, 0.98-(xpos+width+0.03), height]) ax.spines['right'].set_visible(False) ax.spines['top'].set_visible(False) ax.xaxis.set_ticks_position('bottom')