diff --git a/designpattern/lecture/Makefile b/designpattern/lecture/Makefile new file mode 100644 index 0000000..d570705 --- /dev/null +++ b/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 + + diff --git a/designpattern/lecture/designpattern-chapter.tex b/designpattern/lecture/designpattern-chapter.tex new file mode 100644 index 0000000..08ff012 --- /dev/null +++ b/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} diff --git a/designpattern/lecture/designpattern.tex b/designpattern/lecture/designpattern.tex new file mode 100644 index 0000000..c581b9b --- /dev/null +++ b/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} diff --git a/pointprocesses/code/counthist.m b/pointprocesses/code/counthist.m index d64fb23..9abb225 100644 --- a/pointprocesses/code/counthist.m +++ b/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)' ); diff --git a/pointprocesses/code/isihist.m b/pointprocesses/code/isihist.m index a9722fd..afd69d4 100644 --- a/pointprocesses/code/isihist.m +++ b/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.5, 0.5, sprintf( 'std=%.1f ms', 1000.0*sdisi ), 'Units', 'normalized' ) - text( 0.5, 0.4, sprintf( 'CV=%.2f', sdisi/misi ), 'Units', 'normalized' ) + text( 0.95, 0.8, sprintf( 'mean=%.1f ms', 1000.0*misi ), 'Units', 'normalized', 'HorizontalAlignment', 'right' ) + text( 0.95, 0.7, sprintf( 'std=%.1f ms', 1000.0*sdisi ), 'Units', 'normalized', 'HorizontalAlignment', 'right' ) + 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 diff --git a/pointprocesses/code/lifadapt.mat b/pointprocesses/code/lifadapt.mat new file mode 100644 index 0000000..d0272e3 Binary files /dev/null and b/pointprocesses/code/lifadapt.mat differ diff --git a/pointprocesses/code/lifadaptspikes.m b/pointprocesses/code/lifadaptspikes.m index 2ef1874..a186f0a 100644 --- a/pointprocesses/code/lifadaptspikes.m +++ b/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 -% in case of a vector as a stimulus the time step +% tmaxdt: in case of a single value stimulus: the duration of a trial +% 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 diff --git a/pointprocesses/code/lifoustim.mat b/pointprocesses/code/lifoustim.mat new file mode 100644 index 0000000..583f7f4 Binary files /dev/null and b/pointprocesses/code/lifoustim.mat differ diff --git a/pointprocesses/code/lifspikes.m b/pointprocesses/code/lifspikes.m index cfa0f55..dd11408 100644 --- a/pointprocesses/code/lifspikes.m +++ b/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 -% in case of a vector as a stimulus the time step +% tmaxdt: in case of a single value stimulus: the duration of a trial +% in case of a vector as a stimulus: the time step % D: the strength of additive white noise tau = 0.01; diff --git a/pointprocesses/code/lifspikesoustim.m b/pointprocesses/code/lifspikesoustim.m new file mode 100644 index 0000000..fd25766 --- /dev/null +++ b/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 \ No newline at end of file diff --git a/pointprocesses/code/pifou.mat b/pointprocesses/code/pifou.mat new file mode 100644 index 0000000..21d343d Binary files /dev/null and b/pointprocesses/code/pifou.mat differ diff --git a/pointprocesses/code/pifouspikes.m b/pointprocesses/code/pifouspikes.m index b6516cc..b671bd7 100644 --- a/pointprocesses/code/pifouspikes.m +++ b/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 -% in case of a vector as a stimulus the time step +% tmaxdt: in case of a single value stimulus: the duration of a trial +% 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 diff --git a/pointprocesses/code/plotcounthist.m b/pointprocesses/code/plotcounthist.m new file mode 100644 index 0000000..7a0a51d --- /dev/null +++ b/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); \ No newline at end of file diff --git a/pointprocesses/code/plotisih.m b/pointprocesses/code/plotisih.m new file mode 100644 index 0000000..bc3e7c9 --- /dev/null +++ b/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); \ No newline at end of file diff --git a/pointprocesses/code/plotserialcorr.m b/pointprocesses/code/plotserialcorr.m new file mode 100644 index 0000000..05f3ea8 --- /dev/null +++ b/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); \ No newline at end of file diff --git a/pointprocesses/code/plotspikeraster.m b/pointprocesses/code/plotspikeraster.m new file mode 100644 index 0000000..b7607f7 --- /dev/null +++ b/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); \ No newline at end of file diff --git a/pointprocesses/code/poisson.mat b/pointprocesses/code/poisson.mat new file mode 100644 index 0000000..7c2c577 Binary files /dev/null and b/pointprocesses/code/poisson.mat differ diff --git a/pointprocesses/code/poissonspikes.m b/pointprocesses/code/poissonspikes.m index f00291c..6a75f40 100644 --- a/pointprocesses/code/poissonspikes.m +++ b/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 - spikes{k} = find( x < p ) * dt; + x = rand(round(tmax/dt), 1); % uniform random numbers for each bin + spikes{k} = find(x < p) * dt; end end diff --git a/pointprocesses/code/psth.m b/pointprocesses/code/psth.m new file mode 100644 index 0000000..87da768 --- /dev/null +++ b/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 diff --git a/pointprocesses/code/spikeraster.m b/pointprocesses/code/spikeraster.m index f7ff5d4..7d4c12f 100644 --- a/pointprocesses/code/spikeraster.m +++ b/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= 1.0-p ) * dt; + spikes{k} = find( x(k,:) < p ) * dt; end end diff --git a/pointprocesses/exercises/isihist.pdf b/pointprocesses/exercises/isihist.pdf new file mode 100644 index 0000000..1fdab43 Binary files /dev/null and b/pointprocesses/exercises/isihist.pdf differ diff --git a/pointprocesses/exercises/pointprocesses01.tex b/pointprocesses/exercises/pointprocesses01.tex index b4d927c..e537be7 100644 --- a/pointprocesses/exercises/pointprocesses01.tex +++ b/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} -We use the Poisson process to generate spike trains on which we can test and imrpove some -standard analysis functions. - -A homogeneous Poisson process of rate $\lambda$ (measured in Hertz) is a point process -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 -\[ P = \lambda \cdot \Delta t \] -for sufficiently small $\Delta t$. -\begin{parts} + \question \qt{Statistik von Spiketrains} + In Ilias findet ihr die Dateien \code{poisson.mat}, + \code{pifou.mat}, und \code{lifadapt.mat}. Jede dieser Dateien + enth\"alt mehrere Trials von Spiketrains von einer bestimmten Art + von Neuron. Die Spikezeiten sind in Sekunden gemessen. + + Mit den folgenden Aufgaben wollen wir die Statistik der Spiketrains + der drei Neurone miteinander vergleichen. + \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}$ - with rate $\lambda$. - \begin{solution} - \lstinputlisting{hompoissonspikes.m} - \end{solution} + \part Schreibe eine Funktion, die einen einzigen Vektor mit den + Interspikeintervallen aller Trials von Spikezeiten zur\"uckgibt. + \begin{solution} + \lstinputlisting{../code/isis.m} + \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. - \begin{solution} - \lstinputlisting{../code/spikeraster.m} - \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 Schreibe eine Funktion, die die seriellen Korrelationen der + Interspikeintervalle f\"ur Lags bis zu \code{maxlag} berechnet und + plottet. Die Seriellen Korrelationen $\rho_k$ f\"ur Lag $k$ der + Interspikeintervalle $T_i$ sind die Korrelationskoeffizienten + zwischen den Interspikeintervallen $T_i$ und den um das Lag $k$ + verschobenen Intervallen $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) \] + + 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 - from the spike times returned by the first function. - \begin{solution} - \lstinputlisting{../code/isis.m} - \end{solution} + \part Schreibe eine Funktion, die aus Spikezeiten + Histogramme aus der Anzahl von Spikes, die in Fenstern gegebener L\"ange $W$ + gez\"ahlt werden, erzeugt und plottet. + + 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} \ No newline at end of file diff --git a/pointprocesses/exercises/pointprocesses02.tex b/pointprocesses/exercises/pointprocesses02.tex new file mode 100644 index 0000000..c943457 --- /dev/null +++ b/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} diff --git a/pointprocesses/exercises/serialcorr.pdf b/pointprocesses/exercises/serialcorr.pdf new file mode 100644 index 0000000..5ce399b Binary files /dev/null and b/pointprocesses/exercises/serialcorr.pdf differ diff --git a/pointprocesses/exercises/spikeraster.pdf b/pointprocesses/exercises/spikeraster.pdf new file mode 100644 index 0000000..68c6e28 Binary files /dev/null and b/pointprocesses/exercises/spikeraster.pdf differ diff --git a/pointprocesses/lecture/pointprocesses-chapter.tex b/pointprocesses/lecture/pointprocesses-chapter.tex index 8ecbb10..44af27e 100644 --- a/pointprocesses/lecture/pointprocesses-chapter.tex +++ b/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!} \] \ No newline at end of file diff --git a/pointprocesses/lecture/pointprocesses.tex b/pointprocesses/lecture/pointprocesses.tex index b654a20..4777719 100644 --- a/pointprocesses/lecture/pointprocesses.tex +++ b/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} -Scatter plot von aufeinander folgenden Intervallen $(T_{i+k}, T_i)$ getrennt durch das ``lag'' $k$. +\subsection{Korrelationen der Intervalle} +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} -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)} \] -$\rho_0=1$ (Korrelation jedes Intervalls mit sich selber). +Solche Ab\"angigkeiten werden durch die serielle Korrelation der +Intervalle quantifiziert. Das ist der 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). %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \section{Z\"ahlstatistik} -\begin{figure}[t] - \includegraphics[width=0.48\textwidth]{poissoncounthist100hz10ms}\hfill - \includegraphics[width=0.48\textwidth]{poissoncounthist100hz100ms} - \caption{\label{countstatsfig}Count Statistik.} -\end{figure} +% \begin{figure}[t] +% \includegraphics[width=0.48\textwidth]{poissoncounthist100hz10ms}\hfill +% \includegraphics[width=0.48\textwidth]{poissoncounthist100hz100ms} +% \caption{\label{countstatsfig}Count Statistik.} +% \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} - 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.} + \includegraphics[width=1\textwidth]{poissonraster100hz} + \caption{\label{hompoissonfig}Rasterplot von Spikes eine homogenen + Poisson Prozesse mit $\lambda=100$\,Hz.} \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} diff --git a/programming/code/firing_rates.py b/programming/code/firing_rates.py new file mode 100644 index 0000000..b38dac8 --- /dev/null +++ b/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) diff --git a/programming/exercises/plotting_psth.tex b/programming/exercises/plotting_psth.tex new file mode 100644 index 0000000..8dcee53 --- /dev/null +++ b/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} \ No newline at end of file diff --git a/programming/lectures/images/binned_rate.png b/programming/lectures/images/binned_rate.png new file mode 100644 index 0000000..5500a37 Binary files /dev/null and b/programming/lectures/images/binned_rate.png differ diff --git a/programming/lectures/images/convolved_rate.png b/programming/lectures/images/convolved_rate.png new file mode 100644 index 0000000..0af5699 Binary files /dev/null and b/programming/lectures/images/convolved_rate.png differ diff --git a/programming/lectures/images/instantaneous_rate.png b/programming/lectures/images/instantaneous_rate.png new file mode 100644 index 0000000..3613600 Binary files /dev/null and b/programming/lectures/images/instantaneous_rate.png differ diff --git a/programming/lectures/plotting_spike_trains.tex b/programming/lectures/plotting_spike_trains.tex index 13d9e06..0926fbd 100644 --- a/programming/lectures/plotting_spike_trains.tex +++ b/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 - werden k\"onnte. - \item Erzeugt einen Rasterplot. + \item Erzeugt einen sch\"onen Rasterplot der Zellantworten, speichert ihn. + \item Welche Information liefert er, welche Information ist schwer + 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} - \framesubtitle{Feuerrate \"uber - die Zeit} + \frametitle{Spiketrain Analyse} + \framesubtitle{Zeitabh\"angige Feuerrate, PSTH --- Instantane Feuerrate ---} + \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} - \framesubtitle{Feuerrate \"uber - die Zeit} + \frametitle{Spiketrain Analyse} + \framesubtitle{Zeitabh\"angige Feuerrate, PSTH --- Instantane Feuerrate ---} + Berechnung des PSTHs durch die Instantane Feurerrate: + \textbf{Vorteile:} + \begin{enumerate} - \item Z.B. das \textbf{P}eri- \textbf{S}timulus - \textbf{T}ime - - \textbf{H}istogram, \textbf{PSTH} - \item Wird berechnet indem die Zeit in ``bins'' geteilt wird die - Anzahl Spikes pro bin gezaehlt werden. - \item Die Anzahl wird dann in eine Feuerrate umgerechnet. - \end{enumerate}\pause + \item Sehr einfach zu Berechnen. + \item Macht keine Annahmen \"uber ein Zeitraster, oder die Zeitskala + der neuronalen Verarbeitung. + \end{enumerate} + + \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.5\columnwidth]{images/binning} + \includegraphics[width=0.9\columnwidth]{images/binned_rate} \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} - \framesubtitle{Feuerrate \"uber die Zeit} + \frametitle{Spiketrain Analyse} + \framesubtitle{Zeitabh\"angige Feuerrate, PSTH --- Faltungsmethode ---} + \large Algortihmus: + \normalsize \begin{enumerate} - \item Z.B. das \textbf{P}eri- \textbf{S}timulus - \textbf{T}ime - - \textbf{H}istogram, \textbf{PSTH} - \item Wird berechnet indem man die Aktivit\"at bin\"ar ausdr\"uckt. - \item Jede 1 wird dann duch einen ``Kern'' ersetzt. - \item Der Vorgang heisst Faltung (\verb+conv+). - \end{enumerate}\pause + \item Die neuronalen Antworten werden ``bin\"ar'' ausgedr\"uckt. + \item Ein Filterkern wird berechnet, der das Integral 1 hat. + \item Mithilfe der Faltung (\code{conv} Funktion) wird jede 1 durch + den ``Kern'' ersetzt.\\ + \code{conv(x, kern, 'mode', 'same')} + \end{enumerate} +\end{frame} + + +\begin{frame} + \frametitle{Spiketrain Analyse} + \framesubtitle{Feurerrate als Funktion der Zeit --- Faltungsmethode ---} \begin{figure} - \centering - \includegraphics[width=0.5\columnwidth]{images/conv} + \includegraphics[width=0.9\columnwidth]{images/convolved_rate} \end{figure} \end{frame} -\begin{frame} [fragile] - \frametitle{Spiketrain Analyse I} - \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} +\begin{frame} + \frametitle{Spiketrain Analyse} + \framesubtitle{Feurerrate als Funktion der Zeit --- Faltungsmethode ---} + \textbf{Vorteile:} + \begin{enumerate} + \item Sehr ``nat\"urliche'' erscheinende Darstellung. + \item Sehr gutes Verhalten im Fourrier Raum. + \end{enumerate} + + \textbf{Nachteile:} + \begin{enumerate} + \item Relativ rechenintensiv. + \item Macht Annahmen \"uber die Zeitskalen neuronaler Verarbeitung. + \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} \ No newline at end of file diff --git a/statistics/code/boltzmann.m b/regression/code/boltzmann.m similarity index 100% rename from statistics/code/boltzmann.m rename to regression/code/boltzmann.m diff --git a/statistics/code/create_linear_data.m b/regression/code/create_linear_data.m similarity index 100% rename from statistics/code/create_linear_data.m rename to regression/code/create_linear_data.m diff --git a/statistics/code/estimate_regression.m b/regression/code/estimate_regression.m similarity index 100% rename from statistics/code/estimate_regression.m rename to regression/code/estimate_regression.m diff --git a/statistics/code/exponential.m b/regression/code/exponential.m similarity index 100% rename from statistics/code/exponential.m rename to regression/code/exponential.m diff --git a/scientificcomputing-script.tex b/scientificcomputing-script.tex index a053861..65e8b2b 100644 --- a/scientificcomputing-script.tex +++ b/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} diff --git a/statistics/code/lsq_gradient_sigmoid.m b/statistics/code/lsq_gradient_sigmoid.m deleted file mode 100644 index 05e9721..0000000 --- a/statistics/code/lsq_gradient_sigmoid.m +++ /dev/null @@ -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 \ No newline at end of file diff --git a/statistics/code/lsq_sigmoid_error.m b/statistics/code/lsq_sigmoid_error.m deleted file mode 100644 index 5f05f21..0000000 --- a/statistics/code/lsq_sigmoid_error.m +++ /dev/null @@ -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); \ No newline at end of file diff --git a/statistics/code/sigmoidal_gradient_descent.m b/statistics/code/sigmoidal_gradient_descent.m deleted file mode 100644 index 9763f5a..0000000 --- a/statistics/code/sigmoidal_gradient_descent.m +++ /dev/null @@ -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)