From a945268beeb8d2faa68578774bb2c77e0ea5dbc2 Mon Sep 17 00:00:00 2001
From: Jan Benda <jan.benda@uni-tuebingen.de>
Date: Sat, 21 Nov 2015 23:09:07 +0100
Subject: [PATCH] Improved point processes. Added index.
---
Makefile | 4 +
header.tex | 18 ++--
pointprocesses/code/counthist.m | 29 ++---
pointprocesses/code/counthistpoisson.m | 55 ++++++++++
pointprocesses/code/hompoissonspikes.m | 24 +++++
pointprocesses/code/isi_hist.m | 8 +-
pointprocesses/code/isis.m | 2 +-
pointprocesses/code/isiserialcorr.m | 23 ++--
pointprocesses/code/plot_isi_hist.m | 16 ++-
pointprocesses/lecture/isimethod.py | 124 ++++++++++++----------
pointprocesses/lecture/pointprocesses.tex | 122 +++++++++++----------
scientificcomputing-script.tex | 5 +
12 files changed, 272 insertions(+), 158 deletions(-)
create mode 100644 pointprocesses/code/counthistpoisson.m
create mode 100644 pointprocesses/code/hompoissonspikes.m
diff --git a/Makefile b/Makefile
index 2038f42..16301a5 100644
--- a/Makefile
+++ b/Makefile
@@ -11,6 +11,10 @@ chapters :
$(BASENAME).pdf : $(BASENAME).tex header.tex $(SUBTEXS)
export TEXMFOUTPUT=.; pdflatex -interaction=scrollmode $< | tee /dev/stderr | fgrep -q "Rerun to get cross-references right" && pdflatex -interaction=scrollmode $< || true
+index :
+ splitindex $(BASENAME).idx
+ export TEXMFOUTPUT=.; pdflatex -interaction=scrollmode $(BASENAME).tex
+
clean :
rm -f *~ $(BASENAME).aux $(BASENAME).log $(BASENAME).out $(BASENAME).toc
for sd in $(SUBDIRS); do $(MAKE) -C $$sd/lecture clean; done
diff --git a/header.tex b/header.tex
index d636561..d067f4c 100644
--- a/header.tex
+++ b/header.tex
@@ -23,6 +23,12 @@
\setcounter{secnumdepth}{1}
\setcounter{tocdepth}{1}
+%%%%% index %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+\usepackage[makeindex]{splitidx}
+\makeindex
+\newindex[Fachbegriffe]{term}
+\newindex[Code]{code}
+
%%%%% units %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\usepackage[mediumspace,mediumqspace,Gray]{SIunits} % \ohm, \micro
@@ -32,7 +38,7 @@
%%%%%%%%%%%%% colors %%%%%%%%%%%%%%%%
\usepackage{xcolor}
-\definecolor{myyellow}{HTML}{FFCC00}
+\definecolor{myyellow}{HTML}{FFEE00}
\definecolor{myblue}{HTML}{0055DD}
\colorlet{codeback}{myblue!10}
@@ -183,9 +189,9 @@
\newcommand{\koZ}{\mathds{C}}
%%%%% english, german, code and file terms: %%%%%%%%%%%%%%%
-\newcommand{\enterm}[1]{``#1''}
-\newcommand{\determ}[1]{\textit{#1}}
-\newcommand{\codeterm}[1]{\textit{#1}}
+\newcommand{\enterm}[1]{``#1''\sindex[term]{#1}}
+\newcommand{\determ}[1]{\textit{#1}\sindex[term]{#1}}
+\newcommand{\codeterm}[1]{\textit{#1}\sindex[term]{#1}}
\newcommand{\file}[1]{\texttt{#1}}
%%%%% key-shortcuts %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
@@ -196,9 +202,9 @@
%%%%% code/matlab commands: %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\usepackage{textcomp}
-\newcommand{\code}[1]{\setlength{\fboxsep}{0.5ex}\colorbox{codeback}{\texttt{#1}}}
+\newcommand{\code}[1]{\setlength{\fboxsep}{0.5ex}\colorbox{codeback}{\texttt{#1}}\sindex[code]{#1}}
\newcommand{\matlab}{MATLAB$^{\copyright}$}
-\newcommand{\matlabfun}[1]{(\tr{\matlab{}-function}{\matlab-Funktion} \setlength{\fboxsep}{0.5ex}\colorbox{codeback}{\texttt{#1}})}
+\newcommand{\matlabfun}[1]{(\tr{\matlab{}-function}{\matlab-Funktion} \setlength{\fboxsep}{0.5ex}\colorbox{codeback}{\texttt{#1}})\sindex[code]{#1}}
%%%%% exercises environment: %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% usage:
diff --git a/pointprocesses/code/counthist.m b/pointprocesses/code/counthist.m
index 9d9889a..6885b1a 100644
--- a/pointprocesses/code/counthist.m
+++ b/pointprocesses/code/counthist.m
@@ -1,5 +1,5 @@
function [counts, bins] = counthist(spikes, w)
-% computes count histogram and compare with Poisson distribution
+% Compute and plot histogram of spike counts.
%
% [counts, bins] = counthist(spikes, w)
%
@@ -19,36 +19,25 @@ function [counts, bins] = counthist(spikes, w)
times = spikes{k};
% alternative 1: count the number of spikes in each window:
% for tk = 0:w:tmax-w
-% nn = sum( ( times >= tk ) & ( times < tk+w ) );
-% %nn = length( find( ( times >= tk ) & ( times < tk+w ) ) );
-% n = [ n nn ];
+% nn = sum((times >= tk) & (times < tk+w));
+% %nn = length(find((times >= tk) & (times < tk+w)));
+% n = [n nn];
% end
-% alternative 2: use the hist function to do that!
+% alternative 2: use the hist() function to do that!
tbins = 0.5*w:w:tmax-0.5*w;
nn = hist(times, tbins);
- n = [ n nn ];
- % the rate of the spikes:
- rate = (length(times)-1)/(times(end) - times(1));
- r = [ r rate ];
+ n = [n nn];
end
% histogram of spike counts:
- maxn = max( n );
- [counts, bins ] = hist( n, 0:1:maxn+10 );
+ maxn = max(n);
+ [counts, bins] = hist(n, 0:1:maxn+10);
% normalize to probabilities:
- counts = counts / sum( counts );
+ counts = counts / sum(counts);
% plot:
if nargout == 0
bar( bins, counts );
- hold on;
- % Poisson distribution:
- rate = mean( r );
- x = 0:1:maxn+10;
- a = rate*w;
- y = a.^x.*exp(-a)./factorial(x);
- plot( x, y, 'r', 'LineWidth', 3 );
- hold off;
xlabel( 'counts k' );
ylabel( 'P(k)' );
end
diff --git a/pointprocesses/code/counthistpoisson.m b/pointprocesses/code/counthistpoisson.m
new file mode 100644
index 0000000..9d9889a
--- /dev/null
+++ b/pointprocesses/code/counthistpoisson.m
@@ -0,0 +1,55 @@
+function [counts, bins] = counthist(spikes, w)
+% computes count histogram and compare with Poisson distribution
+%
+% [counts, bins] = counthist(spikes, w)
+%
+% Arguments:
+% spikes: a cell array of vectors of spike times in seconds
+% w: observation window duration in seconds for computing the counts
+%
+% Returns:
+% counts: the histogram of counts normalized to probabilities
+% bins: the bin centers for the histogram
+
+ % collect spike counts:
+ tmax = spikes{1}(end);
+ n = [];
+ r = [];
+ for k = 1:length(spikes)
+ times = spikes{k};
+% alternative 1: count the number of spikes in each window:
+% for tk = 0:w:tmax-w
+% nn = sum( ( times >= tk ) & ( times < tk+w ) );
+% %nn = length( find( ( times >= tk ) & ( times < tk+w ) ) );
+% n = [ n nn ];
+% end
+% alternative 2: use the hist function to do that!
+ tbins = 0.5*w:w:tmax-0.5*w;
+ nn = hist(times, tbins);
+ n = [ n nn ];
+ % the rate of the spikes:
+ rate = (length(times)-1)/(times(end) - times(1));
+ r = [ r rate ];
+ end
+
+ % histogram of spike counts:
+ maxn = max( n );
+ [counts, bins ] = hist( n, 0:1:maxn+10 );
+ % normalize to probabilities:
+ counts = counts / sum( counts );
+
+ % plot:
+ if nargout == 0
+ bar( bins, counts );
+ hold on;
+ % Poisson distribution:
+ rate = mean( r );
+ x = 0:1:maxn+10;
+ a = rate*w;
+ y = a.^x.*exp(-a)./factorial(x);
+ plot( x, y, 'r', 'LineWidth', 3 );
+ hold off;
+ xlabel( 'counts k' );
+ ylabel( 'P(k)' );
+ end
+end
diff --git a/pointprocesses/code/hompoissonspikes.m b/pointprocesses/code/hompoissonspikes.m
new file mode 100644
index 0000000..ff14a1e
--- /dev/null
+++ b/pointprocesses/code/hompoissonspikes.m
@@ -0,0 +1,24 @@
+function spikes = hompoissonspikes(rate, trials, tmax)
+% Generate spike times of a homogeneous poisson process
+% using the exponential interspike interval distribution.
+%
+% spikes = hompoissonspikes(rate, trials, tmax)
+%
+% Arguments:
+% rate: the rate of the Poisson process in Hertz
+% trials: number of trials that should be generated
+% tmax: the duration of each trial in seconds
+%
+% Returns:
+% spikes: a cell array of vectors of spike times in seconds
+
+ spikes = cell(trials, 1);
+ mu = 1.0/rate;
+ nintervals = 2*round(tmax/mu);
+ for k=1:trials
+ % exponential random numbers:
+ intervals = random('Exponential', nintervals, 1);
+ times = cumsum(intervals);
+ spikes{k} = times(times<=tmax);
+ end
+end
diff --git a/pointprocesses/code/isi_hist.m b/pointprocesses/code/isi_hist.m
index 6194948..f529f8c 100644
--- a/pointprocesses/code/isi_hist.m
+++ b/pointprocesses/code/isi_hist.m
@@ -9,14 +9,14 @@ function [pdf, centers] = isi_hist(isis, binwidth)
%
% Returns:
% pdf: vector with probability density of interspike intervals in Hz
-% centers: vector with corresponding centers of interspikeintervalls in seconds
+% centers: vector with centers of interspikeintervalls in seconds
if nargin < 2
% compute good binwidth:
- nperbin = 200; % average number of data points per bin
- bins = length( isis )/nperbin; % number of bins
+ nperbin = 200; % average number of data points per bin
+ bins = length( isis )/nperbin; % number of bins
binwidth = max( isis )/bins;
- if binwidth < 5e-4 % half a millisecond
+ if binwidth < 5e-4 % half a millisecond
binwidth = 5e-4;
end
end
diff --git a/pointprocesses/code/isis.m b/pointprocesses/code/isis.m
index e693ea8..869cb60 100644
--- a/pointprocesses/code/isis.m
+++ b/pointprocesses/code/isis.m
@@ -9,7 +9,7 @@ function isivec = isis( spikes )
for k = 1:length(spikes)
difftimes = diff( spikes{k} );
% difftimes(:) ensures a column vector
- % regardless of the type of vector spikes{k}
+ % regardless of the type of vector in spikes{k}
isivec = [ isivec; difftimes(:) ];
end
end
diff --git a/pointprocesses/code/isiserialcorr.m b/pointprocesses/code/isiserialcorr.m
index 8b7aca9..17df321 100644
--- a/pointprocesses/code/isiserialcorr.m
+++ b/pointprocesses/code/isiserialcorr.m
@@ -1,34 +1,35 @@
-function isicorr = isiserialcorr(isis, maxlag)
+function isicorr = isiserialcorr(isivec, maxlag)
% serial correlation of interspike intervals
%
-% isicorr = isiserialcorr(isis, maxlag)
+% isicorr = isiserialcorr(isivec, maxlag)
%
% Arguments:
-% isis: vector of interspike intervals in seconds
+% isivec: vector of interspike intervals in seconds
% maxlag: the maximum lag in seconds
%
% Returns:
% isicorr: vector with the serial correlations for lag 0 to maxlag
lags = 0:maxlag;
- isicorr = zeros( size( lags ) );
+ isicorr = zeros(size(lags));
for k = 1:length(lags)
lag = lags(k);
- if length( isis ) > lag+10 % ensure "enough" data
+ if length(isivec) > lag+10 % ensure "enough" data
% NOTE: the arguments to corr must be column vectors!
- % We insure this in the isis() function that generats the isis.
- isicorr(k) = corr( isis(1:end-lag), isis(lag+1:end) );
+ % We insure this in the isis() function that
+ % generates the isivec.
+ isicorr(k) = corr(isivec(1:end-lag), isivec(lag+1:end));
end
end
if nargout == 0
% plot:
- plot( lags, isicorr, '-b' );
+ plot(lags, isicorr, '-b');
hold on;
- scatter( lags, isicorr, 100.0, 'b', 'filled' );
+ scatter(lags, isicorr, 100.0, 'b', 'filled');
hold off;
- xlabel( 'Lag k' )
- ylabel( '\rho_k')
+ xlabel('Lag k')
+ ylabel('\rho_k')
end
end
diff --git a/pointprocesses/code/plot_isi_hist.m b/pointprocesses/code/plot_isi_hist.m
index d49b5b2..9fb6bf9 100644
--- a/pointprocesses/code/plot_isi_hist.m
+++ b/pointprocesses/code/plot_isi_hist.m
@@ -7,22 +7,28 @@ function plot_isi_hist(isis, binwidth)
% isis: vector of interspike intervals in seconds
% binwidth: optional width in seconds to be used for the isi bins
+ % compute normalized histogram:
if nargin < 2
[pdf, centers] = isi_hist(isis);
else
[pdf, centers] = isi_hist(isis, binwidth);
end
- bar(1000.0*centers, nelements); % milliseconds on x-axis
+ % plot:
+ bar(1000.0*centers, pdf); % milliseconds on x-axis
xlabel('ISI [ms]')
ylabel('p(ISI) [1/s]')
% annotation:
misi = mean(isis);
sdisi = std(isis);
disi = sdisi^2.0/2.0/misi^3;
- text(0.95, 0.8, sprintf('mean=%.1f ms', 1000.0*misi), 'Units', 'normalized', 'HorizontalAlignment', 'right')
- text(0.95, 0.7, sprintf('std=%.1f ms', 1000.0*sdisi), 'Units', 'normalized', 'HorizontalAlignment', 'right')
- text(0.95, 0.6, sprintf('CV=%.2f', sdisi/misi), 'Units', 'normalized', 'HorizontalAlignment', 'right')
- %text(0.5, 0.3, sprintf('D=%.1f Hz', disi), 'Units', 'normalized')
+ text(0.95, 0.8, sprintf('mean=%.1f ms', 1000.0*misi), ...
+ 'Units', 'normalized', 'HorizontalAlignment', 'right')
+ text(0.95, 0.7, sprintf('std=%.1f ms', 1000.0*sdisi), ...
+ 'Units', 'normalized', 'HorizontalAlignment', 'right')
+ text(0.95, 0.6, sprintf('CV=%.2f', sdisi/misi), ...
+ 'Units', 'normalized', 'HorizontalAlignment', 'right')
+ %text(0.95, 0.5, sprintf('D=%.1f Hz', disi), ...
+ % 'Units', 'normalized', 'HorizontalAlignment', 'right')
end
diff --git a/pointprocesses/lecture/isimethod.py b/pointprocesses/lecture/isimethod.py
index 07578f1..b0baa64 100644
--- a/pointprocesses/lecture/isimethod.py
+++ b/pointprocesses/lecture/isimethod.py
@@ -2,9 +2,9 @@ import matplotlib.pyplot as plt
import numpy as np
from IPython import embed
+figsize=(6,3)
+
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
@@ -17,13 +17,17 @@ def set_rc():
plt.rcParams['ytick.direction'] = "out"
-def create_spikes(isi=0.08, duration=0.5, seed=57111):
- times = np.arange(0., duration, isi)
+def create_spikes(nspikes=11, duration=0.5, seed=1000):
rng = np.random.RandomState(seed)
- times += rng.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)
+ x = np.linspace(0.0, 1.0, nspikes)
+ # double gaussian rate profile:
+ rate = np.exp(-0.5*((x-0.35)/0.25)**2.0)
+ rate += 1.*np.exp(-0.5*((x-0.9)/0.05)**2.0)
+ isis = 1.0/rate
+ isis += rng.randn(len(isis))*0.2
+ times = np.cumsum(isis)
+ times *= 1.05*duration/times[-1]
+ times += 0.01
return times
@@ -34,33 +38,38 @@ def gaussian(sigma, dt):
def setup_axis(spikes_ax, rate_ax):
+ spikes_ax.spines["left"].set_visible(False)
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)
- spikes_ax.set_xlim([0, 0.5])
- spikes_ax.set_xticklabels(np.arange(0., 600, 100))
+ spikes_ax.set_yticks([])
+ spikes_ax.set_ylim(-0.2, 1.0)
+ #spikes_ax.set_ylabel("Spikes")
+ spikes_ax.text(-0.1, 0.5, "Spikes", transform=spikes_ax.transAxes, rotation='vertical', va='center')
+ #spikes_ax.text(-0.125, 1.2, "A", transform=spikes_ax.transAxes)
+ spikes_ax.set_xlim(-1, 500)
+ #spikes_ax.set_xticklabels(np.arange(0., 600, 100))
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[ms]', 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)
- rate_ax.set_xlim([0, 0.5])
- rate_ax.set_xticklabels(np.arange(0., 600, 100))
+ rate_ax.set_xlabel('Time [ms]')
+ #rate_ax.set_ylabel('Firing rate [Hz]')
+ rate_ax.text(-0.1, 0.5, "Rate [Hz]", transform=rate_ax.transAxes, rotation='vertical', va='center')
+ #rate_ax.text(-0.125, 1.15, "B", transform=rate_ax.transAxes)
+ rate_ax.set_xlim(0, 500)
+ #rate_ax.set_xticklabels(np.arange(0., 600, 100))
+ rate_ax.set_ylim(0, 60)
+ rate_ax.set_yticks(np.arange(0,65,20))
def plot_bin_method():
dt = 1e-5
duration = 0.5
- spike_times = create_spikes(0.018, duration)
+ spike_times = create_spikes()
t = np.arange(0., duration, dt)
bins = np.arange(0, 0.55, 0.05)
@@ -69,24 +78,25 @@ def plot_bin_method():
plt.xkcd()
set_rc()
fig = plt.figure()
- fig.set_size_inches(5., 2.5)
+ fig.set_size_inches(*figsize)
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 ti in spike_times:
+ ti *= 1000.0
+ spikes.plot([ti, ti], [0., 1.], '-b', lw=2)
+
+ #spikes.vlines(1000.0*spike_times, 0., 1., color="darkblue", lw=1.25)
+ for tb in 1000.0*bins :
+ spikes.plot([tb, tb], [-2.0, 0.75], '-', color="#777777", lw=1, clip_on=False)
+ #spikes.vlines(1000.0*np.hstack((0,bins)), -2.0, 1.25, color="#777777", lw=1, linestyles='-', clip_on=False)
for i,c in enumerate(count):
- spikes.text(bins[i] + bins[1]/2, 1.05, str(c), fontdict={'color':'red', 'size':9})
- spikes.set_xlim([0, duration])
+ spikes.text(1000.0*(bins[i]+0.5*bins[1]), 1.1, str(c), color='#CC0000', ha='center')
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))
+ rate_ax.step(1000.0*bins, np.hstack((rate, rate[-1])), color='#FF9900', where='post')
fig.tight_layout()
fig.savefig("binmethod.pdf")
plt.close()
@@ -95,66 +105,66 @@ def plot_bin_method():
def plot_conv_method():
dt = 1e-5
duration = 0.5
- spike_times = create_spikes(0.05, duration)
- kernel_time, kernel = gaussian(0.02, dt)
+ spike_times = create_spikes()
+ kernel_time, kernel = gaussian(0.015, 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.asarray(np.round(spike_times[:-1]/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_size_inches(*figsize)
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))
+ #spikes.vlines(1000.0*spike_times, 0., 1., color="darkblue", lw=1.5, zorder=2)
+ for ti in spike_times:
+ ti *= 1000.0
+ spikes.plot([ti, ti], [0., 1.], '-b', lw=2)
+ spikes.plot(1000*kernel_time + ti, kernel/np.max(kernel), color='#cc0000', lw=1, zorder=1)
+
+ rate_ax.plot(1000.0*t, rate, color='#FF9900', lw=2, zorder=2)
+ rate_ax.fill_between(1000.0*t, rate, np.zeros(len(rate)), color='#FFFF66')
+ #rate_ax.fill_between(t, rate, np.zeros(len(rate)), color="red", alpha=0.5)
+ #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.055, 0.5, 1000)
+ spike_times = create_spikes()
plt.xkcd()
set_rc()
fig = plt.figure()
- fig.set_size_inches(5., 2.5)
+ fig.set_size_inches(*figsize)
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))
+ spike_times = np.hstack((0.005, spike_times))
+ #spikes.vlines(1000*spike_times, 0., 1., color="blue", lw=2)
for i in range(1, len(spike_times)):
- t_start = spike_times[i-1]
- t = spike_times[i]
+ t_start = 1000*spike_times[i-1]
+ t = 1000*spike_times[i]
+ spikes.plot([t_start, t_start], [0., 1.], '-b', lw=2)
spikes.annotate(s='', xy=(t_start, 0.5), xytext=(t,0.5), arrowprops=dict(arrowstyle='<->'), color='red')
- spikes.text(t_start+0.01, 0.75,
- "{0:.1f}".format((t - t_start)*1000),
- fontdict={'color':'red', 'size':7})
+ spikes.text(0.5*(t_start+t), 0.75,
+ "{0:.0f}".format((t - t_start)),
+ color='#CC0000', ha='center')
+ #spike_times = np.hstack((0, spike_times))
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, 50])
- rate.set_yticks(np.arange(0,75,25))
+ rate.step(1000*spike_times, np.hstack((i_rate, i_rate[-1])),color='#FF9900', lw=2, where="post")
fig.tight_layout()
fig.savefig("isimethod.pdf")
diff --git a/pointprocesses/lecture/pointprocesses.tex b/pointprocesses/lecture/pointprocesses.tex
index 516e118..7249097 100644
--- a/pointprocesses/lecture/pointprocesses.tex
+++ b/pointprocesses/lecture/pointprocesses.tex
@@ -102,8 +102,8 @@ kann mit den \"ublichen Gr\"o{\ss}en beschrieben werden.
= 1$.
\item Mittleres Intervall: $\mu_{ISI} = \langle T \rangle =
\frac{1}{n}\sum\limits_{i=1}^n T_i$.
-\item Varianz der Intervalle: $\sigma_{ISI}^2 = \langle (T - \langle T
- \rangle)^2 \rangle$\vspace{1ex}
+\item Standardabweichung der Intervalle: $\sigma_{ISI} = \sqrt{\langle (T - \langle T
+ \rangle)^2 \rangle}$\vspace{1ex}
\item Variationskoeffizient (\enterm{coefficient of variation}): $CV_{ISI} =
\frac{\sigma_{ISI}}{\mu_{ISI}}$.
\item Diffusions Koeffizient: $D_{ISI} =
@@ -112,7 +112,7 @@ kann mit den \"ublichen Gr\"o{\ss}en beschrieben werden.
\begin{exercise}{isi_hist.m}{}
Schreibe eine Funktion \code{isi\_hist()}, die einen Vektor mit Interspikeintervallen
- entgegennimmt und daraus ein normalisiertes Histogramm der Interspikeintervalle
+ entgegennimmt und daraus ein normiertes Histogramm der Interspikeintervalle
berechnet.
\end{exercise}
@@ -141,7 +141,7 @@ sichtbar.
Solche Ab\"angigkeiten werden durch die serielle Korrelation der
Intervalle quantifiziert. Das ist der Korrelationskoeffizient
-zwischen aufeinander folgenden Intervallen getrennt durch \enterm{lag} $k$:
+zwischen aufeinander folgenden Intervallen getrennt durch lag $k$:
\[ \rho_k = \frac{\langle (T_{i+k} - \langle T \rangle)(T_i - \langle T \rangle) \rangle}{\langle (T_i - \langle T \rangle)^2\rangle} = \frac{{\rm cov}(T_{i+k}, T_i)}{{\rm var}(T_i)}
= {\rm corr}(T_{i+k}, T_i) \]
\"Ublicherweise wird die Korrelation $\rho_k$ gegen den Lag $k$
@@ -215,10 +215,21 @@ unabh\"angig von den Zeitpunkten fr\"uherer Ereignisse
(\figref{hompoissonfig}). Die Wahrscheinlichkeit zu irgendeiner Zeit
ein Ereigniss in einem kleinen Zeitfenster der Breite $\Delta t$ zu
bekommen ist
-\[ P = \lambda \cdot \Delta t \; . \]
+\begin{equation}
+ \label{hompoissonprob}
+ P = \lambda \cdot \Delta t \; .
+\end{equation}
Beim inhomogenen Poisson Prozess h\"angt die Rate $\lambda$ von der
Zeit ab: $\lambda = \lambda(t)$.
+\begin{exercise}{poissonspikes.m}{}
+ Schreibe eine Funktion \code{poissonspikes()}, die die Spikezeiten
+ eines homogenen Poisson-Prozesses mit gegebener Rate in Hertz f\"ur
+ eine Anzahl von trials gegebener maximaler L\"ange in Sekunden in
+ einem \codeterm{cell-array} zur\"uckgibt. Benutze \eqnref{hompoissonprob}
+ um die Spikezeiten zu bestimmen.
+\end{exercise}
+
\begin{figure}[t]
\includegraphics[width=1\textwidth]{poissonraster100hz}
\titlecaption{\label{hompoissonfig}Rasterplot von Spikes eines homogenen
@@ -234,7 +245,11 @@ Zeit ab: $\lambda = \lambda(t)$.
Der homogene Poissonprozess hat folgende Eigenschaften:
\begin{itemize}
-\item Die Intervalle $T$ sind exponentiell verteilt: $p(T) = \lambda e^{-\lambda T}$ (\figref{hompoissonisihfig}).
+\item Die Intervalle $T$ sind exponentiell verteilt (\figref{hompoissonisihfig}):
+ \begin{equation}
+ \label{poissonintervals}
+ p(T) = \lambda e^{-\lambda T} \; .
+ \end{equation}
\item Das mittlere Intervall ist $\mu_{ISI} = \frac{1}{\lambda}$ .
\item Die Varianz der Intervalle ist $\sigma_{ISI}^2 = \frac{1}{\lambda^2}$ .
\item Der Variationskoeffizient ist also immer $CV_{ISI} = 1$ .
@@ -248,19 +263,20 @@ Der homogene Poissonprozess hat folgende Eigenschaften:
\item Der Fano Faktor ist immer $F=1$ .
\end{itemize}
+\begin{exercise}{hompoissonspikes.m}{}
+ Schreibe eine Funktion \code{hompoissonspikes()}, die die Spikezeiten
+ eines homogenen Poisson-Prozesses mit gegebener Rate in Hertz f\"ur
+ eine Anzahl von trials gegebener maximaler L\"ange in Sekunden in
+ einem \codeterm{cell-array} zur\"uckgibt. Benutze die exponentiell-verteilten
+ Interspikeintervalle \eqnref{poissonintervals}, um die Spikezeiten zu erzeugen.
+\end{exercise}
+
\begin{figure}[t]
\includegraphics[width=0.48\textwidth]{poissoncounthistdist100hz10ms}\hfill
\includegraphics[width=0.48\textwidth]{poissoncounthistdist100hz100ms}
\titlecaption{\label{hompoissoncountfig}Z\"ahlstatistik von Poisson Spikes.}{}
\end{figure}
-\begin{exercise}{poissonspikes.m}{}
- Schreibe eine Funktion \code{poissonspikes()}, die die Spikezeiten
- eines homogenen Poisson-Prozesses mit gegebener Rate in Hertz f\"ur
- eine Anzahl von trials gegebener maximaler L\"ange in Sekunden in
- einem \codeterm{cell-array} zur\"uckgibt.
-\end{exercise}
-
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Zeitabh\"angige Feuerraten}
@@ -295,16 +311,14 @@ Abbildung \ref{psthfig} n\"aher erl\"autert.
\subsection{Instantane Feuerrate}
-
\begin{figure}[tp]
\includegraphics[width=\columnwidth]{isimethod}
- \titlecaption{Bestimmung des zeitabh\"angigen Feuerrate aus dem
- Interspike Intervall.}{\textbf{A)} Skizze eines Rasterplots einer
- einzelnen neuronalen Antwort. Jeder vertikale Strich notiert den
- Zeitpunkt eines Aktionspotentials. Die Pfeile zwischen
- aufeinanderfolgenden Aktionspotentialen und die Zahlen
- illustrieren das Interspike Interval. \textbf{B)} Der Kehrwert
- des Interspike Intervalls ist die Feuerrate.}\label{instrate}
+ \titlecaption{Instantane Feuerrate.}{Skizze eines Spiketrains
+ (oben). Jeder vertikale Strich notiert den Zeitpunkt eines
+ Aktionspotentials. Die Pfeile zwischen aufeinanderfolgenden
+ Aktionspotentialen mit den Zahlen in Millisekunden illustrieren
+ die Interspikeintervalle. Der Kehrwert des Interspikeintervalle
+ ergibt die instantane Feuerrate.}\label{instrate}
\end{figure}
Ein sehr einfacher Weg, die zeitabh\"angige Feuerrate zu bestimmen ist
@@ -321,6 +335,11 @@ Feuerrate k\"onnen f\"ur manche Analysen nachteilig sein. Au{\ss}erdem
wird die Feuerrate nie gleich Null, auch wenn lange keine Aktionspotentiale
generiert wurden.
+\begin{exercise}{instantaneousRate.m}{}
+ Implementiere die Absch\"atzung der Feuerrate auf Basis der
+ instantanen Feuerrate. Plotte die Feuerrate als Funktion der Zeit.
+\end{exercise}
+
\subsection{Peri-Stimulus-Zeit-Histogramm}
W\"ahrend die Instantane Rate den Kehrwert der Zeit von einem bis zum
@@ -343,12 +362,17 @@ oder durch Verfaltung mit einem Kern bestimmt werden. Beiden Methoden
gemeinsam ist die Notwendigkeit der Wahl einer zus\"atzlichen Zeitskala,
die der Beobachtungszeit $W$ in \eqnref{psthrate} entspricht.
-\begin{exercise}{instantaneousRate.m}{}
- Implementiere die Absch\"atzung der Feuerrate auf Basis der
- instantanen Feuerrate. Plotte die Feuerrate als Funktion der Zeit.
-\end{exercise}
-
\subsubsection{Binning-Methode}
+
+\begin{figure}[tp]
+ \includegraphics[width=\columnwidth]{binmethod}
+ \titlecaption{Bestimmung des PSTH mit der Binning Methode.}{Der
+ gleiche Spiketrain wie in \figref{instrate}. Die grauen Linien
+ markieren die Grenzen der Bins und die Zahlen geben die Anzahl der Spikes
+ in jedem Bin an (oben). Die Feuerrate ergibt sich aus dem
+ mit der Binbreite normierten Zeithistogramm (unten).}\label{binpsth}
+\end{figure}
+
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 (\figref{binpsth} A). Um diese
@@ -366,23 +390,23 @@ vorkommen k\"onnen nicht aufgl\"ost werden. Mit der Wahl der Binweite
wird somit eine Annahme \"uber die relevante Zeitskala des Spiketrains
gemacht.
-\begin{figure}[tp]
- \includegraphics[width=\columnwidth]{binmethod}
- \titlecaption{Bestimmung des PSTH mit der Binning Methode.}{\textbf{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}
-
\begin{exercise}{binnedRate.m}{}
Implementiere die Absch\"atzung der Feuerrate mit der ``binning''
Methode. Plotte das PSTH.
\end{exercise}
\subsubsection{Faltungsmethode}
+
+\begin{figure}[tp]
+ \includegraphics[width=\columnwidth]{convmethod}
+ \titlecaption{Bestimmung des PSTH mit der Faltungsmethode.}{Der
+ gleiche Spiketrain wie in \figref{instrate}. Bei der Verfaltung
+ des Spiketrains mit einem Faltungskern wird jeder Spike durch den
+ Faltungskern ersetzt (oben). Bei korrekter Normierung des
+ Kerns ergibt sich die Feuerrate direkt aus der \"Uberlagerung der
+ Kerne.}\label{convrate}
+\end{figure}
+
Bei der Faltungsmethode werden die harten Kanten der Bins der
Binning-Methode vermieden. Der Spiketrain wird mit einem Kern
verfaltet, d.h. jedes Aktionspotential wird durch den Kern ersetzt.
@@ -399,22 +423,12 @@ ersetzt (Abbildung \ref{convrate} A). Wenn der Kern richtig normiert
wurde (Integral gleich Eins), ergibt sich die Feuerrate direkt aus der
\"Uberlagerung der Kerne (Abb. \ref{convrate} B).
-\begin{figure}[tp]
- \includegraphics[width=\columnwidth]{convmethod}
- \titlecaption{Schematische Darstellung der Faltungsmethode.}{\textbf{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}
-
Die Faltungsmethode f\"uhrt, anders als die anderen Methoden, zu einer
-kontinuierlichen Funktion was insbesondere f\"ur spektrale Analysen
-von Vorteil sein kann. Die Wahl der Kernbreite bestimmt, \"ahnlich zur
+stetigen Funktion was insbesondere f\"ur spektrale Analysen von
+Vorteil sein kann. Die Wahl der Kernbreite bestimmt, \"ahnlich zur
Binweite, die zeitliche Aufl\"osung von $r(t)$. Die Breite des Kerns
-macht also auch wieder eine Annahme \"uber die relevante Zeitskala des Spiketrains.
+macht also auch wieder eine Annahme \"uber die relevante Zeitskala des
+Spiketrains.
\begin{exercise}{convolutionRate.m}{}
Verwende die Faltungsmethode um die Feuerrate zu bestimmen. Plotte
@@ -440,12 +454,12 @@ ausgeschnitten wird und diese dann gemittelt werde (\figref{stafig}).
\begin{figure}[t]
\includegraphics[width=\columnwidth]{sta}
\titlecaption{Spike-triggered Average eines P-Typ Elektrorezeptors
- und Stimulusrekonstruktion.}{\textbf{A)} Der STA: der Rezeptor
+ und Stimulusrekonstruktion.}{Der STA (links): der Rezeptor
wurde mit einem ``white-noise'' Stimulus getrieben. Zeitpunkt 0
ist der Zeitpunkt des beobachteten Aktionspotentials. Die Kurve
ergibt sich aus dem gemittelten Stimulus je 50\,ms vor und nach
- einem Aktionspotential. \textbf{B)} Stimulusrekonstruktion mittels
- STA. Die Zellantwort wird mit dem STA gefaltet um eine
+ einem Aktionspotential. Stimulusrekonstruktion mittels
+ STA (rechts). Die Zellantwort wird mit dem STA gefaltet um eine
Rekonstruktion des Stimulus zu erhalten.}\label{stafig}
\end{figure}
diff --git a/scientificcomputing-script.tex b/scientificcomputing-script.tex
index 8eaf472..c1a5a92 100644
--- a/scientificcomputing-script.tex
+++ b/scientificcomputing-script.tex
@@ -67,4 +67,9 @@
%\chapter{Cheat-Sheet}
+\addcontentsline{toc}{chapter}{Fachbegriffe}
+\printindex[term]
+\addcontentsline{toc}{chapter}{Code}
+\printindex[code]
+
\end{document}