Merge branch 'master' of raven.am28.uni-tuebingen.de:scientificComputing

designpattern/lecture/Makefile

@@ -0,0 +1,29 @@
+BASENAME=designpattern
+
+PYFILES=$(wildcard *.py)
+PYPDFFILES=$(PYFILES:.py=.pdf)
+
+all : pdf 
+
+# script:
+pdf : $(BASENAME)-chapter.pdf
+
+$(BASENAME)-chapter.pdf : $(BASENAME)-chapter.tex $(BASENAME).tex $(PYPDFFILES)
+	pdflatex -interaction=scrollmode $< | tee /dev/stderr | fgrep -q "Rerun to get cross-references right" && pdflatex -interaction=scrollmode $< || true
+
+$(PYPDFFILES) : %.pdf : %.py
+	python $<
+
+clean :
+	rm -f *~ 
+	rm -f $(BASENAME).aux $(BASENAME).log
+	rm -f $(BASENAME)-chapter.aux $(BASENAME)-chapter.log $(BASENAME)-chapter.out
+	rm -f $(PYPDFFILES) $(GPTTEXFILES)
+
+cleanall : clean
+	rm -f $(BASENAME)-chapter.pdf
+
+watchpdf :
+	while true; do ! make -q pdf && make pdf; sleep 0.5; done
+
+

designpattern/lecture/designpattern-chapter.tex

