From d717b436303e22bae07a17bc349fe039f4248e42 Mon Sep 17 00:00:00 2001
From: Jan Benda <jan.benda@uni-tuebingen.de>
Date: Mon, 28 Nov 2016 17:33:27 +0100
Subject: [PATCH] new statistics exercises
---
pointprocesses/exercises/pointprocesses01.tex | 4 +-
statistics/exercises/centrallimit.m | 48 ++---
statistics/exercises/exercises02.tex | 186 +++++++++++-------
statistics/exercises/normprobs.m | 49 ++++-
4 files changed, 183 insertions(+), 104 deletions(-)
diff --git a/pointprocesses/exercises/pointprocesses01.tex b/pointprocesses/exercises/pointprocesses01.tex
index e537be7..ea8674d 100644
--- a/pointprocesses/exercises/pointprocesses01.tex
+++ b/pointprocesses/exercises/pointprocesses01.tex
@@ -98,8 +98,8 @@ jan.benda@uni-tuebingen.de}
Mit den folgenden Aufgaben wollen wir die Statistik der Spiketrains
der drei Neurone miteinander vergleichen.
\begin{parts}
- \part Lade die Spiketrains aus den drei Dateien. Achte darauf, dass sie verschiedene
- Variablen\-namen bekommen.
+ \part Lade die Spiketrains aus den drei Dateien. Achte darauf,
+ dass sie verschiedene Variablen\-namen bekommen.
\begin{solution}
\begin{lstlisting}
clear all
diff --git a/statistics/exercises/centrallimit.m b/statistics/exercises/centrallimit.m
index 0794335..cfc5f0d 100644
--- a/statistics/exercises/centrallimit.m
+++ b/statistics/exercises/centrallimit.m
@@ -3,36 +3,36 @@ n = 10000;
m = 10; % number of loops
%% (b) a single data set of random numbers:
-x = rand( n, 1 );
+x = rand(n, 1);
%% (c) plot probability density:
-%histogram( x, 'Normalization', 'pdf' );
-[h,b] = hist( x, 20 );
+%histogram(x, 'Normalization', 'pdf');
+[h,b] = hist(x, 20);
h = h/sum(h)/(b(2)-b(1)); % normalization
bar(b, h)
title('A uniform distribution')
xlabel('x')
ylabel('probability density')
-pause( 2.0 )
+pause(2.0)
%% (d) sum of two random numbers:
-y = rand( n, 1 );
+y = rand(n, 1);
x = x + y;
-%histogram( x, 'Normalization', 'pdf' );
-[h,b] = hist( x, 20 );
+%histogram(x, 'Normalization', 'pdf');
+[h,b] = hist(x, 20);
h = h/sum(h)/(b(2)-b(1)); % normalization
bar(b, h)
title('Sum of two uniform distributions')
xlabel('x')
ylabel('probability density')
-pause( 2.0 )
+pause(2.0)
%% (f) sum up more distributions:
-x = zeros( n, 1 );
-means = zeros( m, 1 );
-stds = zeros( m, 1 );
+x = zeros(n, 1);
+means = zeros(m, 1);
+stds = zeros(m, 1);
for i=1:m
- y = rand( n, 1 ); % new uniform distributed numbers
+ y = rand(n, 1); % new uniform distributed numbers
x = x + y; % add them to the sum
mu = mean(x); % compute mean
sd = std(x); % compute standard deviation
@@ -42,34 +42,34 @@ for i=1:m
xx = -1:0.01:i+1; % x-axis values for plot of pdf
p = exp(-0.5*(xx-mu).^2/sd^2)/sqrt(2*pi*sd^2); % pdf
plot(xx, p, 'r', 'linewidth', 3 )
- ns = sprintf( 'N=%d', i );
- text( 0.1, 0.9, ns, 'units', 'normalized' )
+ ns = sprintf('N=%d', i);
+ text(0.1, 0.9, ns, 'units', 'normalized')
hold on
- %histogram( x, 20, 'Normalization', 'pdf' );
- [h,b] = hist( x, 20 );
+ %histogram(x, 20, 'Normalization', 'pdf');
+ [h,b] = hist(x, 20);
h = h/sum(h)/(b(2)-b(1)); % normalization
bar(b, h)
hold off
xlim([-0.5, i+0.5])
- xlabel( 'x' )
- ylabel( 'summed pdf' )
- savefigpdf( gcf, sprintf('centrallimit-hist%02d.pdf', i), 6, 5 );
+ xlabel('x')
+ ylabel('summed pdf')
+ savefigpdf(gcf, sprintf('centrallimit-hist%02d.pdf', i), 6, 5);
if i < 6
- pause( 3.0 )
+ pause(3.0)
end
end
%% (h) mean and standard deviation in dependence on number of summed up distributions:
nx = 1:m;
-plot( nx, means, 'b', 'linewidth', 4 );
+plot(nx, means, 'b', 'linewidth', 4);
hold on
-plot( nx, stds, 'r', 'linewidth', 4 );
+plot(nx, stds, 'r', 'linewidth', 4);
xx = 0:0.01:m;
sdu = 1.0/sqrt(12); % standarad deviation of the uniform distribution
plot( xx, sqrt(xx)*sdu, 'k' )
-legend( 'mean', 'std', 'theory' )
+legend('mean', 'std', 'theory')
xlabel('N')
hold off
-savefigpdf( gcf, 'centrallimit-samples.pdf', 6, 5 );
+savefigpdf(gcf, 'centrallimit-samples.pdf', 6, 5);
diff --git a/statistics/exercises/exercises02.tex b/statistics/exercises/exercises02.tex
index 1aa4b20..82e2528 100644
--- a/statistics/exercises/exercises02.tex
+++ b/statistics/exercises/exercises02.tex
@@ -15,7 +15,7 @@
\else
\newcommand{\stitle}{}
\fi
-\header{{\bfseries\large \"Ubung 6X\stitle}}{{\bfseries\large Statistik}}{{\bfseries\large 22. November, 2016}}
+\header{{\bfseries\large \"Ubung 7\stitle}}{{\bfseries\large Statistik}}{{\bfseries\large 29. November, 2016}}
\firstpagefooter{Prof. Dr. Jan Benda}{Phone: 29 74573}{Email:
jan.benda@uni-tuebingen.de}
\runningfooter{}{\thepage}{}
@@ -78,6 +78,8 @@ jan.benda@uni-tuebingen.de}
\newcommand{\extra}{--- Zusatzaufgabe ---\ \mbox{}}
\newcommand{\code}[1]{\texttt{#1}}
+\graphicspath{{../../pointprocesses/exercises/}}
+
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{document}
@@ -96,13 +98,16 @@ Mittelwert enthalten ist.
\part Erzeuge einen Datensatz $X = (x_1, x_2, ... x_n)$ aus
$n=10000$ normalverteilten Zufallszahlen mit Mittelwert $\mu=0$ und
Standardabweichung $\sigma=1$ (\code{randn() Funktion}).
- \part Bestimme und plotte die Wahrscheinlichkeitsdichte dieser Zufallszahlen (normiertes Histogramm).
- \part Plotte zum Vergleich in den gleichen Plot die Normalverteilung
+ \part Bestimme und plotte die Wahrscheinlichkeitsdichte dieser
+ Zufallszahlen (normiertes Histogramm) und plotte zum Vergleich in
+ den gleichen Plot die Normalverteilung
\[ p_g(x) = \frac{1}{\sqrt{2\pi\sigma^2}}e^{-\frac{1}{2}\left(\frac{x-\mu}{\sigma}\right)^2} \; . \]
- \part \label{onesigma} Wieviele dieser Daten sind maximal eine Standardabweichung vom Mittelwert entfernt?\\
+
+ \part \label{onesigma} Wieviele dieser Daten $X$ sind maximal eine Standardabweichung vom Mittelwert entfernt?\\
D.h. wieviele Datenwerte $x_i$ haben den Wert $-\sigma < x_i < +\sigma$?\\
- Wie gro{\ss} ist also die Wahrscheinlichkeit $P_{\pm\sigma}$ einen
+ Wie gro{\ss} ist dann also die Wahrscheinlichkeit $P_{\pm\sigma}$ einen
Wert in diesem Interval zu erhalten?
+
\part \label{probintegral} Berechne numerisch diese
Wahrscheinlichkeit aus dem entsprechenden Integral
\[ P_{\pm\sigma}=\int_{x=\mu-\sigma}^{x=\mu+\sigma} p_g(x) \, dx \]
@@ -110,16 +115,24 @@ Mittelwert enthalten ist.
\"Uberpr\"ufe zuerst, ob tats\"achlich
\[ \int_{-\infty}^{+\infty} p_g(x) \, dx = 1 \; . \]
Warum muss das so sein?
- \part Welcher Anteil der Daten ist in den Intervallen $\pm2\sigma$ sowie $\pm3\sigma$
- enthalten?
- \part \label{givenfraction} Finde heraus in welchem Interval symmetrisch um den Mittelwert
- 50\,\%, 90\,\%, 95\,\% bzw. 99\,\% der Daten enhalten sind.
- \part \extra Modifiziere den Code der Teilaufgaben \pref{onesigma}
- -- \pref{givenfraction} so, dass er f\"ur Datens\"atze mit
- beliebigen Mittelwerten und Standardabweichungen funktioniert.\\
- Teste den Code mit entsprechenden Zufallszahlen.\\
- Wie bekommt man mit \code{randn()} Zufallszahlen mit beliebiger
- Standardabweichung und Mittelwerten?
+
+ \part Welcher Anteil der Daten ist in den Intervallen $\pm 2\sigma$
+ sowie $\pm 3\sigma$ enthalten?
+
+ Vergleiche die Ergebnisse jeweils mit dem entsprechenden Integral
+ \"uber die Wahrscheinlichkeitsdichte.
+
+ \part \label{givenfraction} Finde durch numerische Integration der
+ Wahrscheinlichkeitsdichte heraus, in welchem Interval symmetrisch um
+ den Mittelwert 50\,\%, 90\,\%, 95\,\% bzw. 99\,\% der Daten enhalten
+ sind.
+
+ % \part \extra Modifiziere den Code der Teilaufgaben \pref{onesigma}
+ % -- \pref{givenfraction} so, dass er f\"ur Datens\"atze mit
+ % beliebigen Mittelwerten und Standardabweichungen funktioniert.\\
+ % Teste den Code mit entsprechenden Zufallszahlen.\\
+ % Wie bekommt man mit \code{randn()} Zufallszahlen mit beliebiger
+ % Standardabweichung und Mittelwerten?
\end{parts}
\begin{solution}
\lstinputlisting{normprobs.m}
@@ -134,30 +147,40 @@ distributed) Zufallsvariablen gegen die Normalverteilung konvergiert.
Den Zentralen Grenzwertsatz wollen wir uns im Folgenden veranschaulichen.
\begin{parts}
- \part Versuche dir klar zu machen, was der Zentrale Grenzwertsatz
- bedeutet, und wie du vorgehen k\"onntest ein Programm zu
- schreiben, das den Grenzwertsatz illustriert.
+ \part Bevor du die weiteren Teilaufgaben liest, versuche dir klar zu
+ machen, was der Zentrale Grenzwertsatz bedeutet, und wie du vorgehen
+ k\"onntest ein Programm zu schreiben, das den Grenzwertsatz
+ illustriert.
+
\part Erzeuge 10000 zwischen 0 und 1 gleichverteilte Zufallszahlen
(Funktion \code{rand}).
+
\part Plotte deren Wahrscheinlichkeitsdichte (normiertes Histogram).
+
\part Erzeuge weitere 10000 gleichverteilte Zufallszahlen und
addiere diese zu den bereits vorhandenen auf.
+
\part Plotte die Wahrscheinlichkeitsdichte der aufsummierten
Zufallszahlen.
+
\part Wiederhole Schritt (d) und (e) viele Male.
+
\part Vergleiche in einer Grafik die Wahrscheinlichkeitsdichte der
aufsummierten Zufallszahlen mit der Gaussfunktion
\[ p_g(x) =
\frac{1}{\sqrt{2\pi\sigma^2}}e^{-\frac{1}{2}\left(\frac{x-\mu}{\sigma}\right)^2}\]
mit dem Mittelwert $\mu$ und der Standardabweichung $\sigma$ der
aufsummierten Zufallszahlen.
+
\part Wie \"andert sich der Mittelwert und die
Standardabweichung/Varianz
- der aufsummierten Zufallszahlen?\\
+ der aufsummierten Zufallszahlen?
+
Wie h\"angen diese mit den Werten der urspr\"unglichen Verteilung
zusammen?
- \part \extra \"Uberpr\"ufe den Grenzwertsatz in gleicher Weise mit exponentiell
- verteilten Zufallszahlen (Funktion \code{rande}).
+
+ \part \extra \"Uberpr\"ufe den Grenzwertsatz in gleicher Weise mit
+ exponentiell verteilten Zufallszahlen (Funktion \code{rande}).
\end{parts}
\begin{solution}
\lstinputlisting{centrallimit.m}
@@ -169,56 +192,81 @@ Den Zentralen Grenzwertsatz wollen wir uns im Folgenden veranschaulichen.
\end{solution}
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-\question \qt{Random Walk}
-Im folgenden wollen wir einige Eigenschaften des Random Walks bestimmen.
-\begin{parts}
- \part Schreibe eine Funktion, die einen einzelnen Random Walk mit
- Startwert 0 f\"ur $n$ Schritte und Wahrscheinlichkeit $p$ f\"ur
- einen positiven Schritt als Vektor zur\"uckgibt.
- \part Visualisiere jeweils 10 Random Walks mit $p=0.5$ zusammen in einem Plot
- f\"ur $n=100$, $n=1000$ und $n=10000$ (drei Plots).\\
- Sch\"atze aus den Abbildungen ab, wie sich der Mittelwert und die Standardabweichung
- des Random Walks mit der Zeit (Schritte) sich entwickelt.
- \part \"Uberpr\"uefe deine Hypothese zum Mittelwert und zur
- Standardabweichung, indem du von $m$ Random Walks ($m \ge 10$) f\"ur
- jeden z.B. zehnten Schritt den Mittelwert und die Standardabweichung
- \"uber die Positionen der $m$ Random Walks berechnest.\\
- Wie h\"angt also die Standardabweichung von der Anzahl der Schritte
- ab? Wie entwickelt sich die Standardabweichung f\"ur eine sehr
- gro{\ss}e Anzahl von Schritten?
- \part \extra Erstelle eine Grafik, die die Verteilung der Position eines Random Walkers
- zu drei verschiedenen Zeitpunkten zeigt.
-\end{parts}
-\begin{solution}
- \lstinputlisting{randomwalk.m}
- \lstinputlisting{randomwalkstatistics.m}
- \includegraphics[width=0.8\textwidth]{randomwalk-traces}\\
- \includegraphics[width=0.5\textwidth]{randomwalk-stdev}
- \includegraphics[width=0.5\textwidth]{randomwalk-hists}
-\end{solution}
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+ \question \qt{Intervallstatistik von Spiketrains}
+ In Ilias findet ihr die Dateien \code{poisson.mat},
+ \code{pifou.mat}, und \code{lifadapt.mat}. Jede dieser Dateien
+ enth\"alt mehrere Trials von Spiketrains von einer bestimmten Art
+ von Neuron. Die Spikezeiten sind in Sekunden gemessen.
+ Mit den folgenden Aufgaben wollen wir die Intervallstatistik der
+ Spiketrains der drei Neurone miteinander vergleichen.
+ \begin{parts}
+ \part Lade die Spiketrains aus den drei Dateien. Achte darauf,
+ dass sie verschiedene Variablen\-namen bekommen.
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-\question \qt{\extra 2D Random Walk}
-Bisher hat sich unser Random Walker nur in einer Dimension bewegt
-(nur vorw\"arts oder r\"uckw\"arts). Er kann aber auch in mehreren Dimensionen laufen!\\
-In zwei Dimensionen wird dazu in jedem Schritt eine weitere
-Zufallszahl gezogen, die bestimmt ob er einen Schritt nach links oder
-rechts gemacht hat. Die Bewegung nach vorne/hinten bzw. links/rechts
-sind unabh\"angig voneinander.
-\begin{parts}
- \part Wie kann unter Verwendung unserer Funktion f\"ur den
- eindimensionalen Random Walk ein zweidimensionaler Random Walk
- simuliert werden?
- \part Erstelle h\"ubsche Bilder, die zweidimensionalen Random
- Walks verschiedener L\"ange (bis zu mindestens $n=1000000$) illustrieren.
- \part Animationen sind auch sch\"on! z.B. mit dem \code{pause} Befehl.
- \part Anstatt einfach den Weg des Random Walks zu zeichnen, kann man
- sich auch merken, wie oft er an jeder Stelle vorbeigekommen ist und
- mit einem Farbcode plotten.
-\end{parts}
+ In welchem Datentyp liegen die Daten vor? Wie kann auf einzelne
+ Spiketrains zugegriffen werden? Wie auf einzelne Spikezeiten?
+ \begin{solution}
+ \begin{lstlisting}
+ clear all
+ % not so good:
+ load poisson.mat
+ whos
+ poissonspikes = spikes;
+ load pifou.mat;
+ pifouspikes = spikes;
+ load lifadapt.mat;
+ lifadaptspikes = spikes;
+ clear spikes;
+ % better:
+ clear all
+ x = load( 'poisson.mat' );
+ poissonspikes = x.spikes;
+ x = load( pifou.mat' );
+ pifouspikes = x.spikes;
+ x = load( 'lifadapt.mat' );
+ lifadaptspikes = x.spikes;
+ \end{lstlisting}
+ \end{solution}
+
+ \part Schreibe eine Funktion, die die Spikezeiten der ersten
+ $t_{max}$ Sekunden in einem Rasterplot visualisiert. In jeder
+ Zeile des Rasterplots wird ein Spiketrain dargestellt. Jeder
+ einzelne Spike wird als senkrechte Linie zu der Zeit des
+ Auftretens des Spikes geplottet. Benutze die Funktion, um die
+ Spikeraster der ersten 1\,s der drei Neurone nebeneinander zu plotten.
+ \begin{solution}
+ \lstinputlisting{../../pointprocesses/code/spikeraster.m}
+ \lstinputlisting{../../pointprocesses/code/plotspikeraster.m}
+ \mbox{}\\[-3ex]
+ \colorbox{white}{\includegraphics[width=1\textwidth]{spikeraster}}
+ \end{solution}
+
+ \part Schreibe eine Funktion, die einen einzigen Vektor mit den
+ Interspikeintervallen aller Trials von Spikezeiten zur\"uckgibt.
+ \begin{solution}
+ \lstinputlisting{../../pointprocesses/code/isis.m}
+ \end{solution}
+
+ \part Schreibe eine Funktion, die ein normiertes Histogramm aus
+ einem Vektor von Interspikeintervallen, gegeben in Sekunden,
+ berechnet und dieses mit richtiger Achsenbeschriftung plottet.
+ Die Interspikeintervalle sollen dabei in Millisekunden angegeben
+ werden. Die Funktion soll zus\"atzlich den Mittelwert, die
+ Standardabweichung, und den Variationskoeffizienten der
+ Interspikeintervalle berechnen und diese im Plot mit angeben.
+
+ Benutze die vorherige und diese Funktion, um die
+ Interspikeintervall Verteilung der drei Neurone zu vergleichen.
+ \begin{solution}
+ \lstinputlisting{../../pointprocesses/code/isihist.m}
+ \lstinputlisting{../../pointprocesses/code/plotisih.m}
+ \mbox{}\\[-3ex]
+ \colorbox{white}{\includegraphics[width=1\textwidth]{isihist}}
+ \end{solution}
+ \end{parts}
\end{questions}
-\end{document}
\ No newline at end of file
+\end{document}
diff --git a/statistics/exercises/normprobs.m b/statistics/exercises/normprobs.m
index 33af699..95f1d0c 100644
--- a/statistics/exercises/normprobs.m
+++ b/statistics/exercises/normprobs.m
@@ -1,14 +1,30 @@
%% (a) generate normal distributed random numbers:
n = 10000;
-x = randn( n, 1 );
+x = randn(n, 1);
-%% (b)
+%% (b) plot histogram and compare with pdf:
+bw = 0.5;
+bins=[-5:bw:5];
+n = hist(x, bins);
+p = n/sum(n)/bw;
+subplot(1, 2, 2);
+bar(bins, p);
+hold on;
+xx = [bins(1):0.01:bins(end)];
+gauss = exp(-0.5*xx.^2.0)/sqrt(2*pi);
+plot(xx, gauss, 'r', 'linewidth', 2);
+hold off;
+xlim([-5 5])
+xlabel('x');
+ylabel('p(x)');
+
+%% (c) fraction of data in +/- 1 sigma:
nsig = sum((x>=-1.0)&(x<=1.0));
Psig = nsig/length(x);
fprintf('%d of %d data elements, i.e. %.2f%% are contained in the interval -1 to +1\n\n', ...
nsig, length(x), 100.0*Psig );
-%% (c)
+%% (d) intgegral over pdf:
dx = 0.01;
xx = -10:dx:10; % x-values
pg = exp(-0.5*xx.^2)/sqrt(2*pi); % y-values Gaussian pdf
@@ -17,10 +33,25 @@ P = sum(pg)*dx;
fprintf( 'Integral over the Gaussian pdf is %.3f\n', P );
% integral from -1 to 1:
P = sum(pg((xx>=-1.0)&(xx<=1.0)))*dx; % we need to use xx, not the random numbers x!
-fprintf( 'Integral over the Gaussian pdf from -1 to 1 is %.4f\n', P );
+fprintf( 'Integral over the Gaussian pdf from -1 to 1 is %.4f\n\n', P );
+
+%% (e)
+for sigma = [1.0, 2.0, 3.0]
+ Pdata = sum((x>=-sigma)&(x<=sigma))/length(x);
+ Ppdf = sum(pg((xx>=-sigma)&(xx<=sigma)))*dx;
+ fprintf( 'Integral over the Gaussian pdf from -%.0f to +%.0f is %.4f\n', sigma, sigma, Ppdf );
+ fprintf( 'Probability of data in range from -%.0f to +%.0f is %.4f\n\n', sigma, sigma, Pdata );
+end
+
+%% (f)
+for P = [0.5, 0.9, 0.95, 0.99]
+ for upper = xx(xx>0.0)
+ Ppdf = sum(pg((xx>=-upper)&(xx<=upper)))*dx;
+ if Ppdf > P
+ fprintf('-%.2f < x < +%.2f: P=%.2f\n', upper, upper, P);
+ break
+ end
+ end
+end
+
-%% (d)
-P = sum(pg((xx>=-2.0)&(xx<=2.0)))*dx;
-fprintf( 'Integral over the Gaussian pdf from -2 to 2 is %.4f\n', P );
-P = sum(pg((xx>=-3.0)&(xx<=3.0)))*dx;
-fprintf( 'Integral over the Gaussian pdf from -3 to 3 is %.4f\n\n', P );