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}