Merge branch 'master' of raven:scientificComputing

programming/code/firing_rates.py

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

programming/exercises/plotting_psth.tex

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

programming/lectures/plotting_spike_trains.tex

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

statistics/code/lsq_gradient_sigmoid.m

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

statistics/code/lsq_sigmoid_error.m

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

statistics/code/sigmoidal_gradient_descent.m

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