diff --git a/bootstrap/lecture/bootstrap-chapter.tex b/bootstrap/lecture/bootstrap-chapter.tex index d185cd4..53099ed 100644 --- a/bootstrap/lecture/bootstrap-chapter.tex +++ b/bootstrap/lecture/bootstrap-chapter.tex @@ -163,9 +163,9 @@ numbers=left, showstringspaces=false, language=Matlab, - commentstyle=\itshape\color{darkgray}, - keywordstyle=\color{blue}, - stringstyle=\color{green}, + commentstyle=\itshape\color{red!60!black}, + keywordstyle=\color{blue!50!black}, + stringstyle=\color{green!50!black}, backgroundcolor=\color{blue!10}, breaklines=true, breakautoindent=true, diff --git a/designpattern/lecture/designpattern-chapter.tex b/designpattern/lecture/designpattern-chapter.tex index 08ff012..e821287 100644 --- a/designpattern/lecture/designpattern-chapter.tex +++ b/designpattern/lecture/designpattern-chapter.tex @@ -163,9 +163,9 @@ numbers=left, showstringspaces=false, language=Matlab, - commentstyle=\itshape\color{darkgray}, - keywordstyle=\color{blue}, - stringstyle=\color{green}, + commentstyle=\itshape\color{red!60!black}, + keywordstyle=\color{blue!50!black}, + stringstyle=\color{green!50!black}, backgroundcolor=\color{blue!10}, breaklines=true, breakautoindent=true, diff --git a/designpattern/lecture/designpattern.tex b/designpattern/lecture/designpattern.tex index c581b9b..93ef588 100644 --- a/designpattern/lecture/designpattern.tex +++ b/designpattern/lecture/designpattern.tex @@ -62,7 +62,7 @@ sigma = 2.3; y = randn(100, 1)*sigma + mu; \end{lstlisting} -Das ist manchmal auch sinnvoll f\"ur \code{zeros} oder \code{ones}: +Das gleiche Prinzip ist manchmal auch sinnvoll f\"ur \code{zeros} oder \code{ones}: \begin{lstlisting} x = -1:0.01:2; % Vektor mit x-Werten plot(x, exp(-x.*x)); @@ -75,11 +75,14 @@ plot(x, ones(size(x))*0.5); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \section{for Schleifen \"uber Vektoren} -Manchmal m\"ochte man doch mit einer for-Schleife \"uber einen Vektor iterieren: +Manchmal m\"ochte man doch mit einer for-Schleife \"uber einen Vektor iterieren. \begin{lstlisting} x = [2:3:20]; % irgendein Vektor -for i=1:length(x) - % Benutze den Wert des Vektors x an der Stelle des Indexes i: +for i=1:length(x) % Mit der for-Schleife "loopen" wir ueber den Vektor + i % das ist der Index der die Elemente des Vektors nacheinander indiziert. + x(i) % das ist der Wert des i-ten Elements des Vektors x. + a = x(i); % die Variable a bekommt den Wert des i-ten Elements des Vektors x zugewiesen. + % Benutze den Wert: do_something( x(i) ); end \end{lstlisting} @@ -89,7 +92,7 @@ sollten wir uns vor der Schleife schon einen Vektor f\"ur die Ergebnisse erstellen: \begin{lstlisting} x = [2:3:20]; % irgendein Vektor -y = zeros(size(x)); % Platz fuer die Ergebnisse +y = zeros(length(x),1); % Platz fuer die Ergebnisse, genauso viele wie Loops der Schleife for i=1:length(x) % Schreibe den Rueckgabewert der Funktion get_something an die i-te % Stelle von y: @@ -125,7 +128,8 @@ for i=1:length(x) % Die Funktion get_something gibt uns einen Vektor zurueck: z = get_something( x(i) ); % dessen Inhalt h\"angen wir an unseren Ergebnissvektor an: - y = [y z(:)]; + y = [y; z(:)]; + % z(:) stellt sicher, das wir auf jeden Fall einen Spaltenvektoren aneinanderreihen. end % jetzt koennen wir dem Ergebnisvektor weiter bearbeiten: mean(y) diff --git a/likelihood/lecture/likelihood-chapter.tex b/likelihood/lecture/likelihood-chapter.tex index 732acbe..2a82f79 100644 --- a/likelihood/lecture/likelihood-chapter.tex +++ b/likelihood/lecture/likelihood-chapter.tex @@ -163,9 +163,9 @@ numbers=left, showstringspaces=false, language=Matlab, - commentstyle=\itshape\color{darkgray}, - keywordstyle=\color{blue}, - stringstyle=\color{green}, + commentstyle=\itshape\color{red!60!black}, + keywordstyle=\color{blue!50!black}, + stringstyle=\color{green!50!black}, backgroundcolor=\color{blue!10}, breaklines=true, breakautoindent=true, diff --git a/pointprocesses/code/counthist.m b/pointprocesses/code/counthist.m index 9abb225..15d431e 100644 --- a/pointprocesses/code/counthist.m +++ b/pointprocesses/code/counthist.m @@ -1,24 +1,36 @@ function [counts, bins] = counthist(spikes, w) -% computes count histogram and compare them with Poisson distribution -% spikes: a cell array of vectors of spike times -% w: observation window duration for computing the counts +% computes count histogram and compare them with Poisson distribution +% +% [counts, bins] = counthist(spikes, w) +% spikes: a cell array of vectors of spike times in seconds +% w: observation window duration in seconds for computing the counts +% 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) - for tk = 0:w:tmax-w - nn = sum( ( spikes{k} >= tk ) & ( spikes{k} < tk+w ) ); - %nn = length( find( ( spikes{k} >= tk ) & ( spikes{k} < tk+w ) ) ); - n = [ n nn ]; - end - rate = (length(spikes{k})-1)/(spikes{k}(end) - spikes{k}(1)); + 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 ); if nargout == 0 bar( bins, counts ); @@ -26,12 +38,11 @@ function [counts, bins] = counthist(spikes, w) % Poisson distribution: rate = mean( r ); x = 0:1:maxn+10; - l = rate*w; - y = l.^x.*exp(-l)./factorial(x); + 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/isihist.m b/pointprocesses/code/isihist.m index afd69d4..17974cc 100644 --- a/pointprocesses/code/isihist.m +++ b/pointprocesses/code/isihist.m @@ -1,7 +1,9 @@ function isihist( isis, binwidth ) -% plot histogram of isis -% isis: vector of interspike intervals -% binwidth: optional width to be used for the isi bins +% plot histogram of interspike intervals +% +% isihist(isis, binwidth) +% isis: vector of interspike intervals in seconds +% binwidth: optional width in seconds to be used for the isi bins if nargin < 2 nperbin = 200; % average number of data points per bin diff --git a/pointprocesses/code/isis.m b/pointprocesses/code/isis.m index b0b95ff..e693ea8 100644 --- a/pointprocesses/code/isis.m +++ b/pointprocesses/code/isis.m @@ -1,15 +1,16 @@ function isivec = isis( spikes ) % returns a single list of isis computed from all trials in spikes -% spikes: a cell array of vectors of spike times +% +% isivec = isis( spikes ) +% spikes: a cell array of vectors of spike times in seconds +% isivec: a column vector with all the interspike intervalls isivec = []; for k = 1:length(spikes) difftimes = diff( spikes{k} ); - if ( size( difftimes, 1 ) == 1 ) - isivec = [ isivec difftimes ]; - elseif ( size( difftimes, 2 ) == 1 ) - isivec = [ isivec difftimes' ]; - end + % difftimes(:) ensures a column vector + % regardless of the type of vector spikes{k} + isivec = [ isivec; difftimes(:) ]; end end diff --git a/pointprocesses/code/isiserialcorr.m b/pointprocesses/code/isiserialcorr.m index d5b44c1..c93fb5e 100644 --- a/pointprocesses/code/isiserialcorr.m +++ b/pointprocesses/code/isiserialcorr.m @@ -1,15 +1,19 @@ function isicorr = isiserialcorr( isis, maxlag ) -% serial correlation of isis -% isis: vector of interspike intervals -% maxlag: the maximum lag +% serial correlation of interspike intervals +% +% isicorr = isiserialcorr( isis, maxlag ) +% isis: vector of interspike intervals in seconds +% maxlag: the maximum lag in seconds +% isicorr: vector with the serial correlations for lag 0 to maxlag lags = 0:maxlag; isicorr = zeros( size( lags ) ); for k = 1:length(lags) lag = lags(k); - if length( isis ) > lag+10 - cc = corrcoef( [ isis(1:end-lag)', isis(1+lag:end)' ] ); - isicorr(k) = cc( 1, 2 ); + if length( isis ) > lag+10 % ensure "enough" data + % DANGER: 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) ); end end diff --git a/pointprocesses/code/plotspikestats.m b/pointprocesses/code/plotspikestats.m new file mode 100644 index 0000000..a664eca --- /dev/null +++ b/pointprocesses/code/plotspikestats.m @@ -0,0 +1,100 @@ +%% load data: +clear all +% alternative 1: +% pro: no structs. contra: global unknown variables +load poisson.mat +whos +poissonspikes = spikes; +load pifou.mat; +pifouspikes = spikes; +load lifadapt.mat; +lifadaptspikes = spikes; +clear spikes; +% alternative 2: +% pro: clean code. contra: structs that we do not really know yet +clear all +x = load( 'poisson.mat' ); +poissonspikes = x.spikes; +x = load( 'pifou.mat' ); +pifouspikes = x.spikes; +x = load( 'lifadapt.mat' ); +lifadaptspikes = x.spikes; + +%% spike raster plots: +tmax = 1.0; +subplot(1, 3, 1); +spikeraster(poissonspikes, tmax); +title('Poisson'); + +subplot(1, 3, 2); +spikeraster(pifouspikes, tmax); +title('PIF OU'); + +subplot(1, 3, 3); +spikeraster(lifadaptspikes, tmax); +title('LIF adapt'); + +%% isi histograms: +maxisi = 300.0; +binwidth = 0.002; +subplot(1, 3, 1); +poissonisis = isis(poissonspikes); +isihist(poissonisis, binwidth); +xlim([0, maxisi]) +title('Poisson'); + +subplot(1, 3, 2); +pifouisis = isis(pifouspikes); +isihist(pifouisis, binwidth); +xlim([0, maxisi]) +title('PIF OU'); + +subplot(1, 3, 3); +lifadaptisis = isis(lifadaptspikes); +isihist(lifadaptisis, binwidth); +xlim([0, maxisi]) +title('LIF adapt'); + +%% serial correlations: +maxlag = 10; +rrange = [-0.5, 1.05]; +subplot(1, 3, 1); +isiserialcorr(poissonisis, maxlag); +ylim(rrange) +title('Poisson'); + +subplot(1, 3, 2); +isiserialcorr(pifouisis, maxlag); +ylim(rrange) +title('PIF OU'); + +subplot(1, 3, 3); +isiserialcorr(lifadaptisis, maxlag); +ylim(rrange) +title('LIF adapt'); + +%% spike counts: +w = 0.1; +cmax = 8; +pmax = 0.5; +subplot(1, 3, 1); +counthist(poissonspikes, w); +xlim([0 cmax]) +set(gca, 'XTick', 0:2:cmax) +ylim([0 pmax]) +title('Poisson'); + +subplot(1, 3, 2); +counthist(pifouspikes, w); +xlim([0 cmax]) +set(gca, 'XTick', 0:2:cmax) +ylim([0 pmax]) +title('PIF OU'); + +subplot(1, 3, 3); +counthist(lifadaptspikes, w); +xlim([0 cmax]) +set(gca, 'XTick', 0:2:cmax) +ylim([0 pmax]) +title('LIF adapt'); +savefigpdf(gcf, 'counthist.pdf', 20, 7); diff --git a/pointprocesses/code/poissonspikes.m b/pointprocesses/code/poissonspikes.m index 6a75f40..37f2234 100644 --- a/pointprocesses/code/poissonspikes.m +++ b/pointprocesses/code/poissonspikes.m @@ -1,9 +1,11 @@ function spikes = poissonspikes( trials, rate, tmax ) % Generate spike times of a homogeneous poisson process -% trials: number of trials that should be generated -% rate: the rate of the Poisson process in Hertz -% tmax: the duration of each trial in seconds -% returns a cell array of vectors of spike times +% +% spikes = poissonspikes( trials, rate, tmax ) +% trials: number of trials that should be generated +% rate: the rate of the Poisson process in Hertz +% tmax: the duration of each trial in seconds +% spikes: a cell array of vectors of spike times in seconds dt = 3.33e-5; p = rate*dt; % probability of event per bin of width dt diff --git a/pointprocesses/code/spikeraster.m b/pointprocesses/code/spikeraster.m index 7d4c12f..cde46dc 100644 --- a/pointprocesses/code/spikeraster.m +++ b/pointprocesses/code/spikeraster.m @@ -1,7 +1,9 @@ function spikeraster(spikes, tmax) % Display a spike raster of the spike times given in spikes. -% spikes: a cell array of vectors of spike times -% tmax: plot spike raster upto tmax seconds +% +% spikeraster(spikes, tmax) +% spikes: a cell array of vectors of spike times in seconds +% tmax: plot spike raster upto tmax seconds ntrials = length(spikes); for k = 1:ntrials @@ -16,8 +18,10 @@ for k = 1:ntrials end if tmax < 1.5 xlabel( 'Time [ms]' ); + xlim([0.0 1000.0*tmax]); else xlabel( 'Time [s]' ); + xlim([0.0 tmax]); end ylabel( 'Trials'); ylim( [ 0.3 ntrials+0.7 ] ) diff --git a/pointprocesses/lecture/pointprocesses-chapter.tex b/pointprocesses/lecture/pointprocesses-chapter.tex index 44af27e..1e68316 100644 --- a/pointprocesses/lecture/pointprocesses-chapter.tex +++ b/pointprocesses/lecture/pointprocesses-chapter.tex @@ -163,9 +163,9 @@ numbers=left, showstringspaces=false, language=Matlab, - commentstyle=\itshape\color{darkgray}, - keywordstyle=\color{blue}, - stringstyle=\color{green}, + commentstyle=\itshape\color{red!60!black}, + keywordstyle=\color{blue!50!black}, + stringstyle=\color{green!50!black}, backgroundcolor=\color{blue!10}, breaklines=true, breakautoindent=true, diff --git a/scientificcomputing-script.tex b/scientificcomputing-script.tex index 65e8b2b..e4049a8 100644 --- a/scientificcomputing-script.tex +++ b/scientificcomputing-script.tex @@ -163,9 +163,9 @@ numbers=left, showstringspaces=false, language=Matlab, - commentstyle=\itshape\color{darkgray}, - keywordstyle=\color{blue}, - stringstyle=\color{green}, + commentstyle=\itshape\color{red!60!black}, + keywordstyle=\color{blue!50!black}, + stringstyle=\color{green!50!black}, backgroundcolor=\color{blue!10}, breaklines=true, breakautoindent=true, diff --git a/statistics/lecture/diehistograms.py b/statistics/lecture/diehistograms.py index a8832e0..c80c68a 100644 --- a/statistics/lecture/diehistograms.py +++ b/statistics/lecture/diehistograms.py @@ -2,8 +2,9 @@ import numpy as np import matplotlib.pyplot as plt # roll the die: -x1 = np.random.random_integers( 1, 6, 100 ) -x2 = np.random.random_integers( 1, 6, 500 ) +rng = np.random.RandomState(57281) +x1 = rng.random_integers( 1, 6, 100 ) +x2 = rng.random_integers( 1, 6, 500 ) bins = np.arange(0.5, 7, 1.0) plt.xkcd() @@ -14,7 +15,10 @@ ax.spines['right'].set_visible(False) ax.spines['top'].set_visible(False) ax.yaxis.set_ticks_position('left') ax.xaxis.set_ticks_position('bottom') +ax.set_xlim(0, 7) +ax.set_xticks( range(1, 7) ) ax.set_xlabel( 'x' ) +ax.set_ylim(0, 98) ax.set_ylabel( 'Frequency' ) ax.hist([x2, x1], bins, color=['#FFCC00', '#FFFF66' ]) @@ -23,9 +27,13 @@ ax.spines['right'].set_visible(False) ax.spines['top'].set_visible(False) ax.yaxis.set_ticks_position('left') ax.xaxis.set_ticks_position('bottom') +ax.set_xlim(0, 7) +ax.set_xticks( range(1, 7) ) ax.set_xlabel( 'x' ) +ax.set_ylim(0, 0.23) ax.set_ylabel( 'Probability' ) -ax.hist([x2, x1], bins, normed=True, color=['#FFCC00', '#FFFF66' ]) +ax.plot([0.2, 6.8], [1.0/6.0, 1.0/6.0], '-r', lw=2, zorder=1) +ax.hist([x2, x1], bins, normed=True, color=['#FFCC00', '#FFFF66' ], zorder=10) plt.tight_layout() fig.savefig( 'diehistograms.pdf' ) #plt.show() diff --git a/statistics/lecture/pdfhistogram.py b/statistics/lecture/pdfhistogram.py index a61460e..1faa52d 100644 --- a/statistics/lecture/pdfhistogram.py +++ b/statistics/lecture/pdfhistogram.py @@ -2,9 +2,10 @@ import numpy as np import matplotlib.pyplot as plt # normal distribution: +rng = np.random.RandomState(6281) x = np.arange( -4.0, 4.0, 0.01 ) g = np.exp(-0.5*x*x)/np.sqrt(2.0*np.pi) -r = np.random.randn( 100 ) +r = rng.randn( 100 ) plt.xkcd() @@ -15,9 +16,10 @@ ax.spines['top'].set_visible(False) ax.yaxis.set_ticks_position('left') ax.xaxis.set_ticks_position('bottom') ax.set_xlabel( 'x' ) +ax.set_xlim(-3.2, 3.2) +ax.set_xticks( np.arange( -3.0, 3.1, 1.0 ) ) ax.set_ylabel( 'Frequency' ) -#ax.set_ylim( 0.0, 0.46 ) -#ax.set_yticks( np.arange( 0.0, 0.45, 0.1 ) ) +ax.set_yticks( np.arange( 0.0, 41.0, 10.0 ) ) ax.hist(r, 5, color='#CC0000') ax.hist(r, 20, color='#FFCC00') @@ -27,11 +29,14 @@ ax.spines['top'].set_visible(False) ax.yaxis.set_ticks_position('left') ax.xaxis.set_ticks_position('bottom') ax.set_xlabel( 'x' ) +ax.set_xlim(-3.2, 3.2) +ax.set_xticks( np.arange( -3.0, 3.1, 1.0 ) ) ax.set_ylabel( 'Probability density p(x)' ) -#ax.set_ylim( 0.0, 0.46 ) -#ax.set_yticks( np.arange( 0.0, 0.45, 0.1 ) ) -ax.hist(r, 5, normed=True, color='#CC0000') -ax.hist(r, 20, normed=True, color='#FFCC00') +ax.set_ylim(0.0, 0.44) +ax.set_yticks( np.arange( 0.0, 0.41, 0.1 ) ) +ax.plot(x, g, '-b', lw=2, zorder=-1) +ax.hist(r, 5, normed=True, color='#CC0000', zorder=-10) +ax.hist(r, 20, normed=True, color='#FFCC00', zorder=-5) plt.tight_layout() fig.savefig( 'pdfhistogram.pdf' ) diff --git a/statistics/lecture/statistics-chapter.tex b/statistics/lecture/statistics-chapter.tex index d43ca3e..1d13265 100644 --- a/statistics/lecture/statistics-chapter.tex +++ b/statistics/lecture/statistics-chapter.tex @@ -163,9 +163,9 @@ numbers=left, showstringspaces=false, language=Matlab, - commentstyle=\itshape\color{darkgray}, - keywordstyle=\color{blue}, - stringstyle=\color{green}, + commentstyle=\itshape\color{red!60!black}, + keywordstyle=\color{blue!50!black}, + stringstyle=\color{green!50!black}, backgroundcolor=\color{blue!10}, breaklines=true, breakautoindent=true, diff --git a/statistics/lecture/statistics.tex b/statistics/lecture/statistics.tex index 85ce58d..2ca6481 100644 --- a/statistics/lecture/statistics.tex +++ b/statistics/lecture/statistics.tex @@ -111,7 +111,8 @@ Wahrscheinlichkeitsverteilung der Messwerte ab. to their sum.}{Histogramme des Ergebnisses von 100 oder 500 mal W\"urfeln. Links: das absolute Histogramm z\"ahlt die Anzahl des Auftretens jeder Augenzahl. Rechts: Normiert auf die Summe des - Histogramms werden die beiden Messungen vergleichbar.}} + Histogramms werden die beiden Messungen untereinander als auch + mit der theoretischen Verteilung $P=1/6$ vergleichbar.}} \end{figure} Bei ganzzahligen Messdaten (z.B. die Augenzahl eines W\"urfels) @@ -142,7 +143,9 @@ Meistens haben wir es jedoch mit reellen Messgr\"o{\ss}en zu tun. normalverteilten Messwerten. Links: Die H\"ohe des absoluten Histogramms h\"angt von der Klassenbreite ab. Rechts: Bei auf das Integral normierten Histogrammen werden auch - unterschiedliche Klassenbreiten vergleichbar.}} + unterschiedliche Klassenbreiten untereinander vergleichbar und + auch mit der theoretischen Wahrschinlichkeitsdichtefunktion + (blau).}} \end{figure} Histogramme von reellen Messwerten m\"ussen auf das Integral 1