Fixed plotting Makefile nad pathes.

Copied psth_sta text and scripts to pointprocess chapter.

Makefile

@@ -1,7 +1,7 @@
 BASENAME=scientificcomputing-script
 
 SUBDIRS=designpattern statistics bootstrap regression likelihood pointprocesses 
-#programming  spike_trains
+#programming
 SUBTEXS=$(foreach subd, $(SUBDIRS), $(subd)/lecture/$(subd).tex)
 
 pdf : chapters $(BASENAME).pdf

header.tex

@@ -170,14 +170,15 @@
 \newcommand{\koZ}{\mathds{C}}
 
 %%%%% code/matlab commands: %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+\usepackage{textcomp}
 \newcommand{\code}[1]{\setlength{\fboxsep}{0.5ex}\colorbox{blue!10}{\texttt{#1}}}
-\newcommand{\matlab}{MATLAB}
+\newcommand{\matlab}{MATLAB$^{\copyright}$}
 \newcommand{\matlabfun}[1]{(\tr{\matlab{}-function}{\matlab-Funktion} \setlength{\fboxsep}{0.5ex}\colorbox{blue!10}{\texttt{#1}})}
 
 %%%%% exercises environment: %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 % usage:
 %
-% \begin{exercise}[somecode.m]
+% \begin{exercise}{somecode.m}{someoutput.out}
 %   Write a function that computes the median.
 % \end{exercise}
 %
@@ -187,6 +188,8 @@
 % Write a function that computes the median.
 % Listing 1: somecode.m
 % the listing
+% Output:
+% content of someoutput.out
 %
 % Innerhalb der exercise Umgebung ist enumerate umdefiniert, um (a), (b), (c), .. zu erzeugen.
 \usepackage{ifthen}
@@ -194,8 +197,7 @@
 \newcounter{maxexercise} 
 \setcounter{maxexercise}{10000}  % show listings up to exercise maxexercise
 \newcounter{exercise}[chapter]
-\setcounter{exercise}{0}
+\renewcommand{\theexercise}{\arabic{exercise}}
-\newcommand{\theexercise}{\arabic{exercise}}
 \newcommand{\codepath}{}
 \newenvironment{exercise}[2]%
 {\newcommand{\exercisesource}{#1}%

plotting/lecture/plotting-chapter.tex

@@ -11,7 +11,7 @@
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 \begin{document} 
 
-\include{psth_sta}
+\include{plotting}
 
 \end{document}
 

pointprocesses/lecture/firingrates.py

@@ -0,0 +1,268 @@
+import numpy as np
+import matplotlib.pyplot as plt
+import scipy.io as spio
+import scipy.stats as spst
+import scipy as sp
+from IPython import embed
+
+
+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)
+    ax1.set_xlim([0, 5])
+    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("figures/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("figures/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 = spst.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("figures/convolved_rate.png")
+    plt.close()
+
+
+def plot_comparison(spike_times, bin_width, sigma, max_t=30., dt=1e-4):
+    times = np.squeeze(spike_times[0][0])
+    time, conv_rate = get_convolved_rate(times, bin_width/np.sqrt(12.))
+    time, inst_rate = get_instantaneous_rate(times)
+    time, binn_rate = get_binned_rate(times, bin_width)
+   
+    fig = plt.figure()
+    ax1 = fig.add_subplot(411)
+    ax2 = fig.add_subplot(412)
+    ax3 = fig.add_subplot(413)
+    ax4 = fig.add_subplot(414)
+    ax1.spines["right"].set_visible(False)
+    ax1.spines["top"].set_visible(False)
+    ax1.yaxis.set_ticks_position('left')
+    ax1.xaxis.set_ticks_position('bottom')
+    ax2.spines["right"].set_visible(False)
+    ax2.spines["top"].set_visible(False)
+    ax2.yaxis.set_ticks_position('left')
+    ax2.xaxis.set_ticks_position('bottom')
+    ax3.spines["right"].set_visible(False)
+    ax3.spines["top"].set_visible(False)
+    ax3.yaxis.set_ticks_position('left')
+    ax3.xaxis.set_ticks_position('bottom')
+    ax4.spines["right"].set_visible(False)
+    ax4.spines["top"].set_visible(False)
+    ax4.yaxis.set_ticks_position('left')
+    ax4.xaxis.set_ticks_position('bottom')
+
+
+
+    ax1.vlines(times[times < (100000*dt)], ymin=0, ymax=1, color="dodgerblue", lw=1.5)
+    ax1.set_ylabel("spikes", fontsize=10)
+    ax1.set_xlim([1.5, 3.5])
+    ax1.set_ylim([0, 1])
+    ax1.set_yticks([0, 1])
+    set_axis_fontsize(ax1, 10)
+    ax1.set_xticklabels([])
+
+    ax2.plot(time, inst_rate, lw=1.5, label="instantaneous rate")
+    ax2.set_ylabel("firing rate [Hz]", fontsize=10)
+    ax2.legend(fontsize=10)
+    ax2.set_xlim([1.5, 3.5])
+    ax2.set_ylim([0, 300])
+    set_axis_fontsize(ax2, 10)
+    ax2.set_xticklabels([])
+
+    ax3.plot(time, binn_rate, lw=1.5, label="binned rate")
+    ax3.set_ylabel("firing rate [Hz]", fontsize=10)
+    ax3.legend(fontsize=10)
+    ax3.set_xlim([1.5, 3.5])
+    ax3.set_ylim([0, 300])
+    set_axis_fontsize(ax3, 10)
+    ax3.set_xticklabels([])
+
+    ax4.plot(time, conv_rate, lw=1.5, label="convolved rate")
+    ax4.set_xlabel("time [s]", fontsize=10)
+    ax4.set_ylabel("firing rate [Hz]", fontsize=10)
+    ax4.legend(fontsize=10)
+    ax4.set_xlim([1.5, 3.5])
+    ax4.set_ylim([0, 300])
+    set_axis_fontsize(ax4, 10)
+
+    fig.set_size_inches(7.5, 5)
+    fig.subplots_adjust(left=0.1, bottom=0.125, top=0.95, right=0.95, )
+    fig.set_facecolor("white")
+    fig.savefig("firingrates.pdf")
+    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)
+    plot_comparison(spike_times, 0.05, 0.025)
+

pointprocesses/lecture/isimethod.py

@@ -0,0 +1,158 @@
+import matplotlib.pyplot as plt
+import numpy as np
+from IPython import embed
+
+def set_rc():
+    plt.rcParams['xtick.labelsize'] = 8
+    plt.rcParams['ytick.labelsize'] = 8
+    plt.rcParams['xtick.major.size'] = 5
+    plt.rcParams['xtick.minor.size'] = 5
+    plt.rcParams['xtick.major.width'] = 2
+    plt.rcParams['xtick.minor.width'] = 2
+    plt.rcParams['ytick.major.size'] = 5
+    plt.rcParams['ytick.minor.size'] = 5
+    plt.rcParams['ytick.major.width'] = 2
+    plt.rcParams['ytick.minor.width'] = 2
+    plt.rcParams['xtick.direction'] = "out"
+    plt.rcParams['ytick.direction'] = "out"
+
+
+def create_spikes(isi=0.08, duration=0.5):
+    times = np.arange(0., duration, isi)
+    times += np.random.randn(len(times)) * (isi / 2.5)
+    times = np.delete(times, np.nonzero(times < 0))
+    times = np.delete(times, np.nonzero(times > duration))
+    times = np.sort(times)
+    return times
+
+
+def gaussian(sigma, dt):
+    x = np.arange(-4*sigma, 4*sigma, dt)
+    y = np.exp(-0.5 * (x / sigma)**2)/np.sqrt(2*np.pi)/sigma; 
+    return x, y
+
+    
+def setup_axis(spikes_ax, rate_ax):
+    spikes_ax.spines["right"].set_visible(False)
+    spikes_ax.spines["top"].set_visible(False)
+    spikes_ax.yaxis.set_ticks_position('left')
+    spikes_ax.xaxis.set_ticks_position('bottom')
+    spikes_ax.set_yticks([0, 1.0])
+    spikes_ax.set_ylim([0, 1.05])
+    spikes_ax.set_ylabel("spikes", fontsize=9)
+    spikes_ax.text(-0.125, 1.2, "A", transform=spikes_ax.transAxes, size=10)
+
+    rate_ax.spines["right"].set_visible(False)
+    rate_ax.spines["top"].set_visible(False)
+    rate_ax.yaxis.set_ticks_position('left')
+    rate_ax.xaxis.set_ticks_position('bottom')
+    rate_ax.set_xlabel('time[s]', fontsize=9)
+    rate_ax.set_ylabel('firing rate [Hz]', fontsize=9)
+    rate_ax.text(-0.125, 1.15, "B", transform=rate_ax.transAxes, size=10)   
+
+
+def plot_bin_method():
+    dt = 1e-5
+    duration = 0.5
+    
+    spike_times = create_spikes(0.018, duration)
+    t = np.arange(0., duration, dt)
+    
+    bins = np.arange(0, 0.55, 0.05)
+    count, _ = np.histogram(spike_times, bins)
+    
+    plt.xkcd()
+    set_rc()
+    fig = plt.figure()
+    fig.set_size_inches(5., 2.5)
+    fig.set_facecolor('white')
+    spikes = plt.subplot2grid((7,1), (0,0), rowspan=3, colspan=1)
+    rate_ax = plt.subplot2grid((7,1), (3,0), rowspan=4, colspan=1)
+    setup_axis(spikes, rate_ax)
+    spikes.set_ylim([0., 1.25])
+
+    spikes.vlines(spike_times, 0., 1., color="darkblue", lw=1.25)    
+    spikes.vlines(np.hstack((0,bins)), 0., 1.25, color="red", lw=1.5, linestyles='--')
+    for i,c in enumerate(count):
+        spikes.text(bins[i] + bins[1]/2, 1.05, str(c), fontdict={'color':'red'})
+    spikes.set_xlim([0, duration])
+
+    rate = count / 0.05
+    rate_ax.step(bins, np.hstack((rate, rate[-1])), where='post')
+    rate_ax.set_xlim([0., duration])
+    rate_ax.set_ylim([0., 100.])
+    rate_ax.set_yticks(np.arange(0,105,25))
+    fig.tight_layout()
+    fig.savefig("binmethod.pdf")
+    plt.close()
+
+
+def plot_conv_method():
+    dt = 1e-5
+    duration = 0.5
+    spike_times = create_spikes(0.05, duration) 
+    kernel_time, kernel = gaussian(0.02, dt)
+    
+    t = np.arange(0., duration, dt)
+    rate = np.zeros(t.shape)
+    rate[np.asarray(np.round(spike_times/dt), dtype=int)] = 1  
+    rate = np.convolve(rate, kernel, mode='same')
+    rate = np.roll(rate, -1)
+
+    plt.xkcd()
+    set_rc()
+    fig = plt.figure()
+    fig.set_size_inches(5., 2.5)
+    fig.set_facecolor('white')
+    spikes = plt.subplot2grid((7,1), (0,0), rowspan=3, colspan=1)
+    rate_ax = plt.subplot2grid((7,1), (3,0), rowspan=4, colspan=1)
+    setup_axis(spikes, rate_ax)
+    
+    spikes.vlines(spike_times, 0., 1., color="darkblue", lw=1.5, zorder=2)
+    for i in spike_times:
+        spikes.plot(kernel_time + i, kernel/np.max(kernel), color="orange", lw=0.75, zorder=1)
+    spikes.set_xlim([0, duration])
+
+    rate_ax.plot(t, rate, color="darkblue", lw=1, zorder=2)
+    rate_ax.fill_between(t, rate, np.zeros(len(rate)), color="red", alpha=0.5)
+    rate_ax.set_xlim([0, duration])
+    rate_ax.set_ylim([0, 50])
+    rate_ax.set_yticks(np.arange(0,75,25))
+    fig.tight_layout()
+    fig.savefig("convmethod.pdf")
+
+
+def plot_isi_method():
+    spike_times = create_spikes(0.09, 0.5)
+    
+    plt.xkcd()
+    set_rc()
+    fig = plt.figure()
+    fig.set_size_inches(5., 2.5)
+    fig.set_facecolor('white')
+    
+    spikes = plt.subplot2grid((7,1), (0,0), rowspan=3, colspan=1)
+    rate = plt.subplot2grid((7,1), (3,0), rowspan=4, colspan=1)
+    setup_axis(spikes, rate)
+    
+    spikes.vlines(spike_times, 0., 1., color="darkblue", lw=1.25)
+    spike_times = np.hstack((0, spike_times))
+    for i in range(1, len(spike_times)):
+        t_start = spike_times[i-1]
+        t = spike_times[i]
+        spikes.annotate(s='', xy=(t_start, 0.5), xytext=(t,0.5), arrowprops=dict(arrowstyle='<->'), color='red')
+
+    i_rate = 1./np.diff(spike_times)
+
+    rate.step(spike_times, np.hstack((i_rate, i_rate[-1])),color="darkblue", lw=1.25, where="post")
+    rate.set_ylim([0, 75])
+    rate.set_yticks(np.arange(0,100,25))
+    
+    fig.tight_layout()
+    fig.savefig("isimethod.pdf")
+    
+    
+if __name__ == '__main__':
+    plot_isi_method()
+    plot_conv_method()
+    plot_bin_method()

pointprocesses/lecture/pointprocesses.tex

@@ -1,6 +1,35 @@
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-\chapter{\tr{Point processes}{Punktprozesse}}
+\chapter{Analyse von Spiketrains}
+
+\begin{figure}[t]
+  \includegraphics[width=1\textwidth]{rasterexamples}
+  \caption{\label{rasterexamplesfig}Raster-Plot von jeweils 10
+    Realisierungen eines station\"arenen Punktprozesses (homogener
+    Poisson Prozess mit Rate $\lambda=20$\;Hz, links) und eines
+    nicht-station\"aren Punktprozesses (perfect integrate-and-fire
+    Neuron getrieben mit Ohrnstein-Uhlenbeck Rauschen mit
+    Zeitkonstante $\tau=100$\,ms, rechts).}
+\end{figure}
+
+Aktionspotentiale (``Spikes'') sind die Tr\"ager der Information in
+Nervensystemen.  Dabei ist in erster Linie nur der Zeitpunkt des
+Auftretens eines Aktionspotentials von Bedeutung. Die genaue Form
+spielt keine oder nur eine untergeordnete Rolle.
+
+Nach etwas Vorverarbeitung haben elektrophysiologischer Messungen
+deshalb Listen von Spikezeitpunkten als Ergebniss - sogenannte
+``Spiketrains''. Diese Messungen k\"onnen wiederholt werden und es
+ergeben sich mehrere ``trials'' von Spiketrains
+(\figref{rasterexamples}).
+
+Spiketrains sind Zeitpunkte von Ereignissen --- den Aktionspotentialen
+--- und deren Analyse f\"allt daher in das Gebiet der Statistik von
+sogenannten Punktprozessen.
+
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+\section{Punktprozesse}
 
 \begin{figure}[t]
   \texpicture{pointprocessscetchB}
@@ -28,16 +57,6 @@ erzeugt. Zum Beispiel:
   Dynamik des Nervensystems und des Muskelapparates erzeugt.
 \end{itemize}
 
-\begin{figure}[t]
-  \includegraphics[width=1\textwidth]{rasterexamples}
-  \caption{\label{rasterexamplesfig}Raster-Plot von jeweils 10
-    Realisierungen eines station\"arenen Punktprozesses (homogener
-    Poisson Prozess mit Rate $\lambda=20$\;Hz, links) und eines
-    nicht-station\"aren Punktprozesses (perfect integrate-and-fire
-    Neuron getrieben mit Ohrnstein-Uhlenbeck Rauschen mit
-    Zeitkonstante $\tau=100$\,ms, rechts).}
-\end{figure}
-
 F\"ur die Neurowissenschaften ist die Statistik der Punktprozesse
 besonders wichtig, da die Zeitpunkte der Aktionspotentiale als
 zeitlicher Punktprozess betrachtet werden k\"onnen und entscheidend
@@ -116,9 +135,9 @@ Intervalls mit sich selber).
 %   \caption{\label{countstatsfig}Count Statistik.}
 % \end{figure}
 
-Die Anzahl der Ereignisse $n_i$ in Zeifenstern $i$ der
+Die Anzahl der Ereignisse (Spikes) $n_i$ in Zeifenstern $i$ der
-L\"ange $W$ ergeben ganzzahlige, positive Zufallsvariablen die meist durch folgende
+L\"ange $W$ ergeben ganzzahlige, positive Zufallsvariablen, die meist
-Sch\"atzer charakterisiert werden:
+durch folgende Sch\"atzer charakterisiert werden:
 \begin{itemize}
 \item Histogramm der counts $n_i$.
 \item Mittlere Anzahl von Ereignissen: $\mu_N = \langle n \rangle$.
@@ -146,7 +165,7 @@ Insbesondere ist die mittlere Rate der Ereignisse $r$ (``Spikes pro Zeit'', Feue
 \section{Homogener Poisson Prozess}
 F\"ur kontinuierliche Me{\ss}gr\"o{\ss}en ist die Normalverteilung
 u.a. wegen dem Zentralen Grenzwertsatz die Standardverteilung. Eine
-\"ahnliche Rolle spilet bei Punktprozessen der ``Poisson Prozess''.
+\"ahnliche Rolle spielt 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
@@ -192,3 +211,176 @@ Der homogne Poissonprozess hat folgende Eigenschaften:
   \includegraphics[width=0.48\textwidth]{poissoncounthistdist100hz100ms}
   \caption{\label{hompoissoncountfig}Z\"ahlstatistik von Poisson Spikes.}
 \end{figure}
+
+
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+\section{Zeitabh\"angiger Feuerraten}
+
+Bisher haben wir station\"are Spiketrains betrachtet, deren Statistik
+innerhlab der Analysezeit nicht ver\"andert.  Meistens jedoch \"andert
+sich die Statistik der Spiketrains eines Neurons mit der
+Zeit. Z.B. kann ein sensorisches Neuron auf einen Reiz hin mit einer
+erh\"ohten Feuerrate antworten. Im folgenden Stellen wir drei Methoden
+vor, mit denen eine sich zeitlich ver\"andernde Rate des Punktprozesses
+abgesch\"atzt werden kann.
+
+Eine klassische Darstellung zeitabh\"angiger neuronaler Aktivit\"at
+ist das sogenannte Peri Stimulus Zeithistogramm (peri stimulus time
+histogram, PSTH). Es wird der zeitliche Verlauf der Feuerrate $r(t)$
+dargestellt. Die Einheit der Feuerrate ist Hertz, das heisst, die
+Anzahl Aktionspotentiale pro Sekunde. Es verschiedene Methoden diese
+zu bestimmen. Drei solcher Methoden sind in Abbildung \ref{psthfig}
+dargestellt. Alle Methoden haben ihre Berechtigung und ihre Vor- und
+Nachteile. Im folgenden werden die drei Methoden aus Abbildung
+\ref{psthfig} n\"aher erl\"autert.
+
+\begin{figure}[t]
+  \includegraphics[width=\columnwidth]{firingrates}
+  \caption{\textbf{Verschiedene Methoden die zeitabh\"angige Feuerrate
+      zu bestimmen. A)} Rasterplot einer einzelnen neuronalen
+    Antwort. Jeder vertikale Strich notiert den Zeitpunkt eines
+    Aktionspotentials. \textbf{B)} Feurerrate aus der instantanen
+    Feuerrate bestimmt. \textbf{C)} klassisches PSTH mit der Binning
+    Methode. \textbf{D)} Feuerrate durch Faltung mit einem Gauss Kern
+    bestimmt.}\label{psthfig}
+\end{figure}
+
+
+\paragraph{Instantane Feuerrate}
+Ein sehr einfacher Weg, die zeitabh\"angige Feuerrate zu bestimmen ist
+die sogenannte \textit{instantane Feuerrate}. Dabei wird die Feuerrate
+aus dem Kehrwert des \textit{Interspike Intervalls}, der Zeit zwischen
+zwei aufeinander folgenden Aktionspotentialen (Abbildung \ref{instrate}
+A), bestimmt. Die abgesch\"atzte Feuerrate (Abbildung \ref{instrate} B)
+ist g\"ultig f\"ur das gesammte Interspike Intervall
+\ref{instrate}. Diese Methode hat den Vorteil, dass sie sehr einfach zu
+berechnen ist und keine Annahme \"uber eine relevante Zeitskala (der
+Kodierung oder des Auslesemechanismus der postsynaptischen Zelle)
+macht. $r(t)$ ist allerdings keine kontinuierliche Funktion, die
+Spr\"unge in der Feuerrate k\"onnen f\"ur manche Analysen nachteilig
+sein. Des Weiteren ist die Feuerrate nie null, auch wenn lange keine
+Aktionspotentiale generiert wurden.
+
+\begin{figure}[!htb]
+  \includegraphics[width=\columnwidth]{isimethod}
+  \caption{\textbf{Bestimmung des zeitabh\"angigen Feuerrate aus dem
+      Interspike Interval. A)} Skizze eines Rasterplots einer
+    einzelnen neuronalen Antwort. Jeder vertikale Strich notiert den
+    Zeitpunkt eines Aktionspotentials. Die Pfeile zwischen
+    aufeinanderfolgenden Aktionspotentialen illustrieren das
+    Interspike Interval. \textbf{B)} Der Kehrwert des Interspike
+    Intervalls ist die Feuerrate.}\label{instrate}
+\end{figure}
+
+
+\paragraph{Binning Methode}
+Bei der Binning Methode wird die Zeitachse in gleichm\"aßige
+Abschnitte (Bins) eingeteilt und die Anzahl Aktionspotentiale, die in
+die jeweiligen Bins fallen, gez\"ahlt (Abbildung \ref{binpsth} A). Um
+diese Z\"ahlungen in die Feuerrate umzurechnen muss noch mit der
+Binweite normiert werden. Die bestimmte Feuerrate gilt f\"ur das
+gesamte Bin (Abbildung \ref{binpsth} B). \textbf{Tipp:} Um die Anzahl
+Spikes pro Bin zu berechnen kann die \code{hist} Funktion benutzt
+werden. Das so berechnete PSTH hat wiederum eine stufige Form, die von
+der Wahl der Binweite anh\"angt. Die Binweite bestimmt die zeitliche
+Aufl\"osung der Darstellung. \"Anderungen in der Feuerrate, die
+innerhalb eines Bins vorkommen k\"onnen nicht aufgl\"ost werden. Die
+Wahl der Binweite stellt somit eine Annahme \"uber die relevante
+Zeitskala der Verarbeitung dar. Auch hier ist $r(t)$ keine
+koninuierliche Funktion.
+
+\begin{figure}[h!]
+  \includegraphics[width=\columnwidth]{binmethod}
+  \caption{\textbf{Bestimmung des PSTH mit der Binning Methode. A)}
+    Skizze eines Rasterplots einer einzelnen neuronalen Antwort. Jeder
+    vertikale Strich notiert den Zeitpunkt eines
+    Aktionspotentials. Die roten gestrichelten Linien stellen die
+    Grenzen der Bins dar und die Zahlen geben den Spike Count pro Bin
+    an. \textbf{B)} Die Feuerrate erh\"alt man indem das
+    Zeithistogramm mit der Binweite normiert.}\label{binpsth}
+\end{figure}
+
+
+\paragraph{Faltungsmethode}
+Bei der Faltungsmethode geht wird etwas anders vorgegangen um die
+Feuerrate zeitaufgel\"ost zu berechnen. Die Aktionspotentialfolge wird
+zun\"achst ``bin\"ar'' dargestellt. Das heisst, dass eine Antwort als
+(Zeit-)Vektor dargestellt wird, in welchem die Zeitpunkte der
+Aktionspotentiale als 1 notiert werden. Alle anderen Elemente des
+Vektors sind 0. Anschlie{\ss}end wir dieser bin\"are Spiketrain mit
+einem Gauss Kern bestimmter Breite gefaltet.
+ 
+\[r(t) = \int_{-\infty}^{\infty}d\tau \omega(\tau)\rho(t-\tau) \],
+wobei $\omega(\tau)$ der Filterkern und $\rho(t)$ die bin\"are Antwort
+ist. Bildlich geprochen wird jede 1 in $rho(t)$ durch den Filterkern
+ersetzt (Abbildung \ref{convrate} A). Wenn der Kern richtig normiert
+wurde (Integral 1), ergibt sich die Feuerrate direkt aus der
+\"Uberlagerung der Kerne (Abb. \ref{convrate} B). Die Faltungsmethode
+f\"uhrt, anders als die anderen Methoden, zu einer kontinuierlichen
+Funktion was f\"ur spektrale Analysen von Vorteil sein kann. Die Wahl
+der Kernbreite bestimmt, \"ahnlich zur Binweite, die zeitliche
+Aufl\"osung von $r(t)$. Man macht also eine Annahme \"uber die
+relevante Zeitskala.
+
+\begin{figure}[h!]
+  \includegraphics[width=\columnwidth]{convmethod}
+  \caption{\textbf{Schematische Darstellung der Faltungsmethode. A)}
+    Rasterplot einer einzelnen neuronalen Antwort. Jeder vertikale
+    Strich notiert den Zeitpunkt eines Aktionspotentials. In der
+    Faltung werden die mit einer 1 notierten Aktionspotential durch
+    den Faltungskern ersetzt. \textbf{B)} Bei korrekter Normierung des
+    Kerns ergibt sich die Feuerrate direkt aus der \"Uberlagerung der
+    Kerne.}\label{convrate}
+\end{figure}
+
+
+
+\section{Spike triggered Average}
+Die graphischer Darstellung der Feuerrate allein reicht nicht aus um
+den Zusammenhang zwischen neuronaler Antwort und einem Stimulus zu
+analysieren.  Eine Methode mehr \"uber diesen Zusammenhang zu erfahren
+ist der Spike triggered average (STA). Der STA ist der mittlere
+Stimulus, der zu einem Aktionspotential in der neuronalen Antwort
+f\"uhrt.
+
+\begin{equation}
+  STA(\tau) = \frac{1}{\langle n \rangle} \left\langle  \displaystyle\sum_{i=1}^{n}{s(t_i - \tau)} \right\rangle
+\end{equation}
+
+Der STA l\"a{\ss}t sich relativ einfach berechnen, indem aus dem
+Stimulus f\"ur jeden beobachteten Spike ein entsprechender Abschnitt
+ausgeschnitten wird und diese dann mittelt (Abbildung
+\ref{stafig}). 
+
+\begin{figure}
+  \includegraphics[width=0.5\columnwidth]{sta}
+  \caption{\textbf{Spike Triggered Average eines P-Typ
+      Elektrorezeptors.} Der Rezeptor wurde mit einem ``white-noise''
+    Stimulus getrieben. Zeitpunkt 0 ist der Zeitpunkt des beobachteten
+    Aktionspotentials.}\label{stafig}
+\end{figure}
+
+Aus dem STA k\"onnen verschiedene Informationen \"uber den
+Zusammenhang zwischen Stimulus und neuronaler Antwort gewonnen
+werden. Die Breite des STA repr\"asentiert die zeitliche Pr\"azision,
+mit der Stimulus und Antwort zusammenh\"angen und wie weit das neuron
+zeitlich integriert. Die Amplitude des STA (gegeben in der gleichen
+Einheit, wie der Stimulus) deutet auf die Empfindlichkeit des Neurons
+bez\"uglich des Stimulus hin. Eine hohe Amplitude besagt, dass es
+einen starken Stimulus ben\"otigt um ein Aktionspotential
+hervorzurufen. Aus dem zeitlichen Versatz des STA kann die Zeit
+abgelesen werden, die das System braucht um auf den Stimulus zu
+antworten.
+
+Der STA kann auch dazu benutzt werden, aus den Antworten der Zelle den
+Stimulus zu rekonstruieren (Abbildung
+\ref{reverse_reconstruct_fig}). Bei der \textit{invertierten
+  Rekonstruktion} wird die Zellantwort mit dem STA gefaltet.
+
+\begin{figure}
+  \includegraphics[width=\columnwidth]{reconstruction}
+  \caption{\textbf{Rekonstruktion des Stimulus mittels STA.} Die
+    Zellantwort wird mit dem STA gefaltet um eine Rekonstruktion des
+    Stimulus zu erhalten.}\label{reverse_reconstruct_fig}
+\end{figure}

