= vthresh :
+ v = vreset
+ times.append(k*dt)
+ spikes.append( times )
+ return spikes
+
+def isis( spikes ) :
+ isi = []
+ for k in xrange(len(spikes)) :
+ isi.extend(np.diff(spikes[k]))
+ return isi
+
+def plotisih( ax, isis, binwidth=None ) :
+ if binwidth == None :
+ nperbin = 200.0 # average number of isis per bin
+ bins = len(isis)/nperbin # number of bins
+ binwidth = np.max(isis)/bins
+ if binwidth < 5e-4 : # half a millisecond
+ binwidth = 5e-4
+ h, b = np.histogram(isis, np.arange(0.0, np.max(isis)+binwidth, binwidth), density=True)
+ ax.text(0.9, 0.85, 'rate={:.0f}Hz'.format(1.0/np.mean(isis)), ha='right', transform=ax.transAxes)
+ ax.text(0.9, 0.75, 'mean={:.0f}ms'.format(1000.0*np.mean(isis)), ha='right', transform=ax.transAxes)
+ ax.text(0.9, 0.65, 'CV={:.2f}'.format(np.std(isis)/np.mean(isis)), ha='right', transform=ax.transAxes)
+ ax.set_xlabel('ISI [ms]')
+ ax.set_ylabel('p(ISI) [1/s]')
+ ax.bar( 1000.0*b[:-1], h, 1000.0*np.diff(b) )
+
+# parameter:
+rate = 20.0
+drate = 50.0
+trials = 10
+duration = 100.0
+dt = 0.001
+tau = 0.1;
+
+# homogeneous spike trains:
+homspikes = hompoisson(rate, trials, duration)
+
+# OU noise:
+rng = np.random.RandomState(54637281)
+time = np.arange(0.0, duration, dt)
+x = np.zeros(time.shape)+rate
+n = rng.randn(len(time))*drate*tau/np.sqrt(dt)+rate
+for k in xrange(1,len(x)) :
+ x[k] = x[k-1] + (n[k]-x[k-1])*dt/tau
+x[x<0.0] = 0.0
+
+# pif spike trains:
+inhspikes = pifspikes(x, trials, dt, D=0.3)
+
+fig = plt.figure( figsize=(9,4) )
+ax = fig.add_subplot(1, 2, 1)
+ax.set_title('stationary')
+ax.set_xlim(0.0, 200.0)
+ax.set_ylim(0.0, 40.0)
+plotisih(ax, isis(homspikes))
+
+ax = fig.add_subplot(1, 2, 2)
+ax.set_title('non-stationary')
+ax.set_xlim(0.0, 200.0)
+ax.set_ylim(0.0, 40.0)
+plotisih(ax, isis(inhspikes))
+
+plt.tight_layout()
+plt.savefig('isihexamples.pdf')
+plt.show()
diff --git a/pointprocesses/lecture/pointprocesses-chapter.tex b/pointprocesses/lecture/pointprocesses-chapter.tex
new file mode 100644
index 0000000..8ecbb10
--- /dev/null
+++ b/pointprocesses/lecture/pointprocesses-chapter.tex
@@ -0,0 +1,271 @@
+\documentclass[12pt]{report}
+
+%%%%% title %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+\title{\tr{Introduction to Scientific Computing}{Einf\"uhrung in die wissenschaftliche Datenverarbeitung}}
+\author{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{pointprocesses}
+
+\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-slides.tex b/pointprocesses/lecture/pointprocesses-slides.tex
new file mode 100644
index 0000000..ff49cb2
--- /dev/null
+++ b/pointprocesses/lecture/pointprocesses-slides.tex
@@ -0,0 +1,412 @@
+\documentclass{beamer}
+
+%%%%% title %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+\title[]{Scientific Computing --- Point Processes}
+\author[]{Jan Benda}
+\institute[]{Neuroethology}
+\date[]{WS 14/15}
+\titlegraphic{\includegraphics[width=0.3\textwidth]{UT_WBMW_Rot_RGB}}
+
+%%%%% beamer %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+\mode = vthresh :
+ v = vreset
+ times.append(k*dt)
+ spikes.append( times )
+ return spikes
+
+# parameter:
+rate = 20.0
+drate = 50.0
+trials = 10
+duration = 2.0
+dt = 0.001
+tau = 0.1;
+
+# homogeneous spike trains:
+homspikes = hompoisson(rate, trials, duration)
+
+# OU noise:
+rng = np.random.RandomState(54637281)
+time = np.arange(0.0, duration, dt)
+x = np.zeros(time.shape)+rate
+n = rng.randn(len(time))*drate*tau/np.sqrt(dt)+rate
+for k in xrange(1,len(x)) :
+ x[k] = x[k-1] + (n[k]-x[k-1])*dt/tau
+x[x<0.0] = 0.0
+
+# inhomogeneous spike trains:
+#inhspikes = inhompoisson(x, trials, dt)
+# pif spike trains:
+inhspikes = pifspikes(x, trials, dt, D=0.3)
+
+fig = plt.figure( figsize=(9,4) )
+ax = fig.add_subplot(1, 2, 1)
+ax.set_title('stationary')
+ax.set_xlim(0.0, duration)
+ax.set_ylim(-0.5, trials-0.5)
+ax.set_xlabel('Time [s]')
+ax.set_ylabel('Trials')
+ax.eventplot(homspikes, colors=[[0, 0, 0]], linelength=0.8)
+
+ax = fig.add_subplot(1, 2, 2)
+ax.set_title('non-stationary')
+ax.set_xlim(0.0, duration)
+ax.set_ylim(-0.5, trials-0.5)
+ax.set_xlabel('Time [s]')
+ax.set_ylabel('Trials')
+ax.eventplot(inhspikes, colors=[[0, 0, 0]], linelength=0.8)
+
+plt.tight_layout()
+plt.savefig('rasterexamples.pdf')
+plt.show()
diff --git a/pointprocesses/lecture/returnmapexamples.py b/pointprocesses/lecture/returnmapexamples.py
new file mode 100644
index 0000000..96803ff
--- /dev/null
+++ b/pointprocesses/lecture/returnmapexamples.py
@@ -0,0 +1,105 @@
+import numpy as np
+import matplotlib.pyplot as plt
+
+def hompoisson(rate, trials, duration) :
+ spikes = []
+ for k in range(trials) :
+ times = []
+ t = 0.0
+ while t < duration :
+ t += np.random.exponential(1/rate)
+ times.append( t )
+ spikes.append( times )
+ return spikes
+
+def inhompoisson(rate, trials, dt) :
+ spikes = []
+ p = rate*dt
+ for k in range(trials) :
+ x = np.random.rand(len(rate))
+ times = dt*np.nonzero(x = vthresh :
+ v = vreset
+ times.append(k*dt)
+ spikes.append( times )
+ return spikes
+
+def isis( spikes ) :
+ isi = []
+ for k in xrange(len(spikes)) :
+ isi.extend(np.diff(spikes[k]))
+ return np.array( isi )
+
+def plotisih( ax, isis, binwidth=None ) :
+ if binwidth == None :
+ nperbin = 200.0 # average number of isis per bin
+ bins = len(isis)/nperbin # number of bins
+ binwidth = np.max(isis)/bins
+ if binwidth < 5e-4 : # half a millisecond
+ binwidth = 5e-4
+ h, b = np.histogram(isis, np.arange(0.0, np.max(isis)+binwidth, binwidth), density=True)
+ ax.text(0.9, 0.85, 'rate={:.0f}Hz'.format(1.0/np.mean(isis)), ha='right', transform=ax.transAxes)
+ ax.text(0.9, 0.75, 'mean={:.0f}ms'.format(1000.0*np.mean(isis)), ha='right', transform=ax.transAxes)
+ ax.text(0.9, 0.65, 'CV={:.2f}'.format(np.std(isis)/np.mean(isis)), ha='right', transform=ax.transAxes)
+ ax.set_xlabel('ISI [ms]')
+ ax.set_ylabel('p(ISI) [1/s]')
+ ax.bar( 1000.0*b[:-1], h, 1000.0*np.diff(b) )
+
+def plotreturnmap(ax, isis, lag=1, max=None) :
+ ax.set_xlabel(r'ISI$_i$ [ms]')
+ ax.set_ylabel(r'ISI$_{i+1}$ [ms]')
+ if max != None :
+ ax.set_xlim(0.0, 1000.0*max)
+ ax.set_ylim(0.0, 1000.0*max)
+ ax.scatter( 1000.0*isis[:-lag], 1000.0*isis[lag:] )
+
+# parameter:
+rate = 20.0
+drate = 50.0
+trials = 10
+duration = 10.0
+dt = 0.001
+tau = 0.1;
+
+# homogeneous spike trains:
+homspikes = hompoisson(rate, trials, duration)
+
+# OU noise:
+rng = np.random.RandomState(54637281)
+time = np.arange(0.0, duration, dt)
+x = np.zeros(time.shape)+rate
+n = rng.randn(len(time))*drate*tau/np.sqrt(dt)+rate
+for k in xrange(1,len(x)) :
+ x[k] = x[k-1] + (n[k]-x[k-1])*dt/tau
+x[x<0.0] = 0.0
+
+# pif spike trains:
+inhspikes = pifspikes(x, trials, dt, D=0.3)
+
+fig = plt.figure( figsize=(9,4) )
+ax = fig.add_subplot(1, 2, 1)
+ax.set_title('stationary')
+plotreturnmap(ax, isis(homspikes), 1, 0.3)
+
+ax = fig.add_subplot(1, 2, 2)
+ax.set_title('non-stationary')
+plotreturnmap(ax, isis(inhspikes), 1, 0.3)
+
+plt.tight_layout()
+plt.savefig('returnmapexamples.pdf')
+#plt.show()
diff --git a/pointprocesses/lecture/serialcorrexamples.py b/pointprocesses/lecture/serialcorrexamples.py
new file mode 100644
index 0000000..e7fff9c
--- /dev/null
+++ b/pointprocesses/lecture/serialcorrexamples.py
@@ -0,0 +1,117 @@
+import numpy as np
+import matplotlib.pyplot as plt
+
+def hompoisson(rate, trials, duration) :
+ spikes = []
+ for k in range(trials) :
+ times = []
+ t = 0.0
+ while t < duration :
+ t += np.random.exponential(1/rate)
+ times.append( t )
+ spikes.append( times )
+ return spikes
+
+def inhompoisson(rate, trials, dt) :
+ spikes = []
+ p = rate*dt
+ for k in range(trials) :
+ x = np.random.rand(len(rate))
+ times = dt*np.nonzero(x = vthresh :
+ v = vreset
+ times.append(k*dt)
+ spikes.append( times )
+ return spikes
+
+def isis( spikes ) :
+ isi = []
+ for k in xrange(len(spikes)) :
+ isi.extend(np.diff(spikes[k]))
+ return np.array( isi )
+
+def plotisih( ax, isis, binwidth=None ) :
+ if binwidth == None :
+ nperbin = 200.0 # average number of isis per bin
+ bins = len(isis)/nperbin # number of bins
+ binwidth = np.max(isis)/bins
+ if binwidth < 5e-4 : # half a millisecond
+ binwidth = 5e-4
+ h, b = np.histogram(isis, np.arange(0.0, np.max(isis)+binwidth, binwidth), density=True)
+ ax.text(0.9, 0.85, 'rate={:.0f}Hz'.format(1.0/np.mean(isis)), ha='right', transform=ax.transAxes)
+ ax.text(0.9, 0.75, 'mean={:.0f}ms'.format(1000.0*np.mean(isis)), ha='right', transform=ax.transAxes)
+ ax.text(0.9, 0.65, 'CV={:.2f}'.format(np.std(isis)/np.mean(isis)), ha='right', transform=ax.transAxes)
+ ax.set_xlabel('ISI [ms]')
+ ax.set_ylabel('p(ISI) [1/s]')
+ ax.bar( 1000.0*b[:-1], h, 1000.0*np.diff(b) )
+
+def plotreturnmap(ax, isis, lag=1, max=None) :
+ ax.set_xlabel(r'ISI$_i$ [ms]')
+ ax.set_ylabel(r'ISI$_{i+1}$ [ms]')
+ if max != None :
+ ax.set_xlim(0.0, 1000.0*max)
+ ax.set_ylim(0.0, 1000.0*max)
+ ax.scatter( 1000.0*isis[:-lag], 1000.0*isis[lag:] )
+
+def plotserialcorr(ax, isis, maxlag=10) :
+ lags = np.arange(maxlag+1)
+ corr = [1.0]
+ for lag in lags[1:] :
+ corr.append(np.corrcoef(isis[:-lag], isis[lag:])[0,1])
+ ax.set_xlabel(r'lag $k$')
+ ax.set_ylabel(r'ISI correlation $\rho_k$')
+ ax.set_xlim(0.0, maxlag)
+ ax.set_ylim(-1.0, 1.0)
+ ax.plot(lags, corr, '.-', markersize=20)
+
+# parameter:
+rate = 20.0
+drate = 50.0
+trials = 10
+duration = 500.0
+dt = 0.001
+tau = 0.1;
+
+# homogeneous spike trains:
+homspikes = hompoisson(rate, trials, duration)
+
+# OU noise:
+rng = np.random.RandomState(54637281)
+time = np.arange(0.0, duration, dt)
+x = np.zeros(time.shape)+rate
+n = rng.randn(len(time))*drate*tau/np.sqrt(dt)+rate
+for k in xrange(1,len(x)) :
+ x[k] = x[k-1] + (n[k]-x[k-1])*dt/tau
+x[x<0.0] = 0.0
+
+# pif spike trains:
+inhspikes = pifspikes(x, trials, dt, D=0.3)
+
+fig = plt.figure( figsize=(9,3) )
+
+ax = fig.add_subplot(1, 2, 1)
+plotserialcorr(ax, isis(homspikes))
+ax.set_ylim(-0.2, 1.0)
+
+ax = fig.add_subplot(1, 2, 2)
+plotserialcorr(ax, isis(inhspikes))
+ax.set_ylim(-0.2, 1.0)
+
+plt.tight_layout()
+plt.savefig('serialcorrexamples.pdf')
+#plt.show()
diff --git a/programming/lectures/images/badbarleft.png b/programming/lectures/images/badbarleft.png
new file mode 100644
index 0000000..667e029
Binary files /dev/null and b/programming/lectures/images/badbarleft.png differ
diff --git a/programming/lectures/images/badbarplot.jpg b/programming/lectures/images/badbarplot.jpg
new file mode 100644
index 0000000..f9327de
Binary files /dev/null and b/programming/lectures/images/badbarplot.jpg differ
diff --git a/programming/lectures/images/badbarright.png b/programming/lectures/images/badbarright.png
new file mode 100644
index 0000000..22f32dc
Binary files /dev/null and b/programming/lectures/images/badbarright.png differ
diff --git a/programming/lectures/images/comparison_properly_improperly_graph.png b/programming/lectures/images/comparison_properly_improperly_graph.png
new file mode 100644
index 0000000..eaae9c2
Binary files /dev/null and b/programming/lectures/images/comparison_properly_improperly_graph.png differ
diff --git a/programming/lectures/images/histogram.png b/programming/lectures/images/histogram.png
new file mode 100755
index 0000000..8eaaa88
Binary files /dev/null and b/programming/lectures/images/histogram.png differ
diff --git a/programming/lectures/images/histogrambad.png b/programming/lectures/images/histogrambad.png
new file mode 100755
index 0000000..3ce2c65
Binary files /dev/null and b/programming/lectures/images/histogrambad.png differ
diff --git a/programming/lectures/images/histogrambad2.png b/programming/lectures/images/histogrambad2.png
new file mode 100755
index 0000000..a11dc47
Binary files /dev/null and b/programming/lectures/images/histogrambad2.png differ
diff --git a/programming/lectures/images/improperly_scaled_graph.png b/programming/lectures/images/improperly_scaled_graph.png
new file mode 100644
index 0000000..ac99274
Binary files /dev/null and b/programming/lectures/images/improperly_scaled_graph.png differ
diff --git a/programming/lectures/images/line_graph1.png b/programming/lectures/images/line_graph1.png
new file mode 100644
index 0000000..afce7a1
Binary files /dev/null and b/programming/lectures/images/line_graph1.png differ
diff --git a/programming/lectures/images/line_graph1_3.png b/programming/lectures/images/line_graph1_3.png
new file mode 100644
index 0000000..3a4a528
Binary files /dev/null and b/programming/lectures/images/line_graph1_3.png differ
diff --git a/programming/lectures/images/line_graph1_4.png b/programming/lectures/images/line_graph1_4.png
new file mode 100644
index 0000000..d13e30f
Binary files /dev/null and b/programming/lectures/images/line_graph1_4.png differ
diff --git a/programming/lectures/images/misleading_pie.png b/programming/lectures/images/misleading_pie.png
new file mode 100644
index 0000000..3dae0cd
Binary files /dev/null and b/programming/lectures/images/misleading_pie.png differ
diff --git a/programming/lectures/images/nobelbad.png b/programming/lectures/images/nobelbad.png
new file mode 100644
index 0000000..4602b7c
Binary files /dev/null and b/programming/lectures/images/nobelbad.png differ
diff --git a/programming/lectures/images/one_d_problem_c.pdf b/programming/lectures/images/one_d_problem_c.pdf
new file mode 100644
index 0000000..31974db
Binary files /dev/null and b/programming/lectures/images/one_d_problem_c.pdf differ
diff --git a/programming/lectures/images/properly_scaled_graph.png b/programming/lectures/images/properly_scaled_graph.png
new file mode 100644
index 0000000..1a5de32
Binary files /dev/null and b/programming/lectures/images/properly_scaled_graph.png differ
diff --git a/programming/lectures/images/sample_pie.png b/programming/lectures/images/sample_pie.png
new file mode 100644
index 0000000..937a231
Binary files /dev/null and b/programming/lectures/images/sample_pie.png differ
diff --git a/programming/lectures/plotting_spike_trains.tex b/programming/lectures/plotting_spike_trains.tex
index cb85c51..13d9e06 100644
--- a/programming/lectures/plotting_spike_trains.tex
+++ b/programming/lectures/plotting_spike_trains.tex
@@ -86,6 +86,8 @@
\end{flushright}
}
+\newcommand{\code}[1]{\texttt{#1}}
+
\input{../../latex/environments.tex}
\makeatother
@@ -103,11 +105,259 @@
\begin{enumerate}
\item Graphische Darstellung von Daten
\item Spiketrain Analyse
- \item \"Ubungen, \"Ubungen, \"Ubungen.
\end{enumerate}
\end{frame}
+\begin{frame}[plain]
+ \huge{1. Graphische Darstellung von Daten}\pause
+ \begin{figure}
+ \includegraphics[width=0.9\columnwidth]{images/convincing}
+ \end{figure}
+\end{frame}
+
+\begin{frame}
+ \frametitle{Graphische Darstellung von Daten}
+ \framesubtitle{Was soll ein Datenplot erreichen?}
+
+ \begin{itemize}
+ \item Ist eine m\"oglichst neutrale Darstellung der Daten.
+ \item Soll dem Leser die Daten greifbar machen und die Aussagen der
+ Analyse darstellen.
+ \item Erlaubt dem Leser die gezeigten Effekte selbst zu beguachten
+ und zu validieren.
+ \item Muss vollst\"andig annotiert sein.
+ \item Folgt dem Prinzip der \textbf{ink minimization}. (Das
+ Verh\"altnis aus Tinte, die f\"ur die Darstellung der Daten
+ gebraucht wird und der Menge Tinte, die f\"ur die Graphik
+ ben\"otigt wird sollte m\"oglichst gro{\ss} sein )
+ \end{itemize}
+\end{frame}
+
+\begin{frame}
+ \frametitle{Graphische Darstellung von Daten}
+ \framesubtitle{Was sollte vermieden werden?}
+
+ \begin{itemize}
+ \item Suggestive oder gar fehlleitende Darstellung.
+ \item Ablenkung durch unruhige oder \"uberm\"a{\ss}ige Effekte.
+ \item Comicartige Effekte...
+ \end{itemize}\pause
+ \begin{figure}
+ \includegraphics[width=0.35\columnwidth]{images/one_d_problem_c}
+ \end{figure}\pause
+ ... aus{\ss}er sie werden rein zur Illustration benutzt ohne einen
+ Anspruch auf Richtigkeit zu erheben.
+\end{frame}
+
+\begin{frame}
+ \frametitle{Graphische Darstellung von Daten}
+ \framesubtitle{Suboptimale Beispiele}
+ \only <1> {
+ \begin{figure}
+ \includegraphics[width=0.5\columnwidth]{images/nobelbad}
+ \end{figure}
+ \vspace{0.25cm}
+ Aus Hafting et al., Nature, 2005
+ }
+ \only <2> {
+ \begin{figure}
+ \includegraphics[width=0.7\columnwidth]{images/misleading_pie}
+ \end{figure}
+ \vspace{0.5cm}
+ \url{https://en.wikipedia.org/wiki/Misleading_graph}
+ }
+ \only <3> {
+ \begin{figure}
+ \includegraphics[width=0.7\columnwidth]{images/sample_pie}
+ \end{figure}
+ \vspace{0.5cm}
+ \url{https://en.wikipedia.org/wiki/Misleading_graph}
+ }
+ \only <4> {
+ \begin{figure}
+ \includegraphics[width=0.4\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}
+ \end{figure}
+ \vspace{0.5cm}
+ \url{https://en.wikipedia.org/wiki/Misleading_graph}
+ }
+ \only <6> {
+ \begin{figure}
+ \includegraphics[width=0.8\columnwidth]{images/badbarplot}
+ \end{figure}
+ \vspace{0.5cm}
+ \url{https://en.wikipedia.org/wiki/Misleading_graph}
+ }
+ \only <7> {
+ Wahl der Zeichenfl\"ache kann den visuellen Eindruck beeinflu{\ss}en.
+ \begin{columns}
+ \begin{column}{4.cm}
+ \begin{figure}
+ \includegraphics[width=0.7\columnwidth]{images/line_graph1}
+ \end{figure}
+ \end{column}
+
+ \begin{column}{4.cm}
+ \begin{figure}
+ \includegraphics[width=0.7\columnwidth]{images/line_graph1_3}
+ \end{figure}
+ \end{column}
+
+ \begin{column}{4.cm}
+ \begin{figure}
+ \includegraphics[width=0.7\columnwidth]{images/line_graph1_4}
+ \end{figure}
+ \end{column}
+ \end{columns}
+ \vspace{0.5cm}
+ \url{https://en.wikipedia.org/wiki/Misleading_graph}
+ }
+ \only <8> {
+ Vorsicht bei der Skalierung von Symbolen!
+ \begin{columns}
+ \begin{column}{4.cm}
+ \begin{figure}
+ \includegraphics[width=0.7\columnwidth]{images/improperly_scaled_graph}
+ \end{figure}
+ \end{column}
+
+ \begin{column}{4.cm}
+ \begin{figure}
+ \includegraphics[width=0.7\columnwidth]{images/properly_scaled_graph}
+ \end{figure}
+ \end{column}
+
+ \begin{column}{4.cm}
+ \begin{figure}
+ \includegraphics[width=0.7\columnwidth]{images/comparison_properly_improperly_graph}
+ \end{figure}
+ \end{column}
+ \end{columns}
+ \vspace{0.5cm}
+ \url{https://en.wikipedia.org/wiki/Misleading_graph}
+ }
+\end{frame}
+
+\begin{frame}[fragile]
+ \frametitle{Graphische Darstellung von Daten}
+ \framesubtitle{Plotting Interfaces in Matlab}
+ Es gibt zwei Wege Graphen zu bearbeiten:
+ \begin{enumerate}
+ \item Interaktiv \"uber das \textit{graphische User Interface}\pause
+ \item Die Kommandozeile bzw. in Skripten und Funktionen.\pause
+ \end{enumerate}
+ Beides hat seine Berechtigung und seine eigenen Vor- und Nachteile. Welche?
+\end{frame}
+
+
+\begin{frame}
+ \frametitle{Graphische Darstellung von Daten}
+ \framesubtitle{Ver\"anderung des Graphen \"uber die Kommandozeile}
+ \begin{itemize}
+ \item Erstellt ein Skript, dass einen Plot erstellt.
+ \item Dieser soll zwei Sinus unterschiedlicher Frequenz darstellen.
+ \end{itemize}
+ Wir werden jetzt die Kommandozeil bzw. das Skript verbessern um den
+ Plot ``sch\"oner'' zu machen.
+\end{frame}
+
+
+\begin{frame}
+ \frametitle{Graphische Darstellung von Daten}
+ \framesubtitle{Ver\"anderung des Graphen \"uber die Kommandozeile}
+ \begin{enumerate}
+ \item Einstellungen der Linienplots:
+ \begin{itemize}
+ \item St\"arke und Farbe.
+ \item Linienstil, Marker.
+ \end{itemize}\pause
+ \item Achsbeschriftung:
+ \begin{itemize}
+ \item \code{xlabel}, \code{ylabel}.
+ \item Schriftart und Gr\"o{\ss}e.
+ \end{itemize}\pause
+ \item Achsenskalierung und Ticks:
+ \begin{itemize}
+ \item Skalierung der Achsen (Minumum und Maxmimum, logarithmisch oder linear).
+ \item Manuelles Setzen der Ticks, ihrer Richtung und Beschriftung.
+ \item Grid or no Grid?
+ \end{itemize}\pause
+ \item Setzen von globalen Parametern:
+ \begin{itemize}
+ \item Einstellung der Papiergr\"o{\ss}e und plzieren der
+ Zeichenfl\"ache.
+ \item Box oder nicht?
+ \item Speichern der Abbildung als pdf.
+ \end{itemize}
+ \end{enumerate}
+\end{frame}
+
+
+\begin{frame} [fragile]
+ \frametitle{Graphische Darstellung von Daten}
+ \framesubtitle{Ver\"andern von Eigenschaften \"uber die Kommandozeile}
+ \vspace{-0.75em}
+ \scriptsize
+ \begin{lstlisting}
+fig = figure();
+set(gcf, 'PaperUnits', 'centimeters', 'PaperSize', [11.7 9.0]);
+set(gcf, 'PaperPosition',[0.0 0.0 11.7 9.0], 'Color', 'white')
+hold on
+plot(time, neuronal_data, 'color', [ 0.2 0.5 0.7], 'linewidth', 1.)
+plot(spike_times, ones(size(spike_times))*threshold, 'ro', 'markersize', 4)
+line([time(1) time(end)], [threshold threshold], 'linestyle', '--',
+ 'linewidth', 0.75, 'color', [0.9 0.9 0.9])
+ylim([0 35])
+xlim([0 2.25])
+box('off')
+xlabel('time [s]', 'fontname', 'MyriadPro-Regular', 'fontsize', 10)
+ylabel('potential [mV]', 'fontname', 'MyriadPro-Regular', 'fontsize', 10)
+title('pyramidal cell', 'fontname', 'MyriadPro-Regular', 'fontsize', 12)
+set(gca, 'TickDir','out', 'linewidth', 1.5, 'fontname', 'MyriadPro-Regular')
+saveas(fig, 'spike_detection.pdf', 'pdf')
+\end{lstlisting}
+\end{frame}
+
+
+\begin{frame} [fragile]
+ \frametitle{Graphische Darstellung von Daten}
+ \framesubtitle{Ver\"andern von Eigenschaften \"uber die Kommandozeile}
+ \begin{figure}
+ \centering
+ \includegraphics[width=0.75\columnwidth]{./images/spike_detection}
+ \end{figure}
+\end{frame}
+
+
+\begin{frame} [fragile]
+ \frametitle{Graphische Darstellung von Daten}
+ \framesubtitle{Welche Art Plot wof\"ur?}
+ \url{http://www.mathworks.de/discovery/gallery.html}
+\end{frame}
+
+
+\begin{frame} [fragile]
+ \frametitle{Graphische Darstellung von Daten}
+ \framesubtitle{Was macht einen guten Abbildung aus?}
+ \begin{enumerate}
+ \item Klarheit.
+ \item Vollstaendige Beschriftung.
+ \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}
\end{frame}
diff --git a/scientificcomputing-script.tex b/scientificcomputing-script.tex
index 306a84a..a053861 100644
--- a/scientificcomputing-script.tex
+++ b/scientificcomputing-script.tex
@@ -59,7 +59,8 @@
% figures:
\setlength{\fboxsep}{0pt}
-\newcommand{\texpicture}[1]{{\sffamily\footnotesize\input{#1.tex}}}
+\newcommand{\texinputpath}{}
+\newcommand{\texpicture}[1]{{\sffamily\footnotesize\input{\texinputpath#1.tex}}}
%\newcommand{\texpicture}[1]{\fbox{\sffamily\footnotesize\input{#1.tex}}}
%\newcommand{\texpicture}[1]{\setlength{\fboxsep}{2mm}\fbox{#1}}
%\newcommand{\texpicture}[1]{}
@@ -204,8 +205,9 @@
\newenvironment{definition}[1][]{\medskip\noindent\textbf{Definition}\ifthenelse{\equal{#1}{}}{}{ #1}:\newline}%
{\medskip}
+%%%%% exercises: %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\newcounter{maxexercise}
-\setcounter{maxexercise}{9} % show listings up to exercise maxexercise
+\setcounter{maxexercise}{10} % show listings up to exercise maxexercise
\newcounter{theexercise}
\setcounter{theexercise}{1}
\newcommand{\codepath}{}
@@ -213,7 +215,7 @@
\arabic{theexercise}:}\newline \newcommand{\exercisesource}{#1}}%
{\ifthenelse{\equal{\exercisesource}{}}{}{\ifthenelse{\value{theexercise}>\value{maxexercise}}{}{\medskip\lstinputlisting{\codepath\exercisesource}}}\medskip\stepcounter{theexercise}}
-\graphicspath{{statistics/lecture/}{statistics/lecture/figures/}{bootstrap/lecture/}{bootstrap/lecture/figures/}{likelihood/lecture/}{likelihood/lecture/figures/}}
+\graphicspath{{statistics/lecture/}{statistics/lecture/figures/}{bootstrap/lecture/}{bootstrap/lecture/figures/}{likelihood/lecture/}{likelihood/lecture/figures/}{pointprocesses/lecture/}{pointprocesses/lecture/figures/}}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
@@ -233,4 +235,8 @@
\renewcommand{\codepath}{likelihood/code/}
\include{likelihood/lecture/likelihood}
+\renewcommand{\codepath}{pointprocesses/code/}
+\renewcommand{\texinputpath}{pointprocesses/lecture/}
+%\include{pointprocesses/lecture/pointprocesses}
+
\end{document}
diff --git a/statistics/code/boltzmann.m b/statistics/code/boltzmann.m
new file mode 100644
index 0000000..8190460
--- /dev/null
+++ b/statistics/code/boltzmann.m
@@ -0,0 +1,7 @@
+function y = boltzmann(parameter, x)
+% parameter 1: alpha
+% parameter 2: k
+% parameter 3: x_0
+% parameter 4: y_0
+
+y = (parameter(1) ./ (1 + exp(-parameter(2) .* (x - parameter(3))))) + parameter(4);
\ No newline at end of file
diff --git a/statistics/code/create_linear_data.m b/statistics/code/create_linear_data.m
new file mode 100644
index 0000000..f929f23
--- /dev/null
+++ b/statistics/code/create_linear_data.m
@@ -0,0 +1,7 @@
+function y = create_linear_data(x)
+
+m = 2.5;
+n = -0.35;
+d = 2.5;
+
+y = x .* m + n + randn(size(x)) .* d;
diff --git a/statistics/code/estimate_regression.m b/statistics/code/estimate_regression.m
new file mode 100644
index 0000000..fd3d2c4
--- /dev/null
+++ b/statistics/code/estimate_regression.m
@@ -0,0 +1,7 @@
+function param = estimate_regression(x,y,p_0)
+
+objective_function = @(p)lsq_error(p, x, y);
+param = fminunc(objective_function, p_0);
+disp(param)
+param1 = fminsearch(objective_function, p_0);
+disp(param1)
\ No newline at end of file
diff --git a/statistics/code/exponential.m b/statistics/code/exponential.m
new file mode 100644
index 0000000..1bb2a04
--- /dev/null
+++ b/statistics/code/exponential.m
@@ -0,0 +1,5 @@
+function y = exponential(parameter, x)
+% Function implements an exponential function with two parameters
+% controlling the amplitude and the time constant.
+
+y = parameter(1) .* exp(x./parameter(2));
\ No newline at end of file
diff --git a/statistics/code/lsq_gradient_sigmoid.m b/statistics/code/lsq_gradient_sigmoid.m
new file mode 100644
index 0000000..05e9721
--- /dev/null
+++ b/statistics/code/lsq_gradient_sigmoid.m
@@ -0,0 +1,9 @@
+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
new file mode 100644
index 0000000..5f05f21
--- /dev/null
+++ b/statistics/code/lsq_sigmoid_error.m
@@ -0,0 +1,8 @@
+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
new file mode 100644
index 0000000..9763f5a
--- /dev/null
+++ b/statistics/code/sigmoidal_gradient_descent.m
@@ -0,0 +1,44 @@
+
+
+%% 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)
diff --git a/statistics/exercises/mlestd.m b/statistics/exercises/mlestd.m
index 61f4b22..de37a56 100644
--- a/statistics/exercises/mlestd.m
+++ b/statistics/exercises/mlestd.m
@@ -27,4 +27,4 @@ subplot(1, 2, 2);
plot(psigs, loglm);
xlabel('standard deviation')
ylabel('log likelihood')
-savefigpdf(gcf, 'mlestd.pdf', 12, 5);
+savefigpdf(gcf, 'mlestd.pdf', 15, 5);
diff --git a/statistics/exercises/mlestd.pdf b/statistics/exercises/mlestd.pdf
index dad420c..f01f927 100644
Binary files a/statistics/exercises/mlestd.pdf and b/statistics/exercises/mlestd.pdf differ
diff --git a/statistics/exercises/statistics04.tex b/statistics/exercises/statistics04.tex
index 276f7ea..1c3dcdf 100644
--- a/statistics/exercises/statistics04.tex
+++ b/statistics/exercises/statistics04.tex
@@ -113,8 +113,10 @@ Absch\"atzung der Standardabweichung verdeutlichen.
\continue
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\question \qt{Maximum-Likelihood-Sch\"atzer einer Ursprungsgeraden}
-In der Vorlesung haben wir eine Gleichung f\"ur die Maximum-Likelihood
-Absch\"atzung der Steigung einer Ursprungsgeraden hergeleitet.
+In der Vorlesung haben wir folgende Formel f\"ur die Maximum-Likelihood
+Absch\"atzung der Steigung $\theta$ einer Ursprungsgeraden durch $n$ Datenpunkte $(x_i|y_i)$ mit Standardabweichung $\sigma_i$ hergeleitet:
+\[\theta = \frac{\sum_{i=1}^n \frac{x_iy_i}{\sigma_i^2}}{ \sum_{i=1}^n
+ \frac{x_i^2}{\sigma_i^2}} \]
\begin{parts}
\part \label{mleslopefunc} Schreibe eine Funktion, die in einem $x$ und einem
$y$ Vektor die Datenpaare \"uberreicht bekommt und die Steigung der
@@ -146,13 +148,12 @@ nicht so einfach wie der Mittelwert und die Standardabweichung einer
Normalverteilung direkt aus den Daten berechnet werden k\"onnen. Solche Parameter
m\"ussen dann aus den Daten mit der Maximum-Likelihood-Methode gefittet werden.
-Um dies zu veranschaulichen ziehen wir uns diesmal Zufallszahlen, die nicht einer
-Normalverteilung entstammen, sonder aus der Gamma-Verteilung.
+Um dies zu veranschaulichen ziehen wir uns diesmal nicht normalverteilte Zufallszahlen, sondern Zufallszahlen aus der Gamma-Verteilung.
\begin{parts}
\part
- Finde heraus welche Funktion die Wahrscheinlichkeitsdichtefunktion
- (probability density function) der Gamma-Verteilung in \code{matlab}
- berechnet.
+ Finde heraus welche \code{matlab} Funktion die
+ Wahrscheinlichkeitsdichtefunktion (probability density function) der
+ Gamma-Verteilung berechnet.
\part
Plotte mit Hilfe dieser Funktion die Wahrscheinlichkeitsdichtefunktion
@@ -169,17 +170,17 @@ Normalverteilung entstammen, sonder aus der Gamma-Verteilung.
\part
Finde heraus mit welcher \code{matlab}-Funktion eine beliebige
- Verteilung (``distribution'') und die Gammaverteilung an die
- Zufallszahlen nach der Maximum-Likelihood Methode gefittet werden
- kann.
+ Verteilung (``distribution'') an die Zufallszahlen nach der
+ Maximum-Likelihood Methode gefittet werden kann. Wie wird diese
+ Funktion benutzt, um die Gammaverteilung an die Daten zu fitten?
\part
- Bestimme mit dieser Funktion die Parameter der
- Gammaverteilung aus den Zufallszahlen.
+ Bestimme mit dieser Funktion die Parameter der Gammaverteilung aus
+ den Zufallszahlen.
\part
- Plotte anschlie{\ss}end
- die Gammaverteilung mit den gefitteten Parametern.
+ Plotte anschlie{\ss}end die Gammaverteilung mit den gefitteten
+ Parametern.
\end{parts}
\begin{solution}
\lstinputlisting{mlepdffit.m}
diff --git a/statistics/lecture/Makefile b/statistics/lecture/Makefile
index 445913e..ec84405 100644
--- a/statistics/lecture/Makefile
+++ b/statistics/lecture/Makefile
@@ -1,22 +1,29 @@
BASENAME=statistics
+
PYFILES=$(wildcard *.py)
PYPDFFILES=$(PYFILES:.py=.pdf)
-pdf : $(BASENAME)-chapter.pdf $(PYPDFFILES)
+all : pdf
+
+# script:
+pdf : $(BASENAME)-chapter.pdf
-$(BASENAME)-chapter.pdf : $(BASENAME)-chapter.tex $(BASENAME).tex
+$(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 *~ $(BASENAME)-chapter.aux $(BASENAME)-chapter.log $(BASENAME)-chapter.out $(BASENAME).aux $(BASENAME).log
+ 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
-watch :
+watchpdf :
while true; do ! make -q pdf && make pdf; sleep 0.5; done
diff --git a/statistics/lecture/boxwhisker.py b/statistics/lecture/boxwhisker.py
index 1209f7e..38a0dcb 100644
--- a/statistics/lecture/boxwhisker.py
+++ b/statistics/lecture/boxwhisker.py
@@ -43,5 +43,5 @@ ax.annotate('maximum',
ax.boxplot( x, whis=100.0 )
plt.tight_layout()
plt.savefig('boxwhisker.pdf')
-plt.show()
+#plt.show()
diff --git a/statistics/lecture/correlation.py b/statistics/lecture/correlation.py
index db317c7..2fa85d4 100644
--- a/statistics/lecture/correlation.py
+++ b/statistics/lecture/correlation.py
@@ -5,7 +5,6 @@ plt.xkcd()
fig = plt.figure( figsize=(6,5) )
n = 200
for k, r in enumerate( [ 1.0, 0.6, 0.0, -0.9 ] ) :
- print r
x = np.random.randn( n )
y = r*x + np.sqrt(1.0-r*r)*np.random.randn( n )
ax = fig.add_subplot( 2, 2, k+1 )
@@ -30,5 +29,4 @@ for k, r in enumerate( [ 1.0, 0.6, 0.0, -0.9 ] ) :
plt.tight_layout()
plt.savefig('correlation.pdf')
-plt.show()
-
+#plt.show()
diff --git a/statistics/lecture/diehistograms.py b/statistics/lecture/diehistograms.py
index d5e0380..a8832e0 100644
--- a/statistics/lecture/diehistograms.py
+++ b/statistics/lecture/diehistograms.py
@@ -28,5 +28,4 @@ ax.set_ylabel( 'Probability' )
ax.hist([x2, x1], bins, normed=True, color=['#FFCC00', '#FFFF66' ])
plt.tight_layout()
fig.savefig( 'diehistograms.pdf' )
-plt.show()
-
+#plt.show()
diff --git a/statistics/lecture/median.py b/statistics/lecture/median.py
index 2bf420c..a1231f9 100644
--- a/statistics/lecture/median.py
+++ b/statistics/lecture/median.py
@@ -29,5 +29,4 @@ ax.plot(x,g, 'b', lw=4)
ax.plot([0.0, 0.0], [0.0, 0.45], 'k', lw=2 )
plt.tight_layout()
fig.savefig( 'median.pdf' )
-plt.show()
-
+#plt.show()
diff --git a/statistics/lecture/nonlincorrelation.py b/statistics/lecture/nonlincorrelation.py
index c0ca723..498d985 100644
--- a/statistics/lecture/nonlincorrelation.py
+++ b/statistics/lecture/nonlincorrelation.py
@@ -39,4 +39,4 @@ ax.scatter( x, z )
plt.tight_layout()
plt.savefig('nonlincorrelation.pdf')
-plt.show()
+#plt.show()
diff --git a/statistics/lecture/pdfhistogram.py b/statistics/lecture/pdfhistogram.py
index 039b524..a61460e 100644
--- a/statistics/lecture/pdfhistogram.py
+++ b/statistics/lecture/pdfhistogram.py
@@ -35,5 +35,5 @@ ax.hist(r, 20, normed=True, color='#FFCC00')
plt.tight_layout()
fig.savefig( 'pdfhistogram.pdf' )
-plt.show()
+#plt.show()
diff --git a/statistics/lecture/pdfprobabilities.py b/statistics/lecture/pdfprobabilities.py
index 6481da1..bcbaa09 100644
--- a/statistics/lecture/pdfprobabilities.py
+++ b/statistics/lecture/pdfprobabilities.py
@@ -32,5 +32,4 @@ ax.fill_between( x[(x>x1)&(x