@@ -0,0 +1,225 @@
+\documentclass[12pt]{report}
+
+%%%%% title %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+\title{\tr{Introduction to Scientific Computing}{Einf\"uhrung in die wissenschaftliche Datenverarbeitung}}
+\author{Jan Grewe \& Jan Benda\\Abteilung Neuroethologie\\[2ex]\includegraphics[width=0.3\textwidth]{UT_WBMW_Rot_RGB}}
+\date{WS 15/16}
+
+%%%% language %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+% \newcommand{\tr}[2]{#1}  % en
+% \usepackage[english]{babel}
+\newcommand{\tr}[2]{#2}  % de
+\usepackage[german]{babel}
+
+%%%%% packages %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+\usepackage{pslatex}   % nice font for pdf file
+\usepackage[breaklinks=true,bookmarks=true,bookmarksopen=true,pdfpagemode=UseNone,pdfstartview=FitH,colorlinks=true,citecolor=blue]{hyperref}
+
+%%%% layout %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+\usepackage[left=25mm,right=25mm,top=20mm,bottom=30mm]{geometry}
+\setcounter{tocdepth}{1}
+
+%%%%% section style %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+\usepackage[sf,bf,it,big,clearempty]{titlesec}
+\setcounter{secnumdepth}{1}
+
+
+%%%%% units %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+\usepackage[mediumspace,mediumqspace,Gray]{SIunits}      % \ohm, \micro
+ 
+
+%%%%% figures %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+\usepackage{graphicx}
+\usepackage{xcolor}
+\pagecolor{white}
+
+\newcommand{\ruler}{\par\noindent\setlength{\unitlength}{1mm}\begin{picture}(0,6)%
+  \put(0,4){\line(1,0){170}}%
+  \multiput(0,2)(10,0){18}{\line(0,1){4}}%
+  \multiput(0,3)(1,0){170}{\line(0,1){2}}%
+  \put(0,0){\makebox(0,0){{\tiny 0}}}%
+  \put(10,0){\makebox(0,0){{\tiny 1}}}%
+  \put(20,0){\makebox(0,0){{\tiny 2}}}%
+  \put(30,0){\makebox(0,0){{\tiny 3}}}%
+  \put(40,0){\makebox(0,0){{\tiny 4}}}%
+  \put(50,0){\makebox(0,0){{\tiny 5}}}%
+  \put(60,0){\makebox(0,0){{\tiny 6}}}%
+  \put(70,0){\makebox(0,0){{\tiny 7}}}%
+  \put(80,0){\makebox(0,0){{\tiny 8}}}%
+  \put(90,0){\makebox(0,0){{\tiny 9}}}%
+  \put(100,0){\makebox(0,0){{\tiny 10}}}%
+  \put(110,0){\makebox(0,0){{\tiny 11}}}%
+  \put(120,0){\makebox(0,0){{\tiny 12}}}%
+  \put(130,0){\makebox(0,0){{\tiny 13}}}%
+  \put(140,0){\makebox(0,0){{\tiny 14}}}%
+  \put(150,0){\makebox(0,0){{\tiny 15}}}%
+  \put(160,0){\makebox(0,0){{\tiny 16}}}%
+  \put(170,0){\makebox(0,0){{\tiny 17}}}%
+  \end{picture}\par}
+
+% figures:
+\setlength{\fboxsep}{0pt}
+\newcommand{\texpicture}[1]{{\sffamily\footnotesize\input{#1.tex}}}
+%\newcommand{\texpicture}[1]{\fbox{\sffamily\footnotesize\input{#1.tex}}}
+%\newcommand{\texpicture}[1]{\setlength{\fboxsep}{2mm}\fbox{#1}}
+%\newcommand{\texpicture}[1]{}
+\newcommand{\figlabel}[1]{\textsf{\textbf{\large \uppercase{#1}}}}
+
+% maximum number of floats:
+\setcounter{topnumber}{2}
+\setcounter{bottomnumber}{0}
+\setcounter{totalnumber}{2}
+
+% float placement fractions:
+\renewcommand{\textfraction}{0.2}
+\renewcommand{\topfraction}{0.8}
+\renewcommand{\bottomfraction}{0.0}
+\renewcommand{\floatpagefraction}{0.5}
+
+% spacing for floats:
+\setlength{\floatsep}{12pt plus 2pt minus 2pt}
+\setlength{\textfloatsep}{20pt plus 4pt minus 2pt}
+\setlength{\intextsep}{12pt plus 2pt minus 2pt}
+
+% spacing for a floating page:
+\makeatletter
+  \setlength{\@fptop}{0pt}
+  \setlength{\@fpsep}{8pt plus 2.0fil}
+  \setlength{\@fpbot}{0pt plus 1.0fil}
+\makeatother
+
+% rules for floats:
+\newcommand{\topfigrule}{\vspace*{10pt}{\hrule height0.4pt}\vspace*{-10.4pt}}
+\newcommand{\bottomfigrule}{\vspace*{-10.4pt}{\hrule height0.4pt}\vspace*{10pt}}
+
+% captions:
+\usepackage[format=plain,singlelinecheck=off,labelfont=bf,font={small,sf}]{caption}
+
+% put caption on separate float:
+\newcommand{\breakfloat}{\end{figure}\begin{figure}[t]}
+
+% references to panels of a figure within the caption:
+\newcommand{\figitem}[1]{\textsf{\bfseries\uppercase{#1}}}
+% references to figures:
+\newcommand{\panel}[1]{\textsf{\uppercase{#1}}}
+\newcommand{\fref}[1]{\textup{\ref{#1}}}
+\newcommand{\subfref}[2]{\textup{\ref{#1}}\,\panel{#2}}
+% references to figures in normal text:
+\newcommand{\fig}{Fig.}
+\newcommand{\Fig}{Figure}
+\newcommand{\figs}{Figs.}
+\newcommand{\Figs}{Figures}
+\newcommand{\figref}[1]{\fig~\fref{#1}}
+\newcommand{\Figref}[1]{\Fig~\fref{#1}}
+\newcommand{\figsref}[1]{\figs~\fref{#1}}
+\newcommand{\Figsref}[1]{\Figs~\fref{#1}}
+\newcommand{\subfigref}[2]{\fig~\subfref{#1}{#2}}
+\newcommand{\Subfigref}[2]{\Fig~\subfref{#1}{#2}}
+\newcommand{\subfigsref}[2]{\figs~\subfref{#1}{#2}}
+\newcommand{\Subfigsref}[2]{\Figs~\subfref{#1}{#2}}
+% references to figures within bracketed text:
+\newcommand{\figb}{Fig.}
+\newcommand{\figsb}{Figs.}
+\newcommand{\figrefb}[1]{\figb~\fref{#1}}
+\newcommand{\figsrefb}[1]{\figsb~\fref{#1}}
+\newcommand{\subfigrefb}[2]{\figb~\subfref{#1}{#2}}
+\newcommand{\subfigsrefb}[2]{\figsb~\subfref{#1}{#2}}
+
+% references to tables:
+\newcommand{\tref}[1]{\textup{\ref{#1}}}
+% references to tables in normal text:
+\newcommand{\tab}{Tab.}
+\newcommand{\Tab}{Table}
+\newcommand{\tabs}{Tabs.}
+\newcommand{\Tabs}{Tables}
+\newcommand{\tabref}[1]{\tab~\tref{#1}}
+\newcommand{\Tabref}[1]{\Tab~\tref{#1}}
+\newcommand{\tabsref}[1]{\tabs~\tref{#1}}
+\newcommand{\Tabsref}[1]{\Tabs~\tref{#1}}
+% references to tables within bracketed text:
+\newcommand{\tabb}{Tab.}
+\newcommand{\tabsb}{Tab.}
+\newcommand{\tabrefb}[1]{\tabb~\tref{#1}}
+\newcommand{\tabsrefb}[1]{\tabsb~\tref{#1}}
+
+
+%%%%% equation references %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%\newcommand{\eqref}[1]{(\ref{#1})}
+\newcommand{\eqn}{\tr{Eq}{Gl}.}
+\newcommand{\Eqn}{\tr{Eq}{Gl}.}
+\newcommand{\eqns}{\tr{Eqs}{Gln}.}
+\newcommand{\Eqns}{\tr{Eqs}{Gln}.}
+\newcommand{\eqnref}[1]{\eqn~\eqref{#1}}
+\newcommand{\Eqnref}[1]{\Eqn~\eqref{#1}}
+\newcommand{\eqnsref}[1]{\eqns~\eqref{#1}}
+\newcommand{\Eqnsref}[1]{\Eqns~\eqref{#1}}
+
+
+%%%%% listings %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+\usepackage{listings}
+\lstset{
+  inputpath=../code,
+  basicstyle=\ttfamily\footnotesize,
+  numbers=left,
+  showstringspaces=false,
+  language=Matlab,
+  commentstyle=\itshape\color{darkgray},
+  keywordstyle=\color{blue},
+  stringstyle=\color{green},
+  backgroundcolor=\color{blue!10},
+  breaklines=true,
+  breakautoindent=true,
+  columns=flexible,
+  frame=single,
+  caption={\protect\filename@parse{\lstname}\protect\filename@base},
+  captionpos=t,
+  xleftmargin=1em,
+  xrightmargin=1em,
+  aboveskip=10pt
+}
+
+%%%%% math stuff: %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+\usepackage{amsmath}
+\usepackage{bm} 
+\usepackage{dsfont}
+\newcommand{\naZ}{\mathds{N}}
+\newcommand{\gaZ}{\mathds{Z}}
+\newcommand{\raZ}{\mathds{Q}}
+\newcommand{\reZ}{\mathds{R}}
+\newcommand{\reZp}{\mathds{R^+}}
+\newcommand{\reZpN}{\mathds{R^+_0}}
+\newcommand{\koZ}{\mathds{C}}
+
+
+%%%%% structure: %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+\usepackage{ifthen}
+
+\newcommand{\code}[1]{\texttt{#1}}
+
+\newcommand{\source}[1]{    
+  \begin{flushright}
+    \color{gray}\scriptsize \url{#1}
+  \end{flushright}
+}
+
+\newenvironment{definition}[1][]{\medskip\noindent\textbf{Definition}\ifthenelse{\equal{#1}{}}{}{ #1}:\newline}%
+  {\medskip}
+
+\newcounter{maxexercise} 
+\setcounter{maxexercise}{9}  % show listings up to exercise maxexercise
+\newcounter{theexercise} 
+\setcounter{theexercise}{1}
+\newenvironment{exercise}[1][]{\medskip\noindent\textbf{\tr{Exercise}{\"Ubung}
+  \arabic{theexercise}:}\newline \newcommand{\exercisesource}{#1}}%
+  {\ifthenelse{\equal{\exercisesource}{}}{}{\ifthenelse{\value{theexercise}>\value{maxexercise}}{}{\medskip\lstinputlisting{\exercisesource}}}\medskip\stepcounter{theexercise}}
+
+\graphicspath{{figures/}}
+
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+\begin{document} 
+
+\include{designpattern}
+
+\end{document}

designpattern/lecture/designpattern.tex

@@ -0,0 +1,166 @@
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+\chapter{Design Pattern}
+
+Beim Programmieren sind sich viel Codes in ihrer Grundstruktur sehr
+\"ahnlich. Viele Konstrukte kommen in den verschiedensten Kontexten
+immer wieder in \"ahnlicher Weise vor. In diesem Kapitel stellen wir
+einige dieser ``Design pattern'' zusammen.
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+\section{Plotten einer mathematischen Funktion}
+Eine mathematische Funktion ordnet einem beliebigen $x$-Wert einen
+$y$-Wert zu. Um eine solche Funktion zeichnen zu k\"onnen, m\"ussen
+wir uns eine Wertetabelle aus vielen $x$-Werten und den
+dazugeh\"origen Funktionswerten $y=f(x)$ erstellen. 
+
+Wir erstellen uns dazu einen Vektor mit geeigneten $x$-Werten, die von
+dem kleinsten bis zu dem gr\"o{\ss}ten $x$-Wert laufen, den wir
+plotten wollen. Die Schrittweite f\"ur die $x$-Werte w\"ahlen wir
+klein genug, um eine sch\"one glatte Kurve zu bekommen. F\"ur jeden
+Wert $x_i$ dieses Vektors berechnen wir den entsprechenden
+Funktionswert und erzeugen damit einen Vektor mit den $y$-Werten. Die
+Werte des $y$-Vektors k\"onnen dann gegen die Werte des $x$-Vektors
+geplottet werden.
+
+Folgende Programme berechnen und plotten die Funktion $f(x)=e^{-x^2}$:
+\begin{lstlisting}
+xmin = -1.0;
+xmax = 2.0;
+dx = 0.01;          % Schrittweite
+x = xmin:dx:xmax;   % Vektor mit x-Werten
+y = exp(-x.*x);     % keine for Schleife! '.*' fuer elementweises multiplizieren
+plot(x, y);
+\end{lstlisting}
+
+\begin{lstlisting}
+x = -1:0.01:2;      % Vektor mit x-Werten
+y = exp(-x.*x);     % keine for Schleife! '.*' fuer elementweises multiplizieren
+plot(x, y);
+\end{lstlisting}
+
+\begin{lstlisting}
+x = -1:0.01:2;      % Vektor mit x-Werten
+plot(x, exp(-x.*x));
+\end{lstlisting}
+
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+\section{Skalieren und Verschieben nicht nur von Zufallszahlen}
+Zufallsgeneratoren geben oft nur Zufallszahlen mit festen Mittelwerten
+und Standardabweichungen (auch Skalierungen) zur\"uck. Multiplikation
+mit einem Faktor skaliert die Standardabweichung und Addition einer Zahl
+verschiebt den Mittelwert.
+
+\begin{lstlisting}
+% 100 random numbers draw from a Gaussian distribution with mean 0 and standard deviation 1.
+x = randn(100, 1);
+
+% 100 random numbers drawn from a Gaussian distribution with mean 4.8 and standard deviation 2.3.
+mu = 4.8;
+sigma = 2.3;
+y = randn(100, 1)*sigma + mu;
+\end{lstlisting}
+
+Das ist manchmal auch sinnvoll f\"ur \code{zeros} oder \code{ones}:
+\begin{lstlisting}
+x = -1:0.01:2;      % Vektor mit x-Werten
+plot(x, exp(-x.*x));
+% Plotte f\"ur die gleichen x-Werte eine Linie mit y=0.8:
+plot(x, zeros(size(x))+0.8);
+% ... Linie mit y=0.5:
+plot(x, ones(size(x))*0.5);
+\end{lstlisting}
+
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+\section{for Schleifen \"uber Vektoren}
+Manchmal m\"ochte man doch mit einer for-Schleife \"uber einen Vektor iterieren:
+\begin{lstlisting}
+x = [2:3:20];   % irgendein Vektor
+for i=1:length(x)
+  % Benutze den Wert des Vektors x an der Stelle des Indexes i:
+  do_something( x(i) );
+end
+\end{lstlisting}
+
+Wenn in der Schleife das Ergebnis in einen Vektor gespeichert werden soll,
+sollten wir uns vor der Schleife schon einen Vektor f\"ur die Ergebnisse
+erstellen:
+\begin{lstlisting}
+x = [2:3:20];        % irgendein Vektor
+y = zeros(size(x));  % Platz fuer die Ergebnisse
+for i=1:length(x)
+  % Schreibe den Rueckgabewert der Funktion get_something an die i-te
+  % Stelle von y:
+  y(i) = get_something( x(i) );
+end
+% jetzt koennen wir den Ergebnisvektor weiter bearbeiten:
+mean(y)
+\end{lstlisting}
+
+Die Berechnungen in der Schleife k\"onnen statt einer Zahl auch einen Vektor
+zur\"uckgeben. Wenn die L\"ange diese Vektors bekannt ist, dann kann vorher
+eine entsprechend gro{\ss}e Matrix angelegt werden:
+\begin{lstlisting}
+x = [2:3:20];        % irgendein Vektor
+y = zeros(length(x),10);  % Platz fuer die Ergebnisse
+for i=1:length(x)
+  % Schreibe den Rueckgabewert der Funktion get_something - jetzt ein
+  % Vektor mit 10 Elementen - in die i-te
+  % Zeile von y:
+  y(i,:) = get_something( x(i) );
+end
+% jetzt koennen wir die Ergebnismatrix weiter bearbeiten:
+mean(y, 1)
+\end{lstlisting}
+
+
+Alternativ k\"onnen die in der Schleife erzeugten Vektoren zu einem
+einzigen, durchgehenden Vektor zusammengestellt werden:
+\begin{lstlisting}
+x = [2:3:20];        % irgendein Vektor
+y = [];  % Leerer Vektor fuer die Ergebnisse
+for i=1:length(x)
+  % Die Funktion get_something gibt uns einen Vektor zurueck:
+  z = get_something( x(i) );
+  % dessen Inhalt h\"angen wir an unseren Ergebnissvektor an:
+  y = [y z(:)];
+end
+% jetzt koennen wir dem Ergebnisvektor weiter bearbeiten:
+mean(y)
+\end{lstlisting}
+
+
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+\section{Normierung von Histogrammen}
+Meistens sollten Histogramme normiert werden, damit sie vergleichbar
+mit anderen Histogrammen oder mit theoretischen
+Wahrscheinlichkeitsverteilungen werden.
+
+Die \code{histogram} Funktion macht das mit den entsprechenden Parametern automatisch:
+\begin{lstlisting}
+x = randn(100, 1);    % irgendwelche reellwertige Daten
+histogram(x, 'Normalization', 'pdf');
+\end{lstlisting}
+
+\begin{lstlisting}
+x = randi(6, 100, 1);    % irgendwelche integer Daten
+histogram(x, 'Normalization', 'probability');
+\end{lstlisting}
+
+So geht es aber auch:
+\begin{lstlisting}
+x = randn(100, 1);         % irgendwelche reellwertige Daten
+[h, b] = hist(x);          % Histogram berechnen
+h = h/sum(h)/(b(2)-b(1));  % normieren zu einer Wahrscheinlichkeitsdichte
+bar(b, h);                 % und plotten.
+\end{lstlisting}
+
+\begin{lstlisting}
+x = randi(6, 100, 1);    % irgendwelche integer Daten
+[h, b] = hist(x);        % Histogram berechnen
+h = h/sum(h);            % normieren zu Wahrscheinlichkeiten
+bar(b, h);               % und plotten.
+\end{lstlisting}

pointprocesses/code/counthist.m

@@ -1,4 +1,4 @@
-function [ counts, bins ] = counthist( spikes, w )
+function [counts, bins] = counthist(spikes, w)
 % computes count histogram and compare them  with Poisson distribution
 % spikes: a cell array of vectors of spike times
 % w: observation window duration for computing the counts
@@ -18,18 +18,17 @@ function [ counts, bins ] = counthist( spikes, w )
     end
     % histogram of spike counts:
     maxn = max( n );
-    [counts, bins ] = hist( n, 0:1:maxn+1 );
+    [counts, bins ] = hist( n, 0:1:maxn+10 );
     counts = counts / sum( counts );
     if nargout == 0
         bar( bins, counts );
         hold on;
         % Poisson distribution:
         rate = mean( r );
-        x = 0:1:20;
+        x = 0:1:maxn+10;
         l = rate*w;
         y = l.^x.*exp(-l)./factorial(x);
         plot( x, y, 'r', 'LineWidth', 3 );
-        xlim( [ 0 20 ] );
         hold off;
         xlabel( 'counts k' );
         ylabel( 'P(k)' );

pointprocesses/code/isihist.m

@@ -24,9 +24,9 @@ function isihist( isis, binwidth )
     misi = mean( isis );
     sdisi = std( isis );
     disi = sdisi^2.0/2.0/misi^3;
-    text( 0.5, 0.6, sprintf( 'mean=%.1f ms', 1000.0*misi ), 'Units', 'normalized' )
+    text( 0.95, 0.8, sprintf( 'mean=%.1f ms', 1000.0*misi ), 'Units', 'normalized', 'HorizontalAlignment', 'right' )
-    text( 0.5, 0.5, sprintf( 'std=%.1f ms', 1000.0*sdisi ), 'Units', 'normalized' )
+    text( 0.95, 0.7, sprintf( 'std=%.1f ms', 1000.0*sdisi ), 'Units', 'normalized', 'HorizontalAlignment', 'right' )
-    text( 0.5, 0.4, sprintf( 'CV=%.2f', sdisi/misi ), 'Units', 'normalized' )
+    text( 0.95, 0.6, sprintf( 'CV=%.2f', sdisi/misi ), 'Units', 'normalized', 'HorizontalAlignment', 'right' )
     %text( 0.5, 0.3, sprintf( 'D=%.1f Hz', disi ), 'Units', 'normalized' )
 end
 

pointprocesses/code/lifadaptspikes.m

@@ -3,8 +3,8 @@ function spikes = lifadaptspikes( trials, input, tmaxdt, D, tauadapt, adaptincr
 % with an adaptation current
 % trials: the number of trials to be generated
 % input: the stimulus either as a single value or as a vector
-% tmaxdt: in case of a single value stimulus the duration of a trial
+% tmaxdt: in case of a single value stimulus: the duration of a trial
-%         in case of a vector as a stimulus the time step
+%         in case of a vector as a stimulus: the time step
 % D: the strength of additive white noise
 % tauadapt: adaptation time constant
 % adaptincr: adaptation strength

pointprocesses/code/lifspikes.m

@@ -2,8 +2,8 @@ function spikes = lifspikes( trials, input, tmaxdt, D )
 % Generate spike times of a leaky integrate-and-fire neuron
 % trials: the number of trials to be generated
 % input: the stimulus either as a single value or as a vector
-% tmaxdt: in case of a single value stimulus the duration of a trial
+% tmaxdt: in case of a single value stimulus: the duration of a trial
-%         in case of a vector as a stimulus the time step
+%         in case of a vector as a stimulus: the time step
 % D: the strength of additive white noise
     
     tau = 0.01;

pointprocesses/code/lifspikesoustim.m

@@ -0,0 +1,20 @@
+function spikes = lifspikesoustim(trials, tmax, D, Iou, Dou, tauou )
+% Generate spike times of a leaky integrate-and-fire neuron with frozen
+% Ohrnstein-Uhlenbeck stimulus
+% trials: the number of trials to be generated
+% tmax: the duration of a trial
+% D: the strength of additive white noise
+% Iou: the mean input
+% Dou: noise strength of the frozen OU noise
+% tauou: time constant of the OU noise
+    
+    dt = 1e-4;
+    input = zeros(round(tmax/dt), 1);
+    n = 0.0;
+    noise = sqrt(2.0*Dou)*randn(length(input), 1)/sqrt(dt);
+    for i=1:length(noise)
+        n = n + ( - n + noise(i))*dt/tauou;
+        input(i) = Iou + n;
+    end
+    spikes = lifspikes(trials, input, dt, D );
+end

pointprocesses/code/pifouspikes.m

@@ -1,9 +1,10 @@
 function spikes = pifouspikes( trials, input, tmaxdt, D, outau )
 % Generate spike times of a perfect integrate-and-fire neuron
+% with Ornstein-Uhlenbeck noise (colored noise).
 % trials: the number of trials to be generated
 % input: the stimulus either as a single value or as a vector
-% tmaxdt: in case of a single value stimulus the duration of a trial
+% tmaxdt: in case of a single value stimulus: the duration of a trial
-%         in case of a vector as a stimulus the time step
+%         in case of a vector as a stimulus: the time step
 % D: the strength of additive white noise
 % outau: time constant of the colored noise
     

pointprocesses/code/plotcounthist.m

@@ -0,0 +1,24 @@
+w = 0.1;
+cmax = 8;
+pmax = 0.5;
+subplot(1, 3, 1);
+counthist(poissonspikes, w);
+xlim([0 cmax])
+set(gca, 'XTick', 0:2:cmax)
+ylim([0 pmax])
+title('Poisson');
+
+subplot(1, 3, 2);
+counthist(pifouspikes, w);
+xlim([0 cmax])
+set(gca, 'XTick', 0:2:cmax)
+ylim([0 pmax])
+title('PIF OU');
+
+subplot(1, 3, 3);
+counthist(lifadaptspikes, w);
+xlim([0 cmax])
+set(gca, 'XTick', 0:2:cmax)
+ylim([0 pmax])
+title('LIF adapt');
+savefigpdf(gcf, 'counthist.pdf', 20, 7);

pointprocesses/code/plotisih.m

@@ -0,0 +1,19 @@
+maxisi = 300.0;
+subplot(1, 3, 1);
+poissonisis = isis(poissonspikes);
+isihist(poissonisis, 0.001);
+xlim([0, maxisi])
+title('Poisson');
+
+subplot(1, 3, 2);
+pifouisis = isis(pifouspikes);
+isihist(pifouisis, 0.001);
+xlim([0, maxisi])
+title('PIF OU');
+
+subplot(1, 3, 3);
+lifadaptisis = isis(lifadaptspikes);
+isihist(lifadaptisis, 0.001);
+xlim([0, maxisi])
+title('LIF adapt');
+savefigpdf(gcf, 'isihist.pdf', 20, 7);

pointprocesses/code/plotserialcorr.m

@@ -0,0 +1,17 @@
+maxlag = 10;
+rrange = [-0.5, 1.05];
+subplot(1, 3, 1);
+isiserialcorr(poissonisis, maxlag);
+ylim(rrange)
+title('Poisson');
+
+subplot(1, 3, 2);
+isiserialcorr(pifouisis, maxlag);
+ylim(rrange)
+title('PIF OU');
+
+subplot(1, 3, 3);
+isiserialcorr(lifadaptisis, maxlag);
+ylim(rrange)
+title('LIF adapt');
+savefigpdf(gcf, 'serialcorr.pdf', 20, 7);

pointprocesses/code/plotspikeraster.m

@@ -0,0 +1,13 @@
+subplot(1, 3, 1);
+spikeraster(poissonspikes, 1.0);
+title('Poisson');
+
+subplot(1, 3, 2);
+spikeraster(pifouspikes, 1.0);
+title('PIF OU');
+
+subplot(1, 3, 3);
+spikeraster(lifadaptspikes, 1.0);
+title('LIF adapt');
+
+savefigpdf(gcf, 'spikeraster.pdf', 15, 5);

pointprocesses/code/poissonspikes.m

@@ -12,9 +12,9 @@ function spikes = poissonspikes( trials, rate, tmax )
         p = 0.1
         dt = p/rate;
     end
-    spikes = cell( trials, 1 );
+    spikes = cell(trials, 1);
     for k=1:trials
-        x = rand( 1, round(tmax/dt) ); % uniform random numbers for each bin
+        x = rand(round(tmax/dt), 1); % uniform random numbers for each bin
-        spikes{k} = find( x < p ) * dt;
+        spikes{k} = find(x < p) * dt;
     end
 end

pointprocesses/code/psth.m

@@ -0,0 +1,14 @@
+function p = psth(spikes, dt, tmax)
+% plots a PSTH of the spikes with binwidth dt
+    t = 0.0:dt:tmax+dt;
+    p = zeros(1, length(t));
+    for k=1:length(spikes)
+        times = spikes{k};
+        [h, b] = hist(times, t);
+        p = p + h;
+    end
+    p = p/length(spikes)/dt;
+    t(end) = [];
+    p(end) = [];
+    plot(t, p);
+end

pointprocesses/code/spikeraster.m

@@ -1,15 +1,24 @@
-function spikeraster( spikes )
+function spikeraster(spikes, tmax)
 % Display a spike raster of the spike times given in spikes.
 % spikes: a cell array of vectors of spike times
+% tmax: plot spike raster upto tmax seconds
 
 ntrials = length(spikes);
 for k = 1:ntrials
-    times = 1000.0*spikes{k};  % conversion to ms
+    times = spikes{k};
+    times = times(times<tmax);
+    if tmax < 1.5
+        times = 1000.0*times; % conversion to ms
+    end
     for i = 1:length( times )
       line([times(i) times(i)],[k-0.4 k+0.4], 'Color', 'k' ); 
     end
 end
-xlabel( 'Time [ms]' );
+if tmax < 1.5
+    xlabel( 'Time [ms]' );
+else
+    xlabel( 'Time [s]' );
+end
 ylabel( 'Trials');
 ylim( [ 0.3 ntrials+0.7 ] )
 

pointprocesses/code/spikerate.m

@@ -0,0 +1,12 @@
+function r = spikerate(spikes, duration)
+% returns the average spike rate of the spikes
+% for the first duration seconds
+% spikes: a cell array of vectors of spike times
+
+    rates = zeros(length(spikes),1);
+    for k = 1:length(spikes)
+        times = spikes{k};
+        rates(k) = sum(times<duration)/duration;
+    end
+    r = mean(rates);
+end

pointprocesses/exercises/hompoissonspikes.m

@@ -14,6 +14,6 @@ function spikes = hompoissonspikes( trials, rate, tmax )
     x = rand( trials, ceil(tmax/dt) );
     spikes = cell( trials, 1 );
     for k=1:trials
-        spikes{k} = find( x(k,:) >= 1.0-p ) * dt;
+        spikes{k} = find( x(k,:) < p ) * dt;
     end
 end

pointprocesses/exercises/pointprocesses01.tex

@@ -11,11 +11,11 @@
 \usepackage[left=20mm,right=20mm,top=25mm,bottom=25mm]{geometry}
 \pagestyle{headandfoot}
 \ifprintanswers
-\newcommand{\stitle}{: L\"osungen}
+\newcommand{\stitle}{L\"osungen}
 \else
-\newcommand{\stitle}{}
+\newcommand{\stitle}{\"Ubung}
 \fi
-\header{{\bfseries\large \"Ubung 6\stitle}}{{\bfseries\large Statistik}}{{\bfseries\large 27. Oktober, 2015}}
+\header{{\bfseries\large \stitle}}{{\bfseries\large Punktprozesse}}{{\bfseries\large 27. Oktober, 2015}}
 \firstpagefooter{Prof. Dr. Jan Benda}{Phone: 29 74573}{Email:
 jan.benda@uni-tuebingen.de}
 \runningfooter{}{\thepage}{}
@@ -89,114 +89,110 @@ jan.benda@uni-tuebingen.de}
 \begin{questions}
 
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-\question \qt{Homogeneous Poisson process}
+  \question \qt{Statistik von Spiketrains}
-We use the Poisson process to generate spike trains on which we can test and imrpove some 
+  In Ilias findet ihr die Dateien \code{poisson.mat},
-standard analysis functions.
+  \code{pifou.mat}, und \code{lifadapt.mat}.  Jede dieser Dateien
-
+  enth\"alt mehrere Trials von Spiketrains von einer bestimmten Art
-A homogeneous Poisson process of rate $\lambda$ (measured in Hertz) is a point process
+  von Neuron. Die Spikezeiten sind in Sekunden gemessen.
-where the probability of an event is independent of time $t$ and independent of previous events.
+
-The probability $P$ of an event within a bin of width $\Delta t$ is 
+  Mit den folgenden Aufgaben wollen wir die Statistik der Spiketrains
-\[ P = \lambda \cdot \Delta t \]
+  der drei Neurone miteinander vergleichen.
-for sufficiently small $\Delta t$.
+  \begin{parts}
-\begin{parts}
+    \part Lade die Spiketrains aus den drei Dateien. Achte darauf, dass sie verschiedene
+    Variablen\-namen bekommen.
+    \begin{solution}
+      \begin{lstlisting}
+        clear all
+        % not so good:
+        load poisson.mat
+        whos
+        poissonspikes = spikes;
+        load pifou.mat;
+        pifouspikes = spikes;
+        load lifadapt.mat;
+        lifadaptspikes = spikes;
+        clear spikes;
+        % better:
+        clear all
+        x = load( 'poisson.mat' );
+        poissonspikes = x.spikes;
+        x = load( pifou.mat' );
+        pifouspikes = x.spikes;
+        x = load( 'lifadapt.mat' );
+        lifadaptspikes = x.spikes;
+      \end{lstlisting}
+    \end{solution}
+    
+    \part Schreibe eine Funktion, die die Spikezeiten der ersten
+    $t_{max}$ Sekunden in einem Rasterplot visualisiert.  In jeder
+    Zeile des Rasterplots wird ein Spiketrain dargestellt. Jeder
+    einzelne Spike wird als senkrechte Linie zu der Zeit des
+    Auftretens des Spikes geplottet. Benutze die Funktion, um die
+    Spikeraster der ersten 1\,s der drei Neurone nebeneinander zu plotten.
+    \begin{solution}
+      \lstinputlisting{../code/spikeraster.m}
+      \lstinputlisting{../code/plotspikeraster.m}
+      \mbox{}\\[-3ex]
+      \colorbox{white}{\includegraphics[width=1\textwidth]{spikeraster}}
+    \end{solution}
   
-  \part Write a function that generates $n$ homogeneous Poisson spike trains of a given duration $T_{max}$
+    \part Schreibe eine Funktion, die einen einzigen Vektor mit den
-  with rate $\lambda$.
+    Interspikeintervallen aller Trials von Spikezeiten zur\"uckgibt.
-  \begin{solution}
+    \begin{solution}
-    \lstinputlisting{hompoissonspikes.m}
+      \lstinputlisting{../code/isis.m}
-  \end{solution}
+    \end{solution}
+    
+    \part Schreibe eine Funktion, die ein normiertes Histogramm aus
+    einem Vektor von Interspikeintervallen, gegeben in Sekunden,
+    berechnet und dieses mit richtiger Achsenbeschriftung plottet.
+    Die Interspikeintervalle sollen dabei in Millisekunden angegeben
+    werden. Die Funktion soll zus\"atzlich den Mittelwert, die
+    Standardabweichung, und den Variationskoeffizienten der
+    Interspikeintervalle berechnen und diese im Plot mit angeben.
+    
+    Benutze die vorherige und diese Funktion, um die
+    Interspikeintervall Verteilung der drei Neurone zu vergleichen.
+    \begin{solution}
+      \lstinputlisting{../code/isihist.m}
+      \lstinputlisting{../code/plotisih.m}
+      \mbox{}\\[-3ex]
+      \colorbox{white}{\includegraphics[width=1\textwidth]{isihist}}
+    \end{solution}
   
-  \part Using this function, generate a few trials and display them in a raster plot.
+    \part Schreibe eine Funktion, die die seriellen Korrelationen der
-  \begin{solution}
+    Interspikeintervalle f\"ur Lags bis zu \code{maxlag} berechnet und
-    \lstinputlisting{../code/spikeraster.m}
+    plottet.  Die Seriellen Korrelationen $\rho_k$ f\"ur Lag $k$ der
-    \begin{lstlisting}
+    Interspikeintervalle $T_i$ sind die Korrelationskoeffizienten
-      spikes = hompoissonspikes( 10, 100.0, 0.5 );
+    zwischen den Interspikeintervallen $T_i$ und den um das Lag $k$
-      spikeraster( spikes )
+    verschobenen Intervallen $T_{i+k}$:
-    \end{lstlisting}
+    \[ \rho_k = \frac{\langle (T_{i+k} - \langle T \rangle)(T_i -
-    \mbox{}\\[-3ex]
+      \langle T \rangle) \rangle}{\langle (T_i - \langle T
-    \colorbox{white}{\includegraphics[width=0.7\textwidth]{poissonraster100hz}}
+      \rangle)^2\rangle} = \frac{{\rm cov}(T_{i+k}, T_i)}{{\rm
-  \end{solution}
+        var}(T_i)} = {\rm corr}(T_{i+k}, T_i) \] 
+
+    Benutze diese Funktion, um die Interspikeintervall-Korrelationen
+    der drei Neurone zu vergleichen.
+    \begin{solution}
+      \lstinputlisting{../code/isiserialcorr.m}
+      \lstinputlisting{../code/plotserialcorr.m}
+      \colorbox{white}{\includegraphics[width=1\textwidth]{serialcorr}}
+    \end{solution}
   
-  \part Write a function that extracts a single vector of interspike intervals
+    \part Schreibe eine Funktion, die aus Spikezeiten
-  from the spike times returned by the first function.
+    Histogramme aus der Anzahl von Spikes, die in Fenstern gegebener L\"ange $W$
-  \begin{solution}
+    gez\"ahlt werden, erzeugt und plottet.
-    \lstinputlisting{../code/isis.m}
+
-  \end{solution}
+    Wende diese Funktion auf die drei
+    Datens\"atze an. Probiere verschiedene Fenstergr\"o{\ss}en $W$ aus.
+    \begin{solution}
+      \lstinputlisting{../code/counthist.m}
+      \lstinputlisting{../code/plotcounthist.m}
+      \colorbox{white}{\includegraphics[width=1\textwidth]{counthist}}
+    \end{solution}
+
+    
+  \end{parts}
   
-  \part Write a function that plots the interspike-interval histogram
-  from a vector of interspike intervals. The function should also
-  compute the mean, the standard deviation, and the CV of the intervals
-  and display the values in the plot.
-  \begin{solution}
-    \lstinputlisting{../code/isihist.m}
-  \end{solution}
-  
-  \part Compute histograms for Poisson spike trains with rate
-  $\lambda=100$\,Hz. Play around with $T_{max}$ and $n$ and the bin width
-  (start with 1\,ms) of the histogram.
-  How many
-  interspike intervals do you approximately need to get a ``nice''
-  histogram? How long do you need to record from the neuron?
-  \begin{solution}
-    About 5000 intervals for 25 bins. This corresponds to a $5000 / 100\,\hertz = 50\,\second$ recording
-    of a neuron firing with 100\,\hertz.
-  \end{solution}
-  
-  \part Compare the histogram with the true distribution of intervals $T$ of the Poisson process
-  \[ p(T) = \lambda e^{-\lambda T} \]
-  for various rates $\lambda$.
-  \begin{solution}
-    \lstinputlisting{hompoissonisih.m}
-    \colorbox{white}{\includegraphics[width=0.48\textwidth]{poissonisih100hz}}
-    \colorbox{white}{\includegraphics[width=0.48\textwidth]{poissonisih20hz}}
-  \end{solution}
-  
-  \part What happens if you make the bin width of the histogram smaller than $\Delta t$ 
-  used for generating the Poisson spikes?
-  \begin{solution}
-    The bins between the discretization have zero entries. Therefore
-    the other ones become higher than they should be.
-  \end{solution}
-  
-  \part Plot the mean interspike interval, the corresponding standard deviation, and the CV
-  as a function of the rate $\lambda$ of the Poisson process.
-  Compare the ../code with the theoretical expectations for the dependence on $\lambda$.
-  \begin{solution}
-    \lstinputlisting{hompoissonisistats.m}
-    \colorbox{white}{\includegraphics[width=0.98\textwidth]{poissonisistats}}
-  \end{solution}
-  
-  \part Write a function that computes serial correlations for the interspike intervals
-  for a range of lags.
-  The serial correlations $\rho_k$ at lag $k$ are defined as
-  \[ \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)} \]
-  Use this function to show that interspike intervals of Poisson spikes are independent.
-  \begin{solution}
-    \lstinputlisting{../code/isiserialcorr.m}
-    \colorbox{white}{\includegraphics[width=0.98\textwidth]{poissonserial100hz}}
-  \end{solution}
-  
-  \part Write a function that generates from spike times 
-  a histogram of spike counts in a count window of given duration $W$.
-  The function should also plot the Poisson distribution
-  \[ P(k) = \frac{(\lambda W)^ke^{\lambda W}}{k!} \]
-  for the rate $\lambda$ determined from the spike trains.
-  \begin{solution}
-    \lstinputlisting{../code/counthist.m}
-    \colorbox{white}{\includegraphics[width=0.48\textwidth]{poissoncounthistdist100hz10ms}}
-    \colorbox{white}{\includegraphics[width=0.48\textwidth]{poissoncounthistdist100hz100ms}}
-  \end{solution}
-  
-  \part Write a function that computes mean count, variance of count and the corresponding Fano factor
-  for a range of count window durations. The function should generate tow plots: one plotting
-  the count variance against the mean, the other one the Fano factor as a function of the window duration.
-  \begin{solution}
-    \lstinputlisting{../code/fano.m}
-    \colorbox{white}{\includegraphics[width=0.98\textwidth]{poissonfano100hz}}
-  \end{solution}
-  
-\end{parts}
-
 \end{questions}
 
 \end{document}

pointprocesses/exercises/pointprocesses02.tex

@@ -0,0 +1,199 @@
+\documentclass[12pt,a4paper,pdftex]{exam}
+
+\usepackage[german]{babel}
+\usepackage{pslatex}
+\usepackage[mediumspace,mediumqspace,Gray]{SIunits}      % \ohm, \micro
+\usepackage{xcolor}
+\usepackage{graphicx}
+\usepackage[breaklinks=true,bookmarks=true,bookmarksopen=true,pdfpagemode=UseNone,pdfstartview=FitH,colorlinks=true,citecolor=blue]{hyperref}
+
+%%%%% layout %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+\usepackage[left=20mm,right=20mm,top=25mm,bottom=25mm]{geometry}
+\pagestyle{headandfoot}
+\ifprintanswers
+\newcommand{\stitle}{L\"osungen}
+\else
+\newcommand{\stitle}{\"Ubung}
+\fi
+\header{{\bfseries\large \stitle}}{{\bfseries\large Punktprozesse 2}}{{\bfseries\large 27. Oktober, 2015}}
+\firstpagefooter{Prof. Dr. Jan Benda}{Phone: 29 74573}{Email:
+jan.benda@uni-tuebingen.de}
+\runningfooter{}{\thepage}{}
+
+\setlength{\baselineskip}{15pt}
+\setlength{\parindent}{0.0cm}
+\setlength{\parskip}{0.3cm}
+\renewcommand{\baselinestretch}{1.15}
+
+%%%%% listings %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+\usepackage{listings}
+\lstset{
+  language=Matlab,
+  basicstyle=\ttfamily\footnotesize,
+  numbers=left,
+  numberstyle=\tiny,
+  title=\lstname,
+  showstringspaces=false,
+  commentstyle=\itshape\color{darkgray},
+  breaklines=true,
+  breakautoindent=true,
+  columns=flexible,
+  frame=single,
+  xleftmargin=1em,
+  xrightmargin=1em,
+  aboveskip=10pt
+}
+
+%%%%% math stuff: %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+\usepackage{amsmath}
+\usepackage{amssymb}
+\usepackage{bm} 
+\usepackage{dsfont}
+\newcommand{\naZ}{\mathds{N}}
+\newcommand{\gaZ}{\mathds{Z}}
+\newcommand{\raZ}{\mathds{Q}}
+\newcommand{\reZ}{\mathds{R}}
+\newcommand{\reZp}{\mathds{R^+}}
+\newcommand{\reZpN}{\mathds{R^+_0}}
+\newcommand{\koZ}{\mathds{C}}
+
+%%%%% page breaks %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+\newcommand{\continue}{\ifprintanswers%
+\else
+\vfill\hspace*{\fill}$\rightarrow$\newpage%
+\fi}
+\newcommand{\continuepage}{\ifprintanswers%
+\newpage
+\else
+\vfill\hspace*{\fill}$\rightarrow$\newpage%
+\fi}
+\newcommand{\newsolutionpage}{\ifprintanswers%
+\newpage%
+\else
+\fi}
+
+%%%%% new commands %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+\newcommand{\qt}[1]{\textbf{#1}\\}
+\newcommand{\pref}[1]{(\ref{#1})}
+\newcommand{\extra}{--- Zusatzaufgabe ---\ \mbox{}}
+\newcommand{\code}[1]{\texttt{#1}}
+
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+\begin{document}
+
+\input{instructions}
+
+
+\begin{questions}
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+\question \qt{Homogener Poisson Prozess}
+Wir wollen den homogenen Poisson Prozess benutzen um Spikes zu generieren,
+mit denen wir die Analysfunktionen des vorherigen \"Ubungszettel \"uberpr\"ufen k\"onnen.
+
+Ein homogener Poisson Prozess mit der Rate $\lambda$ (gemessen in Hertz) ist ein Punktprozess,
+bei dem die Wahrschienlichkeit eines Ereignisses unabh\"angig von der Zeit $t$ und
+unabh\"angig von vorherigen Ereignissen ist.
+Die Wahrscheinlichkeit $P$ eines Ereignisses innerhalb eines Bins der Breite $\Delta t$ ist
+\[ P = \lambda \cdot \Delta t \]
+f\"ur gen\"ugend kleine $\Delta t$.
+\begin{parts}
+  
+  \part Schreibe eine Funktion die $n$ homogene Poisson Spiketrains
+  einer gegebenen Dauer $T_{max}$ mit Rate $\lambda$ erzeugt.
+  \begin{solution}
+    \lstinputlisting{hompoissonspikes.m}
+  \end{solution}
+  
+  \part Benutze diese Funktion um einige Trials von Spikes zu erzeugen
+  und plotte diese als Spikeraster.
+  \begin{solution}
+    \begin{lstlisting}
+      spikes = hompoissonspikes( 10, 100.0, 0.5 );
+      spikeraster( spikes )
+    \end{lstlisting}
+    \mbox{}\\[-3ex]
+    \colorbox{white}{\includegraphics[width=0.7\textwidth]{poissonraster100hz}}
+  \end{solution}
+  
+  \part Berechne Histogramme aus den Interspikeintervallen von $n$
+  Poisson Spiketrains mit der Rate $\lambda=100$\,Hz. Ver\"andere
+  \"uber die Dauer $T_{max}$ der Spiketrains und die Anzahl $n$ der
+  Trials die Anzahl der Intervalle und ver\"andere auch die Binbreite
+  des Histograms (fang mit 1\,ms an). Wieviele Interspikeintervalle
+  werden ben\"otigt, um ein ``sch\"ones'' Histogramm zu erhalten? Wie
+  lange m\"usste man also von dem Neuron ableiten?
+  \begin{solution}
+    About 5000 intervals for 25 bins. This corresponds to a $5000 /
+    100\,\hertz = 50\,\second$ recording of a neuron firing with
+    100\,\hertz.
+  \end{solution}
+  
+  \part Vergleiche Interspike-Intervall Histogramme von Poisson-Spikes
+  verschiedener Raten $\lambda$ mit der theoretisch zu erwartenden Verteilung
+  der Intervalle $T$ des Poisson Prozesses
+  \[ p(T) = \lambda e^{-\lambda T} \; .\]
+  \begin{solution}
+    \lstinputlisting{hompoissonisih.m}
+    \colorbox{white}{\includegraphics[width=0.48\textwidth]{poissonisih100hz}}
+    \colorbox{white}{\includegraphics[width=0.48\textwidth]{poissonisih20hz}}
+  \end{solution}
+  
+  \part \extra Was pasiert mit den Histogrammen, wenn die Binbreite
+  der Histogramme kleiner als das bei der Erzeugung der Poisson
+  Spiketrains verwendete $\Delta t$ ist?
+  \begin{solution}
+    Die Bins zwischen der durch $\Delta t$ vorgegebenen
+    Diskretisierung haben den Wert 0. Dadurch werden aber die anderen
+    durch die Normierung h\"oher als sie sein sollten.
+  \end{solution}
+  
+  \part Plotte den Mittelwert der Interspikeintervalle, die
+  dazugeh\"orige Standardabweichung und den Variationskoeffizienten
+  als Funktion der Rate $\lambda$ des Poisson Prozesses. Vergleiche
+  die Ergebnisse mit den theoretischen Erwartungen (siehe Vorlesungsskript).
+  \begin{solution}
+    \lstinputlisting{hompoissonisistats.m}
+    \colorbox{white}{\includegraphics[width=0.98\textwidth]{poissonisistats}}
+  \end{solution}
+  
+  \part Plotte die seriellen Korrelationen von Poisson-Spiketrains und
+  erkl\"are kurz das Ergebniss.
+  \begin{solution}
+    \mbox{}\\[-2ex]\hspace*{2cm}
+    \colorbox{white}{\includegraphics[width=0.8\textwidth]{poissonserial100hz}}\\
+    Alle Korrelationen zwischen Interspikeintervallen sind Null, da
+    beim Poisson Prozess das Auftreten jedes Spikes unabh\"angig von
+    den vorherigen Spikes ist.
+  \end{solution}
+  
+  \part Vergleiche Histogramme von Spikecounts gez\"ahlt in Fenstern
+  der Breite $W$ mit der Poisson Verteilung
+  \[ P(k) = \frac{(\lambda W)^ke^{\lambda W}}{k!} \; , \] 
+  wobei die Rate $\lambda$ aus den Daten bestimmt werden
+  soll. Hinweis: es gibt eine \code{matlab} Funktion, die die
+  Fakult\"at $n!$ berechnet.
+  \begin{solution}
+    \lstinputlisting{../code/counthist.m}
+    \colorbox{white}{\includegraphics[width=0.48\textwidth]{poissoncounthistdist100hz10ms}}
+    \colorbox{white}{\includegraphics[width=0.48\textwidth]{poissoncounthistdist100hz100ms}}
+  \end{solution}
+  
+  \part Schreibe eine Funktion, die die mittlere Anzahl, die Varianz
+  und den Fano-Faktor der Anzahl der Spikes in einem Fenster der
+  Breite $W$ bestimmt. Benutze die Funktion, um diese Parameter f\"ur
+  verschiedene Fensterbreiten $W$ zu bestimmen. Zwei Plots sollen aus
+  den Ergebnissen angefertigt werden: (i) Varianz gegen Mittelwert der counts.
+  (ii) Fano Faktor als Funktion der Fensterbreite.
+  \begin{solution}
+    \lstinputlisting{../code/fano.m}
+    \colorbox{white}{\includegraphics[width=0.98\textwidth]{poissonfano100hz}}
+  \end{solution}
+  
+\end{parts}
+
+\end{questions}
+
+\end{document}

pointprocesses/lecture/pointprocesses-chapter.tex

@@ -224,48 +224,3 @@
 
 \end{document}
 
-
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-\section{\tr{Homogeneous Poisson process}{Homogener Poisson Prozess}}
-
-\begin{figure}[t]
-  \includegraphics[width=1\textwidth]{poissonraster100hz}
-  \caption{\label{hompoissonfig}Rasterplot von Poisson-Spikes.}
-\end{figure}
-
-The probability $p(t)\delta t$ of an event occuring at time $t$
-is independent of $t$ and independent of any previous event
-(independent of event history).
-
-The probability $P$ for an event occuring within a time bin of width $\Delta t$
-is
-\[ P=\lambda \cdot \Delta t \]
-for a Poisson process with rate $\lambda$.
-
-\subsection{Statistics of homogeneous Poisson process}
-
-\begin{figure}[t]
-  \includegraphics[width=0.45\textwidth]{poissonisihexp20hz}\hfill
-  \includegraphics[width=0.45\textwidth]{poissonisihexp100hz}
-  \caption{\label{hompoissonisihfig}Interspike interval histograms of poisson spike train.}
-\end{figure}
-
-\begin{itemize}
-\item Exponential distribution of intervals $T$: $p(T) = \lambda e^{-\lambda T}$
-\item Mean interval $\mu_{ISI} = \frac{1}{\lambda}$
-\item Variance of intervals $\sigma_{ISI}^2 = \frac{1}{\lambda^2}$
-\item Coefficient of variation $CV_{ISI} = 1$
-\item Serial correlation $\rho_k =0$ for $k>0$ (renewal process!)   
-\item Fano factor $F=1$
-\end{itemize}
-
-\subsection{Count statistics of Poisson process}
-
-\begin{figure}[t]
-  \includegraphics[width=0.48\textwidth]{poissoncounthistdist100hz10ms}\hfill
-  \includegraphics[width=0.48\textwidth]{poissoncounthistdist100hz100ms}
-  \caption{\label{hompoissoncountfig}Count statistics of poisson spike train.}
-\end{figure}
-
-Poisson distribution:
-\[ P(k) = \frac{(\lambda W)^ke^{\lambda W}}{k!} \]

pointprocesses/lecture/pointprocesses.tex

@@ -56,29 +56,38 @@ erzeugt. Zum Beispiel:
 \item Diffusions Koeffizient $D_{ISI} = \frac{\sigma_{ISI}^2}{2\mu_{ISI}^3}$.
 \end{itemize}
 
-\subsection{Interval return maps}
+\subsection{Korrelationen der Intervalle}
-Scatter plot von aufeinander folgenden Intervallen $(T_{i+k}, T_i)$ getrennt durch das ``lag'' $k$.
+In ``return maps'' werden die um das ``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}
-  \caption{\label{returnmapfig}Interspike-Intervall return maps and serial correlations.}
+  \caption{\label{returnmapfig}Interspike-Intervall return maps und
+    serielle Korrelationen zwischen aufeinander folgenden Intervallen
+    im Abstand des Lags $k$.}
 \end{figure}
 
-\subsection{Serielle Korrelationen der Intervalle}
+Solche Ab\"angigkeiten werden durch die serielle Korrelation der
-Korrelationskoeffizient zwischen aufeinander folgenden Intervallen getrennt durch ``lag'' $k$:
+Intervalle quantifiziert.  Das ist der Korrelationskoeffizient
-\[ \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)} \]
+zwischen aufeinander folgenden Intervallen getrennt durch ``Lag'' $k$:
-$\rho_0=1$ (Korrelation jedes Intervalls mit sich selber).
+\[ \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).
 
 
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 \section{Z\"ahlstatistik}
 
-\begin{figure}[t]
+% \begin{figure}[t]
-  \includegraphics[width=0.48\textwidth]{poissoncounthist100hz10ms}\hfill
+%   \includegraphics[width=0.48\textwidth]{poissoncounthist100hz10ms}\hfill
-  \includegraphics[width=0.48\textwidth]{poissoncounthist100hz100ms}
+%   \includegraphics[width=0.48\textwidth]{poissoncounthist100hz100ms}
-  \caption{\label{countstatsfig}Count Statistik.}
+%   \caption{\label{countstatsfig}Count Statistik.}
-\end{figure}
+% \end{figure}
 
 Statistik der Anzahl der Ereignisse $N_i$ innerhalb von Beobachtungsfenstern $i$ der Breite $W$.
 \begin{itemize}
@@ -87,21 +96,69 @@ Statistik der Anzahl der Ereignisse $N_i$ innerhalb von Beobachtungsfenstern $i$
 \item Varianz der Anzahl: $\sigma_N^2 = \langle (N - \langle N \rangle)^2 \rangle$.
 \item 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'', Feuerrate) gemessen in Hertz
 \[ r = \frac{\langle N \rangle}{W} \; . \]
 
+% \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}
+%   \caption{\label{fanofig}Fano factor.}
+% \end{figure}
+
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+\section{Homogener Poisson Prozess}
+F\"ur kontinuierliche Me{\ss}gr\"o{\ss}en ist die Normalverteilung
+u.a. wegem dem Zentralen Grenzwertsatz die Standardverteilung. Eine
+\"ahnliche Rolle spilet bei Punktprozessen der ``Poisson Prozess''.
+
+Beim 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. Die
+Wahrscheinlichkeit zu irgendeiner Zeit ein Ereigniss in einem kleinen
+Zeitfenster der Breite $\Delta t$ zu bekommen ist
+\[ P = \lambda \cdot \Delta t \; . \]
+
 \begin{figure}[t]
-  \begin{minipage}[t]{0.49\textwidth}
+  \includegraphics[width=1\textwidth]{poissonraster100hz}
-    Poisson process $\lambda=100$\,Hz:\\
+  \caption{\label{hompoissonfig}Rasterplot von Spikes eine homogenen
-    \includegraphics[width=1\textwidth]{poissonfano100hz}
+    Poisson Prozesse mit $\lambda=100$\,Hz.}
-  \end{minipage}
-  \hfill
-  \begin{minipage}[t]{0.49\textwidth}
-    LIF $I=10$, $\tau_{adapt}=100$\,ms:\\
-    \includegraphics[width=1\textwidth]{lifadaptfano10-100ms}
-  \end{minipage}
-  \caption{\label{fanofig}Fano factor.}
 \end{figure}
 
+Beim inhomogenen Poisson Prozess h\"angt die Rate $\lambda$ von der
+Zeit ab: $\lambda = \lambda(t)$.
 
+\begin{figure}[t]
+  \includegraphics[width=0.45\textwidth]{poissonisihexp20hz}\hfill
+  \includegraphics[width=0.45\textwidth]{poissonisihexp100hz}
+  \caption{\label{hompoissonisihfig}Interspikeintervallverteilungen
+    zweier Poissonprozesse.}
+\end{figure}
+
+Der homogne Poissonprozess hat folgende Eigenschaften:
+\begin{itemize}
+\item Die Intervalle $T$ sind exponentiell verteilt: $p(T) = \lambda e^{-\lambda T}$ .
+\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 seriellen Korrelationen $\rho_k =0$ for $k>0$, da das
+  Auftreten der Ereignisse unabh\"angig von der Vorgeschichte ist. Ein
+  solcher Prozess wird auch Erneuerungsprozess genannt (``renewal
+  process'').
+\item Die Anzahl der Ereignisse $k$ innerhalb eines Fensters der L\"ange W ist Poissonverteilt:
+\[ P(k) = \frac{(\lambda W)^ke^{\lambda W}}{k!} \]
+\item Der Fano Faktor ist immer $F=1$ .
+\end{itemize}
+
+\begin{figure}[t]
+  \includegraphics[width=0.48\textwidth]{poissoncounthistdist100hz10ms}\hfill
+  \includegraphics[width=0.48\textwidth]{poissoncounthistdist100hz100ms}
+  \caption{\label{hompoissoncountfig}Z\"ahlstatistik von Poisson Spikes.}
+\end{figure}

programming/code/firing_rates.py

@@ -0,0 +1,198 @@
+import numpy as np
+import matplotlib.pyplot as plt
+import scipy.io as spio
+import scipy as sp
+import seaborn as sb
+from IPython import embed
+sb.set_context("paper")
+
+
+def set_axis_fontsize(axis, label_size, tick_label_size=None, legend_size=None):
+    """
+        Sets axis, tick label and legend font sizes to the desired size.
+
+    :param axis: the axes object
+    :param label_size: the size of the axis label
+    :param tick_label_size: the size of the tick labels. If None, lable_size is used
+    :param legend_size: the size of the font used in the legend.If None, the label_size is used.
+
+    """
+    if not tick_label_size:
+        tick_label_size = label_size
+    if not legend_size:
+        legend_size = label_size
+    axis.xaxis.get_label().set_fontsize(label_size)
+    axis.yaxis.get_label().set_fontsize(label_size)
+    for tick in axis.xaxis.get_major_ticks() + axis.yaxis.get_major_ticks():
+        tick.label.set_fontsize(tick_label_size)
+
+    l = axis.get_legend()
+    if l:
+        for t in l.get_texts():
+            t.set_fontsize(legend_size)
+
+
+def get_instantaneous_rate(times, max_t=30., dt=1e-4):
+    time = np.arange(0., max_t, dt)
+    indices = np.asarray(times / dt, dtype=int)
+    intervals = np.diff(np.hstack(([0], times)))
+    inst_rate = np.zeros(time.shape)
+    
+    for i, index in enumerate(indices[1:]):
+        inst_rate[indices[i-1]:indices[i]] = 1/intervals[i]
+    return time, inst_rate
+
+
+def plot_isi_rate(spike_times, max_t=30, dt=1e-4):
+    times = np.squeeze(spike_times[0][0])[:50000]
+    time, rate = get_instantaneous_rate(times, max_t=50000*dt)
+    
+    rates = np.zeros((len(rate), len(spike_times)))
+    for i in range(len(spike_times)):
+        _, rates[:, i] = get_instantaneous_rate(np.squeeze(spike_times[i][0])[:50000],
+                                                max_t=50000*dt)
+    avg_rate = np.mean(rates, axis=1)
+    rate_std = np.std(rates, axis=1)
+    
+    fig = plt.figure()
+    ax1 = fig.add_subplot(311)
+    ax2 = fig.add_subplot(312)
+    ax3 = fig.add_subplot(313)
+    
+    ax1.vlines(times[times < (50000*dt)], ymin=0, ymax=1, color="dodgerblue", lw=1.5)
+    ax1.set_ylabel("skpikes", fontsize=12)
+    set_axis_fontsize(ax1, 12)
+
+    ax2.plot(time, rate, label="instantaneous rate, trial 1")
+    ax2.set_ylabel("firing rate [Hz]", fontsize=12)
+    ax2.legend(fontsize=12)
+    set_axis_fontsize(ax2, 12)
+
+    ax3.fill_between(time, avg_rate+rate_std, avg_rate-rate_std, color="dodgerblue",
+                     alpha=0.5, label="standard deviation")
+    ax3.plot(time, avg_rate, label="average rate")
+    ax3.set_xlabel("times [s]", fontsize=12)
+    ax3.set_ylabel("firing rate [Hz]", fontsize=12)
+    ax3.legend(fontsize=12)
+    ax3.set_ylim([0, 450])
+    set_axis_fontsize(ax3, 12)
+    
+    fig.set_size_inches(15, 10)
+    fig.subplots_adjust(left=0.1, bottom=0.125, top=0.95, right=0.95)
+    fig.set_facecolor("white")
+    fig.savefig("../lectures/images/instantaneous_rate.png")
+    plt.close()
+
+
+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
+    hist, _ = sp.histogram(times, bins)
+    rate = np.zeros(time.shape)
+    
+    for i, b in enumerate(bin_indices[1:]):
+        rate[bin_indices[i-1]:b] = hist[i-1]/bin_width
+    return time, rate
+
+
+def plot_bin_rate(spike_times, bin_width, max_t=30, dt=1e-4):
+    times = np.squeeze(spike_times[0][0])
+    time, rate = get_binned_rate(times)
+    rates = np.zeros((len(rate), len(spike_times)))
+    for i in range(len(spike_times)):
+        _, rates[:, i] = get_binned_rate(np.squeeze(spike_times[i][0]))
+    avg_rate = np.mean(rates, axis=1)
+    rate_std = np.std(rates, axis=1)
+    
+    fig = plt.figure()
+    ax1 = fig.add_subplot(311)
+    ax2 = fig.add_subplot(312)
+    ax3 = fig.add_subplot(313)
+
+    ax1.vlines(times[times < (100000*dt)], ymin=0, ymax=1, color="dodgerblue", lw=1.5)
+    ax1.set_ylabel("skpikes", fontsize=12)
+    ax1.set_xlim([0, 5])
+    set_axis_fontsize(ax1, 12)
+
+    ax2.plot(time, rate, label="binned rate, trial 1")
+    ax2.set_ylabel("firing rate [Hz]", fontsize=12)
+    ax2.legend(fontsize=12)
+    ax2.set_xlim([0, 5])
+    set_axis_fontsize(ax2, 12)
+
+    ax3.fill_between(time, avg_rate+rate_std, avg_rate-rate_std, color="dodgerblue",
+                     alpha=0.5, label="standard deviation")
+    ax3.plot(time, avg_rate, label="average rate")
+    ax3.set_xlabel("times [s]", fontsize=12)
+    ax3.set_ylabel("firing rate [Hz]", fontsize=12)
+    ax3.legend(fontsize=12)
+    ax3.set_xlim([0, 5])
+    ax3.set_ylim([0, 450])
+    set_axis_fontsize(ax3, 12)
+
+    fig.set_size_inches(15, 10)
+    fig.subplots_adjust(left=0.1, bottom=0.125, top=0.95, right=0.95)
+    fig.set_facecolor("white")
+    fig.savefig("../lectures/images/binned_rate.png")
+    plt.close()
+  
+
+def get_convolved_rate(times, sigma, max_t=30., dt=1.e-4):
+    time = np.arange(0., max_t, dt)
+    kernel = sp.stats.norm.pdf(np.arange(-8*sigma, 8*sigma, dt),loc=0,scale=sigma)
+    indices = np.asarray(times/dt, dtype=int)
+    rate = np.zeros(time.shape)
+    rate[indices] = 1.;
+    conv_rate = np.convolve(rate, kernel, mode="same")
+    return time, conv_rate
+
+
+def plot_conv_rate(spike_times, sigma=0.05, max_t=30, dt=1e-4):
+    times = np.squeeze(spike_times[0][0])
+    time, rate = get_convolved_rate(times, sigma)
+   
+    rates = np.zeros((len(rate), len(spike_times)))
+    for i in range(len(spike_times)):
+        _, rates[:, i] = get_convolved_rate(np.squeeze(spike_times[i][0]), sigma)
+    avg_rate = np.mean(rates, axis=1)
+    rate_std = np.std(rates, axis=1)
+    
+    fig = plt.figure()
+    ax1 = fig.add_subplot(311)
+    ax2 = fig.add_subplot(312)
+    ax3 = fig.add_subplot(313)
+
+    ax1.vlines(times[times < (100000*dt)], ymin=0, ymax=1, color="dodgerblue", lw=1.5)
+    ax1.set_ylabel("skpikes", fontsize=12)
+    ax1.set_xlim([0, 5])
+    set_axis_fontsize(ax1, 12)
+
+    ax2.plot(time, rate, label="convolved rate, trial 1")
+    ax2.set_ylabel("firing rate [Hz]", fontsize=12)
+    ax2.legend(fontsize=12)
+    ax2.set_xlim([0, 5])
+    set_axis_fontsize(ax2, 12)
+
+    ax3.fill_between(time, avg_rate+rate_std, avg_rate-rate_std, color="dodgerblue",
+                     alpha=0.5, label="standard deviation")
+    ax3.plot(time, avg_rate, label="average rate")
+    ax3.set_xlabel("times [s]", fontsize=12)
+    ax3.set_ylabel("firing rate [Hz]", fontsize=12)
+    ax3.legend(fontsize=12)
+    ax3.set_xlim([0, 5])
+    ax3.set_ylim([0, 450])
+    set_axis_fontsize(ax3, 12)
+
+    fig.set_size_inches(15, 10)
+    fig.subplots_adjust(left=0.1, bottom=0.125, top=0.95, right=0.95)
+    fig.set_facecolor("white")
+    fig.savefig("../lectures/images/convolved_rate.png")
+    plt.close()
+
+
+if __name__ == "__main__":
+    spike_times =  spio.loadmat('lifoustim.mat')["spikes"]
+    plot_isi_rate(spike_times)
+    plot_bin_rate(spike_times, 0.05)
+    plot_conv_rate(spike_times, 0.025)

programming/exercises/plotting_psth.tex

@@ -0,0 +1,78 @@
+\documentclass[12pt,a4paper,pdftex]{exam}
+
+\usepackage[german]{babel}
+\usepackage{natbib}
+\usepackage{graphicx}
+\usepackage[small]{caption}
+\usepackage{sidecap}
+\usepackage{pslatex}
+\usepackage{amsmath}
+\usepackage{amssymb}
+\setlength{\marginparwidth}{2cm}
+\usepackage[breaklinks=true,bookmarks=true,bookmarksopen=true,pdfpagemode=UseNone,pdfstartview=FitH,colorlinks=true,citecolor=blue]{hyperref}
+
+%%%%% text size %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+\usepackage[left=20mm,right=20mm,top=25mm,bottom=25mm]{geometry}
+\pagestyle{headandfoot} \header{{\bfseries\large \"Ubung
+    }}{{\bfseries\large Peri Stimulus Time Histogram}}{{\bfseries\large 28. Oktober, 2015}}
+\firstpagefooter{Dr. Jan Grewe}{Phone: 29 74588}{Email:
+  jan.grewe@uni-tuebingen.de} \runningfooter{}{\thepage}{}
+
+\setlength{\baselineskip}{15pt}
+\setlength{\parindent}{0.0cm}
+\setlength{\parskip}{0.3cm}
+\renewcommand{\baselinestretch}{1.15}
+
+\newcommand{\code}[1]{\texttt{#1}}
+\renewcommand{\solutiontitle}{\noindent\textbf{L\"osung:}\par\noindent}
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
+\begin{document}
+
+\vspace*{-6.5ex}
+\begin{center}
+  \textbf{\Large Einf\"uhrung in die wissenschaftliche Datenverarbeitung}\\[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} \\
+\end{center}
+
+\begin{questions}
+  \question Graphische Darstellung der zeitabh\"angigen Antwort eines
+  Neurons. PSTH auf Basis der instantanen Feuerrate. Verwendet den Datensatz \code{}
+  \begin{parts}
+    \part Schreibt eine Funktion, die einen Vektor mit Spikezeiten,
+    die Dauer des Trials, und die zeitliche Aufl\"osung entgegennimmt
+    und die Zeitachse sowie die Feuerrate zur\"uckgiebt.
+    \part Benutzt diese Funktion in einem Skript und stellt das PSTH
+    eines einzelnen Trials sowie den Mittelwert \"uber alle Trials
+    dar.
+    \part Erweitert das Programm so, dass die Abbildung den Standards
+    z.B. vom \textit{Journal of Neuroscience} entspricht
+    (Schriftgr\"o{\ss}e, Abbildungsgr\"o{\ss}e).
+    \part Die Abbildung sollte als pdf gespeichert werden.
+  \end{parts}
+
+  \question Wie zuvor nur unter Verwendung der Binning Methode.
+
+  \question Wie zuvor nur unter Verwendung der Faltungsmethode.
+
+  \question Entscheidet euch f\"ur eine Varainte und erweitert das
+  entsprechende Skript, sodass auch die Interspikeintervallverteilung
+  und die Verteilung der Spikecounts dargestellt werden. Die
+  Abbildungen sollten nat\"urlich ``publikationsreif'' sein und
+  gespeichert werden.
+  
+  \question Einige trials sind anders als die \"Ubrigen. Benutzt den
+  Rasterplot um sie zu finden. Speichert alle Abbildungen in
+  ``publikationsreifer'' Form.
+  \begin{parts}
+    \part Benutzt den Rasterplot um sie zu finden.
+    \part Plottet die Verteilung der Spike counts.
+    \part Filtert all die trials heraus, deren spike count mehr als
+    $2\sigma$ vom Mittelwert abweicht.
+    \part Plottet das PSTH vor und nach dem Filtern.
+  \end{parts}
+
+\end{questions}
+
+\end{document}

programming/lectures/plotting_spike_trains.tex

@@ -176,14 +176,14 @@
   }
   \only <4> {
     \begin{figure}
-      \includegraphics[width=0.4\columnwidth]{images/badbarright}
+      \includegraphics[width=0.3\columnwidth]{images/badbarright}
     \end{figure}
     \vspace{0.5cm}
     \url{https://en.wikipedia.org/wiki/Misleading_graph}
   }
   \only <5> {
     \begin{figure}
-      \includegraphics[width=0.4\columnwidth]{images/badbarleft}
+      \includegraphics[width=0.3\columnwidth]{images/badbarleft}
     \end{figure}
     \vspace{0.5cm}
     \url{https://en.wikipedia.org/wiki/Misleading_graph}
@@ -352,19 +352,18 @@ saveas(fig, 'spike_detection.pdf', 'pdf')
  \item Deutliche Unterscheidbarkeit von Kurven.
  \item Keine suggestive Darstellung. 
  \item Ausgewogenheit von Linienst\"arken Schrift- und Plotgr\"o{\ss}e.
- \item Vermeidung von Suggestiven Darstellungen.
  \item Fehlerbalken, wenn sie angebracht sind.
  \end{enumerate}
 \end{frame}
 
 
 \begin{frame}[plain]
-\huge{2. Spiketrain Analyse I}
+\huge{2. Spiketrain Analyse}
 \end{frame}
 
 
 \begin{frame}
-  \frametitle{Spiketrain Analyse I}
+  \frametitle{Spiketrain Analyse}
   \framesubtitle{Rasterplot}
   \begin{figure}
     \centering
@@ -374,92 +373,221 @@ saveas(fig, 'spike_detection.pdf', 'pdf')
 
 
 \begin{frame}
-  \frametitle{Spiketrain Analyse I}
+  \frametitle{Spiketrain Analyse}
   \framesubtitle{Rasterplot}
   \"Ubung:
   \begin{enumerate}
-  \item Ladet die Datei: electrorecptor\_spike\_times.mat aus dem
+  \item Ladet die Datei: \code{lifoustim.mat} aus dem
     Ilias Ordner.
   \item Der Datensatz enth\"alt die Zeiten von Aktionspotentialen.
-  \item Schaut euch den Inhalt und skizziert wie das Problem gel\"ost
+  \item Erzeugt einen sch\"onen Rasterplot der Zellantworten, speichert ihn.
-    werden k\"onnte.
+  \item Welche Information liefert er, welche Information ist schwer
-  \item Erzeugt einen Rasterplot.
+    abzulesen?
+  \end{enumerate}
+\end{frame}
+
+\begin{frame}
+  \frametitle{Spiketrain Analyse}
+  \framesubtitle{Zeitabh\"angige Feuerrate, PSTH} 
+  
+  Darstellung der Feurreate eines Neuron als Funktion der Zeit. Es
+  gibt verschiedene Methoden dieses \textbf{P}eri \textbf{S}timulus
+  \textbf{T}ime \textbf{H}istogram zu erstellen.
+  \begin{enumerate}
+  \item Auf Basis der \textit{instantanen} Feuerrate.
+  \item Auf Basis des Zeithistogramms.
+  \item Durch Faltung der Zellantwort mit einem Gauss Kern.
   \end{enumerate}
 \end{frame}
 
 
 \begin{frame}
-  \frametitle{Spiketrain Analyse I}
+  \frametitle{Spiketrain Analyse}
-  \framesubtitle{Feuerrate \"uber
+  \framesubtitle{Zeitabh\"angige Feuerrate, PSTH --- Instantane Feuerrate ---}
-    die Zeit}
+  \Large{Berechnung des PSTHs durch die Instantane Feurerrate:}
+  \normalsize
   \begin{enumerate}
   \item Die Feuerrate kann aus dem Abstand zwischen zwei
     aufeinanderfolgenden Aktionspotentialen
     (\textbf{Interspikeinterval}) berechnet werden. \pause
-  \item Auch als \textbf{Instantane Feuerrate} bezeichnet.
+  \item Die \textbf{Instantane Feuerrate} wird aus dem Kehrwert des
+    Interspikeintervals berechnet.
   \end{enumerate}
+\end{frame}
+
+
+\begin{frame}
+  \frametitle{Spiketrain Analyse}
+  \framesubtitle{Zeitabh\"angige Feuerrate, PSTH --- Instantane Feuerrate ---}
   \begin{figure}
     \centering
-    \includegraphics[width=0.5\columnwidth]{images/isi}
+    \includegraphics[width=0.9\columnwidth]{images/instantaneous_rate}
   \end{figure}
 \end{frame}
 
 
 \begin{frame}
-  \frametitle{Spiketrain Analyse I}
+  \frametitle{Spiketrain Analyse}
-  \framesubtitle{Feuerrate \"uber
+  \framesubtitle{Zeitabh\"angige Feuerrate, PSTH --- Instantane Feuerrate ---}
-    die Zeit}
+  Berechnung des PSTHs durch die Instantane Feurerrate:
+  \textbf{Vorteile:}
+
   \begin{enumerate}
-  \item Z.B. das \textbf{P}eri- \textbf{S}timulus - \textbf{T}ime -
+  \item Sehr einfach zu Berechnen.
-    \textbf{H}istogram, \textbf{PSTH}
+  \item Macht keine Annahmen \"uber ein Zeitraster, oder die Zeitskala
-  \item Wird berechnet indem die Zeit in ``bins'' geteilt wird die
+    der neuronalen Verarbeitung.
-    Anzahl Spikes pro bin gezaehlt werden.
+  \end{enumerate}
-  \item Die Anzahl wird dann in eine Feuerrate umgerechnet.
+
-  \end{enumerate}\pause
+  \textbf{Nachteile:}
+  \begin{enumerate}
+  \item Die Feuerrate ist nie null, auch wenn f\"ur lange Zeit kein
+    Aktionspotential auftritt.
+  \item Verh\"alt sich im Fourrier Raum nicht sehr sch\"on.
+  \end{enumerate}
+\end{frame}
+
+
+\begin{frame}
+  \frametitle{Spiketrain Analyse}
+  \framesubtitle{Zeitabh\"angige Feuerrate, PSTH --- Binning Methode ---}
+  \Large{Binning Methode:}
+  \normalsize
+  \begin{enumerate}
+  \item Die Zeitachse wird in gleich gro{\ss}e Abschnitte ``bins''
+    unterteilt.
+  \item F\"ur jedes ``bin'' wird die Anzahl vorkommender
+    Aktionspotentiale gez\"ahlt.
+  \item Der Spike-count pro bin muss nun noch in die Rate umgerechnet
+    werden.
+  \end{enumerate}
+\end{frame}
+
+
+\begin{frame}
+  \frametitle{Spiketrain Analyse}
+  \framesubtitle{Zeitabh\"angige Feuerrate, PSTH --- Binning Methode ---}
   \begin{figure}
-    \centering
+    \includegraphics[width=0.9\columnwidth]{images/binned_rate}
-    \includegraphics[width=0.5\columnwidth]{images/binning}
   \end{figure}
 \end{frame}
 
 
+\begin{frame}
+  \frametitle{Spiketrain Analyse}
+  \framesubtitle{Zeitabh\"angige Feuerrate, PSTH --- Binning Methode ---}
+
+  Berechnung des PSTHs durch die Binning Methode:
+
+  \textbf{Vorteile:}
+  \begin{enumerate}
+  \item Sehr einfach zu Berechnen.
+  \item Zeigt nur da Aktivit\"at an, wo auch Aktionspotentiale
+    generiert wurden.
+  \end{enumerate}
+
+  \textbf{Nachteile:}
+  \begin{enumerate}
+  \item Mach Annahmen \"uber die relevante Zeitskala neuronaler
+    Verarbeitung.
+  \item Die Zeitachse wird diskretisiert.
+  \item Verh\"alt sich im Fourrier Raum nicht sehr sch\"on.
+  \end{enumerate}
+\end{frame}
+
+
+\begin{frame}[fragile]
+  \frametitle{Spiketrain Analyse}
+  \framesubtitle{Zeitabh\"angige Feuerrate, PSTH --- Faltungsmethode ---}
+  \Large{Faltung mit einem Gauss Kern:}
+  \normalsize
+  \[r_{est}(t) = \int_{-\infty}^{\infty}d\tau \omega(\tau)\rho(t-\tau) \]
+
+  wobei $\omega(\tau)$ der Gauss Kern und $\rho(t)$ die Antwortfunktion ist.
+  
+  Gl\"ucklicherweise m\"ussen wir das nicht selbst implementieren... 
+\end{frame}
+
+
 \begin{frame}[fragile]
-  \frametitle{Spiketrain Analyse I}
+  \frametitle{Spiketrain Analyse}
-  \framesubtitle{Feuerrate \"uber die Zeit}
+  \framesubtitle{Zeitabh\"angige Feuerrate, PSTH --- Faltungsmethode ---}
+  \large Algortihmus:
+  \normalsize
   \begin{enumerate}
-  \item Z.B. das \textbf{P}eri- \textbf{S}timulus - \textbf{T}ime -
+  \item Die neuronalen Antworten werden  ``bin\"ar'' ausgedr\"uckt.
-    \textbf{H}istogram, \textbf{PSTH}
+  \item Ein Filterkern wird berechnet, der das Integral 1 hat.
-  \item Wird berechnet indem man die Aktivit\"at bin\"ar ausdr\"uckt.
+  \item Mithilfe der Faltung (\code{conv} Funktion) wird jede 1 durch
-  \item Jede 1 wird dann duch einen ``Kern'' ersetzt.
+    den ``Kern'' ersetzt.\\
-  \item Der Vorgang heisst Faltung (\verb+conv+).
+    \code{conv(x, kern, 'mode', 'same')}
-  \end{enumerate}\pause
+  \end{enumerate}
+\end{frame}
+
+
+\begin{frame}
+  \frametitle{Spiketrain Analyse}
+  \framesubtitle{Feurerrate als Funktion der Zeit --- Faltungsmethode ---}
   \begin{figure}
-    \centering
+    \includegraphics[width=0.9\columnwidth]{images/convolved_rate}
-    \includegraphics[width=0.5\columnwidth]{images/conv}
   \end{figure}
 \end{frame}
 
 
-\begin{frame} [fragile]
+\begin{frame}
- \frametitle{Spiketrain Analyse I}
+  \frametitle{Spiketrain Analyse}
- \framesubtitle{\"Ubung}
+  \framesubtitle{Feurerrate als Funktion der Zeit --- Faltungsmethode ---}
- \begin{enumerate}
+  \textbf{Vorteile:}
- \item Berechnet die Feuerrate eines Neurons mit einer der drei Methoden.
+  \begin{enumerate}
- \item Die Abbildung soll f\"ur eine einspaltige Abbildung im
+  \item Sehr ``nat\"urliche'' erscheinende Darstellung.
-   \textit{Journal of Neuroscience} geeignet sein
+  \item Sehr gutes Verhalten im Fourrier Raum.
-   (\url{http://www.jneurosci.org/site/misc/ifa_illustrations.xhtml}).
+  \end{enumerate}
- \item Erzeugt/ver\"andert/erweitert das Programm zum plotten so, dass
+  
-   die Abbildung automatisch erstellt und gespeichert wird.
+  \textbf{Nachteile:}
- \item Speichert die Abbildung als pdf.
+  \begin{enumerate}
- \item Ladet den Stimulus aus dem Ilias Ordner und benutzt die
+  \item Relativ rechenintensiv.
-   \verb+subplot+ Funktion um den Stimulus zu der neuronalen
+  \item Macht Annahmen \"uber die Zeitskalen neuronaler Verarbeitung.
-   Aktivit\"at zu plotten.
+  \end{enumerate}
- \end{enumerate}
+\end{frame}
+
+\begin{frame}[plain]
+  \huge{3. Analyse der Beziehung zwischen Stimulus und Antwort}
+\end{frame}
+
+
+\begin{frame}[fragile]
+  \frametitle{Spiketrain Analyse}
+  \framesubtitle{Spike-Triggered-Average}
+  \begin{itemize}
+  \item[] Die Antworten darzustellen ist gut und sch\"on, aber was sagt es uns?\pause
+  \item[] Idealerweise wollen wir die Antworten in Beziehung zum
+    hervorrufenden Stimulus setzen.
+  \item[] Eine Methode ist der sogenannte \textbf{Spike-Triggered-Average} (STA).
+  \end{itemize}
 \end{frame}
 
+
 \begin{frame}[fragile]
   \frametitle{Spiketrain Analyse}
   \framesubtitle{Spike-Triggered-Average}
+  
+  Der STA stellt den (mittleren) Stimulus dar, der zu einem
+  Aktionspotential gef\"uhrt hat:
+
+  \begin{equation}
+    STA(\tau) = \frac{1}{\langle n \rangle} \left\langle  \displaystyle\sum_{i=1}^{n}{s(t_i - \tau)} \right\rangle
+  \end{equation}
+  Wobei: $\tau$ ist eine bestimmte Zeit relativ zur Zeit eines
+  Aktionspotentials, $t_i$ ist der Zeitpunkt eines APs, $s(t)$ ist der
+  Stimulus.\\
+
+  Leider m\"ussen wir das selbst implementieren...
+\end{frame}
+
+
+\begin{frame}[fragile]
+   \frametitle{Spiketrain Analyse}
+  \framesubtitle{Spike-Triggered-Average}
+  \large{Algorithmus:}
+  \normalsize
   \begin{enumerate}
   \item Der \textbf{STA} ist der mittlere Stimulus, der zu einem Aktionspotential f\"uhrt. 
   \item F\"ur jeden Spike wird ein entsprechender Abschnitt um die Zeit des Spikes herausgeschnitten.
@@ -471,4 +599,45 @@ saveas(fig, 'spike_detection.pdf', 'pdf')
   \end{figure}
 \end{frame}
 
+\begin{frame}[fragile] 
+  \frametitle{Spiketrain Analyse}
+  \framesubtitle{Spike-Triggered-Average}
+  \vspace{-1em}
+  \begin{figure}
+    \centering
+    \includegraphics[width=0.6\columnwidth]{images/sta}
+  \end{figure}
+  \pause
+  \vspace{-0.5em}
+  Welche Information liefert der \textbf{STA}?
+  \small
+  \begin{enumerate}
+  \item Gibt es eine Beziehung zwischen Stimulus und Antwort?\pause
+  \item Gibt es eine Verz\"ogerung zwischen Stimulus und Antwort? Wie
+    gro{\ss} ist diese?\pause
+  \item Wie weit h\"angt das Auftreten eines Aktionspotentials von der
+    Vergangenheit ab? \pause
+  \item Kann die Zelle in die Zukunft sehen?
+  \end{enumerate}
+\end{frame}
+
 \end{document}
+
+
+
+\begin{frame} [fragile]
+ \frametitle{Spiketrain Analyse}
+ \framesubtitle{\"Ubung}
+ \begin{enumerate}
+ \item Berechnet die Feuerrate eines Neurons mit einer der drei Methoden.
+ \item Die Abbildung soll f\"ur eine einspaltige Abbildung im
+   \textit{Journal of Neuroscience} geeignet sein
+   (\url{http://www.jneurosci.org/site/misc/ifa_illustrations.xhtml}).
+ \item Erzeugt/ver\"andert/erweitert das Programm zum plotten so, dass
+   die Abbildung automatisch erstellt und gespeichert wird.
+ \item Speichert die Abbildung als pdf.
+ \item Ladet den Stimulus aus dem Ilias Ordner und benutzt die
+   \verb+subplot+ Funktion um den Stimulus zu der neuronalen
+   Aktivit\"at zu plotten.
+ \end{enumerate}
+\end{frame}

scientificcomputing-script.tex

@@ -237,6 +237,9 @@
 
 \renewcommand{\codepath}{pointprocesses/code/}
 \renewcommand{\texinputpath}{pointprocesses/lecture/}
-%\include{pointprocesses/lecture/pointprocesses}
+\include{pointprocesses/lecture/pointprocesses}
+
+\renewcommand{\codepath}{designpattern/code/}
+\include{designpattern/lecture/designpattern}
 
 \end{document}

statistics/code/lsq_gradient_sigmoid.m

@@ -1,9 +0,0 @@
-function gradient = lsq_gradient_sigmoid(parameter, x, y)
-h = 1e-6;
-
-gradient = zeros(size(parameter));
-for i = 1:length(parameter)
-    parameter_h = parameter;
-    parameter_h(i) = parameter_h(i) + h;
-    gradient(i) = (lsq_sigmoid_error(parameter_h, x, y) - lsq_sigmoid_error(parameter, x, y)) / h;
-end

statistics/code/lsq_sigmoid_error.m

@@ -1,8 +0,0 @@
-function error = lsq_sigmoid_error(parameter, x, y)
-% p(1) the amplitude
-% p(2) the slope
-% p(3) the x-shift
-% p(4) the y-shift
-
-y_est = parameter(1)./(1+ exp(-parameter(2) .* (x - parameter(3)))) + parameter(4);
-error = mean((y_est - y).^2);

statistics/code/sigmoidal_gradient_descent.m

@@ -1,44 +0,0 @@
-
-
-%% fit the sigmoid
-
-clear 
-close all
-
-load('iv_curve.mat')
-
-figure()
-plot(voltage, current, 'o')
-xlabel('voltate [mV]')
-ylabel('current [pA]')
-
-% amplitude, slope, x-shift, y-shift
-%parameter = [10 0.25 -50, 2.5];
-parameter = [20 0.5 -50, 2.5];
-
-eps = 0.1;
-% do the descent
-gradient = [];
-steps = 0;
-error = [];
-
-while isempty(gradient) || norm(gradient) > 0.01
-    steps = steps + 1; 
-    gradient = lsq_gradient_sigmoid(parameter, voltage, current);
-    error(steps) = lsq_sigmoid_error(parameter, voltage, current);
-    parameter = parameter - eps .* gradient; 
-end
-plot(1:steps, error)
-
-disp('gradient descent done!')
-disp(strcat('final position: ', num2str(parameter)))
-disp(strcat('final error: ', num2str(error(end))))
-
-%% use fminsearch
-parameter = [10 0.5 -50, 2.5];
-
-objective_function = @(p)lsq_sigmoid_error(p, voltage, current);
-param = fminunc(objective_function, parameter);
-disp(param)
-param1 = fminsearch(objective_function, parameter);
-disp(param1)