pointprocesses/lecture/sta.py

@@ -0,0 +1,88 @@
+import numpy as np
+import matplotlib.pyplot as plt
+import scipy.io as scio
+from IPython import embed
+
+
+def plot_sta(times, stim, dt, t_min=-0.1, t_max=.1):
+    count = 0
+    sta = np.zeros((abs(t_min) + abs(t_max))/dt)
+    time = np.arange(t_min, t_max, dt)
+    if len(stim.shape) > 1 and stim.shape[1] > 1:
+        stim = stim[:,1]
+    for i in range(len(times[0])):
+        times = np.squeeze(spike_times[0][i])
+        for t in times:
+            if (int((t + t_min)/dt) < 0) or ((t + t_max)/dt > len(stim)):
+                continue;
+            
+            min_index = int(np.round((t+t_min)/dt))
+            max_index = int(np.round((t+t_max)/dt))
+            snippet = np.squeeze(stim[ min_index : max_index])            
+            sta += snippet
+            count += 1
+    sta /= count
+    
+    fig = plt.figure()
+    fig.set_size_inches(5, 5)
+    fig.subplots_adjust(left=0.15, bottom=0.125, top=0.95, right=0.95, )
+    fig.set_facecolor("white")
+    
+    ax = fig.add_subplot(111)
+    ax.plot(time, sta, color="darkblue", lw=1)
+    ax.set_xlabel("time [s]")
+    ax.set_ylabel("stimulus")
+    ax.xaxis.grid('off')
+    ax.spines["right"].set_visible(False)
+    ax.spines["top"].set_visible(False)
+    ax.yaxis.set_ticks_position('left')
+    ax.xaxis.set_ticks_position('bottom')
+
+    ylim = ax.get_ylim()
+    xlim = ax.get_xlim()
+    ax.plot(list(xlim), [0., 0.], zorder=1, color='darkgray', ls='--')
+    ax.plot([0., 0.], list(ylim), zorder=1, color='darkgray', ls='--')
+    ax.set_xlim(list(xlim))
+    ax.set_ylim(list(ylim))
+    fig.savefig("sta.pdf")
+    plt.close()
+    return sta
+
+
+def reconstruct_stimulus(spike_times, sta, stimulus, t_max=30., dt=1e-4):
+    s_est = np.zeros((spike_times.shape[1], len(stimulus)))
+    for i in range(10):
+        times = np.squeeze(spike_times[0][i])
+        indices = np.asarray((np.round(times/dt)), dtype=int)
+        y = np.zeros(len(stimulus))
+        y[indices] = 1
+        s_est[i, :] = np.convolve(y, sta, mode='same')
+
+    plt.plot(np.arange(0, t_max, dt), stimulus[:,1], label='stimulus', color='darkblue', lw=2.)
+    plt.plot(np.arange(0, t_max, dt), np.mean(s_est, axis=0), label='reconstruction', color='gray', lw=1.5)
+    plt.xlabel('time[s]')
+    plt.ylabel('stimulus')
+    plt.xlim([0.0, 0.25])
+    plt.ylim([-1., 1.])
+    plt.legend()
+    plt.plot([0.0, 0.25], [0., 0.], color="darkgray", lw=1.5, ls='--', zorder=1)
+    plt.gca().spines["right"].set_visible(False)
+    plt.gca().spines["top"].set_visible(False)
+    plt.gca().yaxis.set_ticks_position('left')
+    plt.gca().xaxis.set_ticks_position('bottom')
+
+    fig = plt.gcf()
+    fig.set_size_inches(7.5, 5)
+    fig.subplots_adjust(left=0.15, bottom=0.125, top=0.95, right=0.95, )
+    fig.set_facecolor("white")
+    fig.savefig('reconstruction.pdf')
+    plt.close()
+
+
+if __name__ == "__main__":
+    punit_data =  scio.loadmat('p-unit_spike_times.mat')
+    punit_stim = scio.loadmat('p-unit_stimulus.mat')
+    spike_times = punit_data["spike_times"]
+    stimulus = punit_stim["stimulus"]
+    sta = plot_sta(spike_times, stimulus, 5e-5, -0.05, 0.05)
+    reconstruct_stimulus(spike_times, sta, stimulus, 10, 5e-5)

