Merge branch 'master' of raven.am28.uni-tuebingen.de:scientificComputing

bootstrap/lecture/bootstrap-chapter.tex

@@ -163,9 +163,9 @@
   numbers=left,
   showstringspaces=false,
   language=Matlab,
-  commentstyle=\itshape\color{darkgray},
+  commentstyle=\itshape\color{red!60!black},
-  keywordstyle=\color{blue},
+  keywordstyle=\color{blue!50!black},
-  stringstyle=\color{green},
+  stringstyle=\color{green!50!black},
   backgroundcolor=\color{blue!10},
   breaklines=true,
   breakautoindent=true,

designpattern/lecture/designpattern-chapter.tex

@@ -163,9 +163,9 @@
   numbers=left,
   showstringspaces=false,
   language=Matlab,
-  commentstyle=\itshape\color{darkgray},
+  commentstyle=\itshape\color{red!60!black},
-  keywordstyle=\color{blue},
+  keywordstyle=\color{blue!50!black},
-  stringstyle=\color{green},
+  stringstyle=\color{green!50!black},
   backgroundcolor=\color{blue!10},
   breaklines=true,
   breakautoindent=true,

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)
+for i=1:length(x)     % Mit der for-Schleife "loopen" wir ueber den Vektor
-  % Benutze den Wert des Vektors x an der Stelle des Indexes i:
+  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)

likelihood/lecture/likelihood-chapter.tex

@@ -163,9 +163,9 @@
   numbers=left,
   showstringspaces=false,
   language=Matlab,
-  commentstyle=\itshape\color{darkgray},
+  commentstyle=\itshape\color{red!60!black},
-  keywordstyle=\color{blue},
+  keywordstyle=\color{blue!50!black},
-  stringstyle=\color{green},
+  stringstyle=\color{green!50!black},
   backgroundcolor=\color{blue!10},
   breaklines=true,
   breakautoindent=true,

pointprocesses/code/counthist.m

@@ -1,24 +1,36 @@
 function [counts, bins] = counthist(spikes, w)
-% computes count histogram and compare them  with Poisson distribution
+% 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
+% [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
+        times = spikes{k};
-            nn = sum( ( spikes{k} >= tk ) & ( spikes{k} < tk+w ) );
+%       alternative 1: count the number of spikes in each window:
-            %nn = length( find( ( spikes{k} >= tk ) & ( spikes{k} < tk+w ) ) );
+%         for tk = 0:w:tmax-w
-            n = [ n nn ];
+%             nn = sum( ( times >= tk ) & ( times < tk+w ) );
-        end
+%             %nn = length( find( ( times >= tk ) & ( times < tk+w ) ) );
-        rate = (length(spikes{k})-1)/(spikes{k}(end) - spikes{k}(1));
+%             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;
+        a = rate*w;
-        y = l.^x.*exp(-l)./factorial(x);
+        y = a.^x.*exp(-a)./factorial(x);
         plot( x, y, 'r', 'LineWidth', 3 );
         hold off;
         xlabel( 'counts k' );
         ylabel( 'P(k)' );
     end
 end
-

pointprocesses/code/isihist.m

@@ -1,7 +1,9 @@
 function isihist( isis, binwidth )
-% plot histogram of isis
+% plot histogram of interspike intervals
-% isis: vector of interspike intervals
+%
-% binwidth: optional width to be used for the isi bins
+% 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

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 )
+        % difftimes(:) ensures a column vector
-            isivec = [ isivec difftimes ];
+        % regardless of the type of vector spikes{k}
-        elseif ( size( difftimes, 2 ) == 1 )
+        isivec = [ isivec; difftimes(:) ];
-            isivec = [ isivec difftimes' ];
-        end
     end
 end
 

pointprocesses/code/isiserialcorr.m

@@ -1,15 +1,19 @@
 function isicorr = isiserialcorr( isis, maxlag )
-% serial correlation of isis
+% serial correlation of interspike intervals
-% isis: vector of interspike intervals
+%
-% maxlag: the maximum lag
+% 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
+        if length( isis ) > lag+10        % ensure "enough" data
-            cc = corrcoef( [ isis(1:end-lag)', isis(1+lag:end)' ] );
+	    % DANGER: the arguments to corr must be column vectors!
-            isicorr(k) = cc( 1, 2 );
+            % We insure this in the isis() function that generats the isis.
+            isicorr(k) = corr( isis(1:end-lag), isis(lag+1:end) );
         end
     end
     

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);

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
+% spikes = poissonspikes( trials, rate, tmax )
-% tmax: the duration of each trial in seconds
+%   trials: number of trials that should be generated
-% returns a cell array of vectors of spike times
+%   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

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 ] )

pointprocesses/lecture/pointprocesses-chapter.tex

@@ -163,9 +163,9 @@
   numbers=left,
   showstringspaces=false,
   language=Matlab,
-  commentstyle=\itshape\color{darkgray},
+  commentstyle=\itshape\color{red!60!black},
-  keywordstyle=\color{blue},
+  keywordstyle=\color{blue!50!black},
-  stringstyle=\color{green},
+  stringstyle=\color{green!50!black},
   backgroundcolor=\color{blue!10},
   breaklines=true,
   breakautoindent=true,

scientificcomputing-script.tex

@@ -163,9 +163,9 @@
   numbers=left,
   showstringspaces=false,
   language=Matlab,
-  commentstyle=\itshape\color{darkgray},
+  commentstyle=\itshape\color{red!60!black},
-  keywordstyle=\color{blue},
+  keywordstyle=\color{blue!50!black},
-  stringstyle=\color{green},
+  stringstyle=\color{green!50!black},
   backgroundcolor=\color{blue!10},
   breaklines=true,
   breakautoindent=true,

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 )
+rng = np.random.RandomState(57281)
-x2 = np.random.random_integers( 1, 6, 500 )
+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()

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, 41.0, 10.0 ) )
-#ax.set_yticks( np.arange( 0.0, 0.45, 0.1 ) )
 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_ylim(0.0, 0.44)
-#ax.set_yticks( np.arange( 0.0, 0.45, 0.1 ) )
+ax.set_yticks( np.arange( 0.0, 0.41, 0.1 ) )
-ax.hist(r, 5, normed=True, color='#CC0000')
+ax.plot(x, g, '-b', lw=2, zorder=-1)
-ax.hist(r, 20, normed=True, color='#FFCC00')
+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' )

statistics/lecture/statistics-chapter.tex

@@ -163,9 +163,9 @@
   numbers=left,
   showstringspaces=false,
   language=Matlab,
-  commentstyle=\itshape\color{darkgray},
+  commentstyle=\itshape\color{red!60!black},
-  keywordstyle=\color{blue},
+  keywordstyle=\color{blue!50!black},
-  stringstyle=\color{green},
+  stringstyle=\color{green!50!black},
   backgroundcolor=\color{blue!10},
   breaklines=true,
   breakautoindent=true,

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