regression/lecture/Makefile

@@ -1,22 +1,29 @@
 BASENAME=regression
+
 PYFILES=$(wildcard *.py)
 PYPDFFILES=$(PYFILES:.py=.pdf)
 
-pdf : $(BASENAME).pdf $(PYPDFFILES)
+all : pdf 
+
+# script:
+pdf : $(BASENAME)-chapter.pdf
 
-$(BASENAME).pdf : $(BASENAME)-chapter.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).aux $(BASENAME).log $(BASENAME).out $(BASENAME).toc $(BASENAME).nav $(BASENAME).snm $(BASENAME).vrb
+	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).pdf
+	rm -f $(BASENAME)-chapter.pdf
 
-watch :
+watchpdf :
 	while true; do ! make -q pdf && make pdf; sleep 0.5; done
 
 

scientificcomputing-script.tex

@@ -19,10 +19,13 @@
 \lstset{inputpath=programming/code}
 \include{programming/lectures/programming}
 
-\chapter{Graphische Darstellung}
 
-\graphicspath{{designpattern/lecture}{designpattern/lecture/figures}}
+\graphicspath{{plotting/lecture/}{plotting/lecture/images/}}
-\lstset{inputpath=designpattern/code/}
+\lstset{inputpath=plotting/code/}
+\include{plotting/lecture/plotting}
+
+\graphicspath{{designpattern/lecture/}{designpattern/lecture/figures/}}
+\lstset{inputpath=designpattern/code}
 \include{designpattern/lecture/designpattern}
 
 \chapter{Cheat-Sheet}
@@ -52,8 +55,4 @@
 \renewcommand{\texinputpath}{pointprocesses/lecture/}
 \include{pointprocesses/lecture/pointprocesses}
 
-\graphicspath{{spike_trains/lecture/images/}}
-\lstset{inputpath=spike_trains/code/}
-\include{spike_trains/lecture/psth_sta}
-
 \end{document}

spike_trains/lecture/psth_sta.tex

@@ -20,6 +20,8 @@ werden. Dar\"uber hinaus wird versucht die Beziehung zwischen der
 zeitabh\"angigen neuronalen Antwort und dem zugrundeliegenden
 Stimulus zu analysieren.
 
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%% The following moved to the pointprocesses chapter! 
 \section{Darstellung der zeitabh\"angigen Feuerrate}
 
 Eine klassische Darstellung zeitabh\"angiger neuronaler Aktivit\"at