diff --git a/Makefile b/Makefile index 74de0a4..aa43ca0 100644 --- a/Makefile +++ b/Makefile @@ -2,7 +2,7 @@ BASENAME=scientificcomputing-script pdf : $(BASENAME).pdf -$(BASENAME).pdf : $(BASENAME).tex +$(BASENAME).pdf : $(BASENAME).tex header.tex export TEXMFOUTPUT=.; pdflatex -interaction=scrollmode $< | tee /dev/stderr | fgrep -q "Rerun to get cross-references right" && pdflatex -interaction=scrollmode $< || true clean : diff --git a/bootstrap/code/bootstrapsem.out b/bootstrap/code/bootstrapsem.out new file mode 100644 index 0000000..ddb7462 --- /dev/null +++ b/bootstrap/code/bootstrapsem.out @@ -0,0 +1,5 @@ +>> bootstrapsem + mean stdev sem + single SRS -0.13 0.97 0.10 + bootstrapped distribution -0.13 0.10 - + sampling distribution -0.00 0.10 - diff --git a/bootstrap/code/correlationsignificance.m b/bootstrap/code/correlationsignificance.m new file mode 100644 index 0000000..f5ccb1c --- /dev/null +++ b/bootstrap/code/correlationsignificance.m @@ -0,0 +1,38 @@ +% generate correlated data: +n=200; +a=0.2; +x = randn(n, 1); +y = randn(n, 1) + a*x; + +% correlation coefficient: +rd = corr(x, y); +fprintf('correlation coefficient of data r = %.2f\n', rd ); + +% distribution of null hypothesis by permutation: +nperm = 1000; +rs = zeros(nperm,1); +for i=1:nperm + xr=x(randperm(length(x))); % shuffle x + yr=y(randperm(length(y))); % shuffle y + rs(i) = corr(xr, yr); +end +[h,b] = hist(rs, 20 ); +h = h/sum(h)/(b(2)-b(1)); % normalization + +% significance: +rq = quantile(rs, 0.95); +fprintf('correlation coefficient of null hypothesis at 5%% significance = %.2f\n', rq ); +if rd >= rq + fprintf('--> correlation r=%.2f is significant\n', rd); +else + fprintf('--> r=%.2f is not a significant correlation\n', rd); +end + +% plot: +bar(b, h, 'facecolor', 'b'); +hold on; +bar(b(b>=rq), h(b>=rq), 'facecolor', 'r'); +plot( [rd rd], [0 4], 'r', 'linewidth', 2 ); +xlabel('Correlation coefficient'); +ylabel('Probability density of H0'); +hold off; diff --git a/bootstrap/code/correlationsignificance.out b/bootstrap/code/correlationsignificance.out new file mode 100644 index 0000000..b3d84b5 --- /dev/null +++ b/bootstrap/code/correlationsignificance.out @@ -0,0 +1,4 @@ +>> correlationsignificance +correlation coefficient of data r = 0.23 +correlation coefficient of null hypothesis at 5% significance = 0.12 +--> correlation r=0.23 is significant \ No newline at end of file diff --git a/bootstrap/lecture/bootstrap.tex b/bootstrap/lecture/bootstrap.tex index e5d63aa..5b62f89 100644 --- a/bootstrap/lecture/bootstrap.tex +++ b/bootstrap/lecture/bootstrap.tex @@ -5,32 +5,51 @@ Beim Bootstrap erzeugt man sich die Verteilung von Statistiken durch Resampling aus der Stichprobe. Das hat mehrere Vorteile: \begin{itemize} -\item Weniger Annahmen (z.B. muss eine Stichprobe nicht Normalverteilt sein). +\item Weniger Annahmen (z.B. muss eine Stichprobe nicht normalverteilt sein). \item H\"ohere Genauigkeit als klassische Methoden. \item Allgemeing\"ultigkeit: Bootstrap Methoden sind sich sehr \"ahnlich f\"ur viele verschiedene Statistiken und ben\"otigen nicht f\"ur jede Statistik eine andere Formel. \end{itemize} -\begin{figure}[t] +\begin{figure}[tp] \includegraphics[width=0.8\textwidth]{2012-10-29_16-26-05_771}\\[2ex] \includegraphics[width=0.8\textwidth]{2012-10-29_16-41-39_523}\\[2ex] \includegraphics[width=0.8\textwidth]{2012-10-29_16-29-35_312} - \caption{\tr{Why can we only measure a sample of the + \caption{\label{statisticalpopulationfig}\tr{Why can we only measure a sample of the population?}{Warum k\"onnen wir nur eine Stichprobe der Grundgesamtheit messen?}} \end{figure} -\begin{figure}[t] +Zur Erinnerung: In der Statistik interessieren wir uns f\"ur +Eigenschaften einer Grundgesamtheit. z.B. die mittlere L\"ange von +sauren Gurken (\figref{statisticalpopulationfig}). Aus der +Grundgesamtheit wird eine Stichprobe (simple random sample, SRS) +gezogen, da niemals die gesamte Grundgesamtheit gemessen werden kann. +Dann wird aus dieser einzigen Stichprobe die gew\"unschte Gr\"o{\ss}e +berechnet (die mittlere Gr\"o{\ss}e der sauren Gurken) und man hofft, +dass die erhaltene Zahl an der entsprechenden unbekannten Gr\"o{\ss}e +der Grundgesamtheit (der Populationsparameter) m\"oglichst nah dran +ist. Eine Aufgabe der Statistik ist es, herauszubekommen wie gut der +Populationsparameter abgesch\"atzt worden ist. + +Wenn wir viele Stichproben ziehen w\"urden, dann k\"onnte man f\"ur +jede Stichprobe den gew\"unschten Parameter berechnen, und von diesen +die Wahrscheinlichkeitsverteilung \"uber ein Histogramm bestimmen --- +die ``Stichprobenverteilung'' (sampling distribution, +\subfigref{bootstrapsamplingdistributionfig}{a}). + +\begin{figure}[tp] \includegraphics[height=0.2\textheight]{srs1}\\[2ex] \includegraphics[height=0.2\textheight]{srs2}\\[2ex] \includegraphics[height=0.2\textheight]{srs3} - \caption{Bootstrap der Stichprobenvertielung (a) Von der - Grundgesamtheit (population) mit unbekanntem Parameter - (z.B. Mittelwert $\mu$) zieht man Stichproben (SRS: simple random - samples). Die Statistik (hier Bestimmung von $\bar x$) kann f\"ur - jede Stichprobe berechnet werden. Die erhaltenen Werte entstammen - der Stichprobenverteilung. Meisten wird aber nur eine Stichprobe + \caption{\label{bootstrapsamplingdistributionfig}Bootstrap der + Stichprobenvertielung (a) Von der Grundgesamtheit (population) mit + unbekanntem Parameter (z.B. Mittelwert $\mu$) zieht man + Stichproben (SRS: simple random samples). Die Statistik (hier + Bestimmung von $\bar x$) kann f\"ur jede Stichprobe berechnet + werden. Die erhaltenen Werte entstammen der + Stichprobenverteilung. Meisten wird aber nur eine Stichprobe gezogen! (b) Mit bestimmten Annahmen und Theorien kann man auf die Stichprobenverteilung schlie{\ss}en ohne sie gemessen zu haben. (c) Alternativ k\"onnen aus der einen Stichprobe viele @@ -40,25 +59,138 @@ aus der Stichprobe. Das hat mehrere Vorteile: Permuation Tests} \end{figure} -\section{Bootstrap des Standardfehlers} +In Wirklichkeit haben wir aber nur eine Stichprobe. Wir behelfen uns +dann mit Theorien, die meistens bestimmte Annahmen \"uber die Daten +machen (z.B. Normalverteilung), und uns erlauben etwas \"uber die +Genaugigkeit unserer Sch\"atzung aus der Stichprobe auszusagen +(z.B. die Formel $\sigma/\sqrt{n}$ f\"ur den Standardfehler des +Mittelwerts, die uns die Standardabweichung angibt, mit dem die +Mittelwerte der Stichproben um den Populationsmittelwert streuen +\subfigref{bootstrapsamplingdistributionfig}{b}). + +Wir k\"onnen aber auch aus der einen Stichprobe die wir haben durch +Resampling viele neue Stichproben generieren (Bootstrap). Von diesen +k\"onnen wir jeweils die gew\"unschte Gr\"o{\ss}e berechnen und ihre +Verteilung bestimmen (Bootstrap Verteilung, +\subfigref{bootstrapsamplingdistributionfig}{c}). Diese Verteilung ist +interessanterweise in ihrer Breite und Form der Stichprobenverteilung +sehr \"ahnlich. Nur streut sie nicht um den Populationswert sonder um +die Sch\"atzung aus der Stichprobe. Wir k\"onnen die +Bootstrapverteilung aber benutzen um Aussagen \"uber die Genauigkeit +unserer Sch\"atzung zu treffen (z.B. Standardfehler, +Konfidenzintervalle). Beim Bootstrap erzeugen wir durch Resampling neue Stichproben und -benutzen diese um die Stichprobenverteilung einer Statistik zu +benutzen diese, um die Stichprobenverteilung einer Statistik zu berechnen. Die Bootstrap Stichproben haben jeweils den gleichen Umfang wie die urspr\"unglich gemessene Stichprobe und werden durch Ziehen mit Zur\"ucklegen gewonnen. Jeder Wert der urspr\"unglichen Stichprobe kann also einmal, mehrmals oder gar nicht in einer Bootstrap Stichprobe vorkommen. -\begin{exercise}[bootstrapsem.m] - Ziehe 1000 normalverteilte Zufallszahlen und berechne deren Mittelwert, - Standardabweichung und Standardfehler ($\sigma/\sqrt{n}$). - Resample die Daten 1000 mal (Ziehen mit Zur\"ucklegen) und berechne jeweils - den Mittelwert. +\section{Bootstrap des Standardfehlers} + +Am besten l\"asst sich die Bootstrap Methode am Beispiel des +Standardfehlers des Mittelwertes veranschaulichen. Aus der Stichprobe +k\"onnen wir den Mittelwert berechnen. Der Standardfehler des +Mittelwerts gibt die Standardabweichung an, mit der wir erwarten, dass +der gemessene Mittelwert um den Populationsmittelwert streut. + +\begin{figure}[tp] + \includegraphics[width=1\textwidth]{bootstrapsem} + \caption{\label{bootstrapsemfig}Bootstrap des Standardfehlers des + Mittelwertes. Die --- normalerweise unbekannte --- + Stichprobenverteilung des Mittelwerts (rot) ist um den + Populationsmittelwert bei $\mu=0$ zentriert. Die + Bootstrap-Verteilung (blau) die durch Resampling aus einer + Stichprobe gewonnen worden ist hat die gleiche Form und Breite wie + die Stichprobenverteilung, ist aber um den Mittelwert berechnet + aus der Stichprobe zentriert. Die Standardabweichung der + Bootstrapverteilung kann also als Sch\"atzer f\"ur den + Standardfehler des Mittelwertes verwendet werden.} +\end{figure} - Plotte ein Histogramm dieser Mittelwerte, sowie deren Mittelwert und - die Standardabweichung. +Durch Bootstrap k\"onnen wir unsere Stichprobe resamplen und dadurch +eine ganze Verteilung von Mittelwerten generieren +(\figref{bootstrapsemfig}). Die Standardabweichung dieser Verteilung +ist dann der gesuchte Standardfehler des Mittelwerts. - Was hat das mit dem Standardfehler zu tun? +\begin{exercise}{bootstrapsem.m}{bootstrapsem.out} + Erzeuge die Verteilung der Mittelwerte einer Stichprobe durch Bottstrapping, + um daraus den Standardfehler des Mittelwerts zu bestimmen. + \begin{enumerate} + \item Ziehe 1000 normalverteilte Zufallszahlen und berechne deren + Mittelwert, Standardabweichung und Standardfehler + ($\sigma/\sqrt{n}$). + \item Resample die Daten 1000 mal (Ziehen mit Zur\"ucklegen) und + berechne jeweils den Mittelwert. + \item Plotte ein Histogramm dieser Mittelwerte, berechne deren + Mittelwert und Standardabweichung und vergleiche mit den Werten + der Grundgesamtheit und der Stichprobe. + \end{enumerate} \end{exercise} + + +\section{Permutationstests} +Bei statistischen Tests wird nach der Wahrscheinlichkeit, ob die +beobachtete Me{\ss}gr\"o{\ss}e einer Stichprobe aus der Nullhypothese +kommt, gefragt. Ist diese Wahrscheinlichkeit kleiner als das +Signifikanzniveau, kann die Nullhypothese verworfen werden. + +Traditionell werden diese Wahrscheinlichkeiten \"uber theoretisch +hergeleitete Wahrscheinlichkeitsverteilungen berechnet. Dabei gehen +immer gewisse Annahmen \"uber die Daten ein und es mu{\ss} der zu den +Daten passende Test ausgew\"ahlt werden. + +Alternativ kann die Wahrscheinlichkeits(dichte)verteilung der +Nullhypothese aus den Daten selbst gewonnen werden. Dabei m\"ussen die +Daten entsprechend der Nullhypothese neu aus der Stichprobe gezogen +werden. + +Diese ``Permutationstests'' haben den Vorteil, dass nur die +Eigenschaft von Interesse zerst\"ort wird, um die Nullhypothese zu +generieren. Alle anderen Eigenschaften der Daten bleiben erhalten. + +\begin{figure}[tp] + \includegraphics[width=1\textwidth]{permutecorrelation} + \caption{\label{permutecorrelationfig}Permutationstest f\"ur + Korrelationen. Der Korrelationskoeffizient eines Datensatzes mit + 200 Datenpaaren ist $\rho=0.21$. Die Nullhypothesenverteilung der + aus den permutierten, unkorrelierten Datens\"atzen berechneten + Korrelationskoeffizienten ergibt die gelbe Verteilung, die um Null + streut. Der gemessene Korrelationskoeffizient ist deutlich + gr\"o{\ss}er als das 95\,\%-Perzentil der + Nullhypoothesenverteilung und darum eine signifikante + Korrelation.} +\end{figure} + +Sehr sch\"on lassen sich Permutationstest am Beispiel von +Korrelationen veranschaulichen. Gegeben sind Datenpaare $(x_i, y_i)$. +Daraus k\"onnen wir den Korrelationskoeffizienten berechnen. Wir +wissen dann aber noch nicht, ob der berechnete Wert tats\"achlich eine +Korrelation anzeigt. Die Nullhypothese ist, dass die Daten nicht +miteinander korreliert sind. Indem wir die $x$-Werte und die $y$-Werte +unabh\"angig voneinander permutieren (ihre Reihenfolge zuf\"allig neu +anordnen), werden die Korrelationen der Datenpaare zerst\"ort. Wenn +wir das viele Male wiederholen, bekommen wir die Verteilung der +Korrelationskoeffizienten f\"ur nichtkorrelierte Daten. Aus dieser +Verteilung der Nullhypothese k\"onnen wir dann dann die Signifikanz +der tats\"achlich gemessenen Korrelation bestimmen. + +\begin{exercise}{correlationsignificance.m}{correlationsignificance.out} +Bestimme die Signifikanz eines Korrelationskoeffizienten. +\begin{enumerate} +\item Erzeuge korrelierte Daten indem zu zuf\"allig gezogenen + $x$-Werten $y$-Werte gem\"a{\ss} $y=0.2 \cdot x$ berechnet werden, + zu denen weitere normalverteilte Zufallszahlen addiert werden. +\item Berechne den Korrelationskoeffizient dieser Datenpaare. +\item Generiere die Verteilung der Nullhypothese ``unkorrelierte + Daten'' indem die $x$- und $y$-Daten 1000-mal unabh\"angig + permutiert werden \matlabfun{randperm} und jeweils der + Korrelationskoeffizient berechnet wird. +\item Bestimme aus den Nullhypothesendaten das 95\,\%-Perzentil und + vergleiche es mit dem tats\"achlichen Korrelationskoeffizienten. +\end{enumerate} +\end{exercise} + diff --git a/bootstrap/lecture/bootstrapsem.py b/bootstrap/lecture/bootstrapsem.py new file mode 100644 index 0000000..7b9ac22 --- /dev/null +++ b/bootstrap/lecture/bootstrapsem.py @@ -0,0 +1,53 @@ +import numpy as np +import matplotlib.pyplot as plt + +plt.xkcd() +fig = plt.figure( figsize=(6,4) ) +rng = np.random.RandomState(637281) + +nsamples = 100 +nresamples = 1000 +db = 0.05; +bins = np.arange(-0.4, 0.41, db) + +# draw a SRS (simple random sample, "Stichprobe") from the population: +x = rng.randn(nsamples) + +# bootstrap the mean: +mus = [] +for i in xrange(nresamples) : + mus.append(np.mean(x[rng.randint(0, nsamples, nsamples)])) +hmus, _ = np.histogram(mus, bins, density=True) + +# many SRS: +musrs = [] +for i in xrange(nresamples) : + musrs.append(np.mean(rng.randn(nsamples))) +hmusrs, _ = np.histogram(musrs, bins, density=True) + +ax = fig.add_subplot(1, 1, 1) +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_xlabel('Mean') +ax.set_xlim(-0.4, 0.4) +ax.set_ylabel('Probability density') +ax.set_ylim(0.0, 4.5) +ax.set_yticks(np.arange(0.0, 4.5, 1.0)) +ax.annotate('sampling\ndistribution', + xy=(-0.12, 3.1), xycoords='data', + xytext=(-0.18, 3.5), textcoords='data', ha='right', + arrowprops=dict(arrowstyle="->", relpos=(1.0,0.5), + connectionstyle="angle3,angleA=20,angleB=120") ) +ax.annotate('bootstrap\ndistribution', + xy=(0.13, 3.3), xycoords='data', + xytext=(0.25, 4), textcoords='data', + arrowprops=dict(arrowstyle="->", relpos=(0.0,0.5), + connectionstyle="angle3,angleA=20,angleB=60") ) +ax.bar(bins[:-1]-0.25*db, hmusrs, 0.5*db, color='r') +ax.bar(bins[:-1]+0.25*db, hmus, 0.5*db, color='b') + +plt.tight_layout() +plt.savefig('bootstrapsem.pdf') +#plt.show(); diff --git a/bootstrap/lecture/permutecorrelation.py b/bootstrap/lecture/permutecorrelation.py new file mode 100644 index 0000000..580e12b --- /dev/null +++ b/bootstrap/lecture/permutecorrelation.py @@ -0,0 +1,60 @@ +import numpy as np +import matplotlib.pyplot as plt + +plt.xkcd() +fig = plt.figure( figsize=(6,4) ) +rng = np.random.RandomState(637281) + +# generate correlated data: +n = 200 +a = 0.2 +x = rng.randn(n); +y = rng.randn(n) + a*x; + +rd = np.corrcoef(x, y)[0, 1] +print rd + +# permutation: +nperm = 1000 +rs = [] +for i in xrange(nperm) : + xr=rng.permutation(x) + yr=rng.permutation(y) + rs.append( np.corrcoef(xr, yr)[0, 1] ) + +# pdf of the correlation coefficients: +h, b = np.histogram(rs, 20, density=True) + +# significance: +rq = np.percentile(rs, 95.0); + +ax = fig.add_subplot(1, 1, 1) +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.annotate('Measured\ncorrelation\nis significant!', + xy=(rd, 1.1), xycoords='data', + xytext=(rd, 2.2), textcoords='data', ha='left', + arrowprops=dict(arrowstyle="->", relpos=(0.2,0.0), + connectionstyle="angle3,angleA=10,angleB=80") ) +ax.annotate('95% percentile', + xy=(0.14, 0.9), xycoords='data', + xytext=(0.2, 4.0), textcoords='data', ha='left', + arrowprops=dict(arrowstyle="->", relpos=(0.1,0.0), + connectionstyle="angle3,angleA=30,angleB=70") ) +ax.annotate('Distribution of\nuncorrelated\nsamples', + xy=(-0.08, 3.6), xycoords='data', + xytext=(-0.22, 5.0), textcoords='data', ha='left', + arrowprops=dict(arrowstyle="->", relpos=(0.5,0.0), + connectionstyle="angle3,angleA=150,angleB=100") ) +ax.bar(b[:-1], h, width=b[1]-b[0], color='#ffff66') +ax.bar(b[:-1][b[:-1]>=rq], h[b[:-1]>=rq], width=b[1]-b[0], color='#ff9900') +ax.plot( [rd, rd], [0, 1], 'b', linewidth=4 ) +ax.set_xlim(-0.25, 0.35) +ax.set_xlabel('Correlation coefficient') +ax.set_ylabel('Probability density of H0') + +plt.tight_layout() +plt.savefig('permutecorrelation.pdf') +#plt.show(); diff --git a/header.tex b/header.tex index 401a3ff..131b443 100644 --- a/header.tex +++ b/header.tex @@ -1,8 +1,8 @@ %%%%% title %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -\title{\tr{Introduction to Scientific Computing}{Einf\"uhrung in die wissenschaftliche Datenverarbeitung}} -\author{Jan Benda, Jan Grewe\\Abteilung Neuroethologie\\[2ex]\includegraphics[width=0.3\textwidth]{UT_WBMW_Rot_RGB}} -\date{WS 15/16} +\title{\textbf{\huge\sffamily\tr{Introduction to\\[1ex] Scientific Computing}{Einf\"uhrung in die\\[1ex] wissenschaftliche Datenverarbeitung}}} +\author{{\LARGE Jan Grewe \& Jan Benda}\\[4ex]Abteilung Neuroethologie\\[2ex]\includegraphics[width=0.3\textwidth]{UT_WBMW_Rot_RGB}} +\date{WS 15/16\\\vfill \includegraphics[width=1\textwidth]{announcements/correlationcartoon}} %%%% language %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % \newcommand{\tr}[2]{#1} % en @@ -139,7 +139,7 @@ backgroundcolor=\color{blue!10}, breaklines=true, breakautoindent=true, - columns=flexible, +% columns=flexible, % this kills the monospace font frame=single, caption={\protect\filename@parse{\lstname}\protect\filename@base}, captionpos=t, @@ -172,19 +172,38 @@ {\medskip} %%%%% exercises environment: %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% usage: +% +% \begin{exercise}[somecode.m] +% Write a function that computes the median. +% \end{exercise} +% +% result: +% +% Exercise 1: +% Write a function that computes the median. +% Listing 1: somecode.m +% the listing +% +% Innerhalb der exercise Umgebung ist enumerate umdefiniert, um (a), (b), (c), .. zu erzeugen. \usepackage{framed} \newcounter{maxexercise} \setcounter{maxexercise}{10000} % show listings up to exercise maxexercise \newcounter{theexercise} \setcounter{theexercise}{1} \newcommand{\codepath}{} -\newenvironment{exercise}[1][]% +\newenvironment{exercise}[2]% {\newcommand{\exercisesource}{#1}% + \newcommand{\exerciseoutput}{#2}% \setlength{\fboxsep}{2mm}% + \newcommand{\saveenumi}{\theenumi}\renewcommand{\labelenumi}{(\alph{enumi})}% \renewcommand{\FrameCommand}{\colorbox{yellow!15}}% \MakeFramed{\advance\hsize-\width \FrameRestore}% \noindent\textbf{\tr{Exercise}{\"Ubung} \arabic{theexercise}:}\newline}% {\ifthenelse{\equal{\exercisesource}{}}{}% {\ifthenelse{\value{theexercise}>\value{maxexercise}}{}% - {\lstinputlisting[belowskip=0pt]{\codepath\exercisesource}}}% - \endMakeFramed\stepcounter{theexercise}} + {\lstinputlisting[belowskip=0pt]{\codepath\exercisesource}% + \ifthenelse{\equal{\exerciseoutput}{}}{}% + {\addtocounter{lstlisting}{-1}\lstinputlisting[language={},title={\textbf{\tr{Output}{Ausgabe}:}},belowskip=0pt]{\codepath\exerciseoutput}}}}% + \endMakeFramed\stepcounter{theexercise}% + \renewcommand{\theenumi}{\saveenumi}} diff --git a/likelihood/code/mlegammafit.out b/likelihood/code/mlegammafit.out new file mode 100644 index 0000000..45176e0 --- /dev/null +++ b/likelihood/code/mlegammafit.out @@ -0,0 +1,3 @@ +>> mlegammafit +shape=2.10 +scale=0.95 diff --git a/likelihood/code/mlemean.out b/likelihood/code/mlemean.out new file mode 100644 index 0000000..6a9feaf --- /dev/null +++ b/likelihood/code/mlemean.out @@ -0,0 +1,5 @@ +>> mlemean +standard deviation of the data is 1.93 + mean of the data is 2.98 + maximum of likelihood is at 2.98 + maximum of log-likelihood is at 2.98 diff --git a/likelihood/lecture/likelihood.tex b/likelihood/lecture/likelihood.tex index 5e73626..d806e7f 100644 --- a/likelihood/lecture/likelihood.tex +++ b/likelihood/lecture/likelihood.tex @@ -105,7 +105,7 @@ x$ der Daten. D.h. das arithmetische Mittel maximiert die Wahrscheinlichkeit, dass die Daten aus einer Normalverteilung mit diesem Mittelwert gezogen worden sind (\figref{mlemeanfig}). -\begin{exercise}[mlemean.m] +\begin{exercise}{mlemean.m}{mlemean.out} Ziehe $n=50$ normalverteilte Zufallsvariablen mit einem Mittelwert $\ne 0$ und einer Standardabweichung $\ne 1$. @@ -244,7 +244,7 @@ gesuchten Wahrscheinlichkeitsdichtefunktion bei der die Log-Likelihood nichtlinieares Optimierungsproblem, das mit numerischen Verfahren, wie z.B. dem Gradientenabstieg, gel\"ost wird \matlabfun{mle}. -\begin{exercise}[mlegammafit.m] +\begin{exercise}{mlegammafit.m}{mlegammafit.out} Erzeuge Gammaverteilte Zufallszahlen und benutze Maximum-Likelihood, um die Parameter der Gammafunktion aus den Daten zu bestimmen. \newpage diff --git a/likelihood/lecture/mlecoding.py b/likelihood/lecture/mlecoding.py index ba86bfd..ceac2a1 100644 --- a/likelihood/lecture/mlecoding.py +++ b/likelihood/lecture/mlecoding.py @@ -94,7 +94,7 @@ ax.annotate('', xytext=(maxp, -30), textcoords='data', arrowprops=dict(arrowstyle="->", relpos=(0.5,0.5), connectionstyle="angle3,angleA=80,angleB=90") ) -ax.text(maxp+0.05, -1300, 'most\nlikely\norientation') +ax.text(maxp+0.05, -1100, 'most likely\norientation\ngiven the responses') ax.plot(phases, loglikelihood, '-b') plt.savefig('mlecoding.pdf') diff --git a/pointprocesses/lecture/pointprocesses.tex b/pointprocesses/lecture/pointprocesses.tex index 4777719..de9c1b6 100644 --- a/pointprocesses/lecture/pointprocesses.tex +++ b/pointprocesses/lecture/pointprocesses.tex @@ -38,22 +38,49 @@ erzeugt. Zum Beispiel: Zeitkonstante $\tau=100$\,ms, rechts).} \end{figure} +F\"ur die Neurowissenschaften ist die Statistik der Punktprozesse +besonders wichtig, da die Zeitpunkte der Aktionspotentiale als +zeitlicher Punktprozess betrachtet werden k\"onnen und entscheidend +f\"ur die Informations\"ubertragung sind. + +Bei Punktprozessen k\"onnen wir die Zeitpunkte $t_i$ ihres Auftretens, +die Intervalle zwischen diesen Zeitpunkten $T_i=t_{i+1}-t_i$, sowie +die Anzahl der Ereignisse $n_i$ bis zu einer bestimmten Zeit betrachten +(\figref{pointprocessscetchfig}). + +Zwei Punktprozesse mit verschiedenen Eigenschaften sind in +\figref{rasterexamplesfig} als Rasterplot dargestellt, bei dem die +Zeitpunkte der Ereignisse durch senkrechte Striche markiert werden. + + %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -\section{Intervall Statistik} +\section{Intervallstatistik} + +Die Intervalle $T_i=t_{i+1}-t_i$ zwischen aufeinanderfolgenden +Ereignissen sind reelle, positive Zahlen. Bei Aktionspotentialen +(``Spikes'') heisen die Intervalle auch +``Interspikeintervalle''. Deren Statistik kann mit den \"ublichen +Gr\"o{\ss}en beschrieben werden. \begin{figure}[t] \includegraphics[width=1\textwidth]{isihexamples}\hfill - \caption{\label{isihexamplesfig}Interspike-Intervall Histogramme der in + \caption{\label{isihexamplesfig}Interspikeintervall Histogramme der in \figref{rasterexamplesfig} gezeigten Spikes.} \end{figure} -\subsection{(Interspike) Intervall Statistik erster Ordnung} +\subsection{Intervallstatistik erster Ordnung} \begin{itemize} -\item Histogramm $p(T)$ der Intervalle $T$. Normiert auf $\int_0^{\infty} p(T) \; dT = 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 Variationskoeffizient (``Coefficient of variation'') $CV_{ISI} = \frac{\sigma_{ISI}}{\mu_{ISI}}$. -\item Diffusions Koeffizient $D_{ISI} = \frac{\sigma_{ISI}^2}{2\mu_{ISI}^3}$. +\item Histogramm $p(T)$ der Intervalle $T$ + (\figref{isihexamplesfig}). Normiert auf $\int_0^{\infty} p(T) \; dT + = 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 Variationskoeffizient (``Coefficient of variation'') $CV_{ISI} = + \frac{\sigma_{ISI}}{\mu_{ISI}}$. +\item Diffusions Koeffizient $D_{ISI} = + \frac{\sigma_{ISI}^2}{2\mu_{ISI}^3}$. \end{itemize} \subsection{Korrelationen der Intervalle} @@ -65,7 +92,7 @@ sichtbar. \begin{figure}[t] \includegraphics[width=1\textwidth]{returnmapexamples} \includegraphics[width=1\textwidth]{serialcorrexamples} - \caption{\label{returnmapfig}Interspike-Intervall return maps und + \caption{\label{returnmapfig}Interspikeintervall return maps und serielle Korrelationen zwischen aufeinander folgenden Intervallen im Abstand des Lags $k$.} \end{figure} @@ -89,15 +116,17 @@ Intervalls mit sich selber). % \caption{\label{countstatsfig}Count Statistik.} % \end{figure} -Statistik der Anzahl der Ereignisse $N_i$ innerhalb von Beobachtungsfenstern $i$ der Breite $W$. +Die Anzahl der Ereignisse $n_i$ in Zeifenstern $i$ der +L\"ange $W$ ergeben ganzzahlige, positive Zufallsvariablen die meist durch folgende +Sch\"atzer charakterisiert werden: \begin{itemize} -\item Histogramm der counts $N_i$. -\item Mittlere Anzahl von Ereignissen: $\mu_N = \langle N \rangle$. -\item Varianz der Anzahl: $\sigma_N^2 = \langle (N - \langle N \rangle)^2 \rangle$. -\item Fano Faktor (Varianz geteilt durch Mittelwert): $F = \frac{\sigma_N^2}{\mu_N}$. +\item Histogramm der counts $n_i$. +\item Mittlere Anzahl von Ereignissen: $\mu_N = \langle n \rangle$. +\item Varianz der Anzahl: $\sigma_n^2 = \langle (n - \langle n \rangle)^2 \rangle$. +\item Fano Faktor (Varianz geteilt durch Mittelwert): $F = \frac{\sigma_n^2}{\mu_n}$. \end{itemize} Insbesondere ist die mittlere Rate der Ereignisse $r$ (``Spikes pro Zeit'', Feuerrate) gemessen in Hertz -\[ r = \frac{\langle N \rangle}{W} \; . \] +\[ r = \frac{\langle n \rangle}{W} \; . \] % \begin{figure}[t] % \begin{minipage}[t]{0.49\textwidth} @@ -116,25 +145,25 @@ Insbesondere ist die mittlere Rate der Ereignisse $r$ (``Spikes pro Zeit'', Feue %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \section{Homogener Poisson Prozess} F\"ur kontinuierliche Me{\ss}gr\"o{\ss}en ist die Normalverteilung -u.a. wegem dem Zentralen Grenzwertsatz die Standardverteilung. Eine +u.a. wegen dem Zentralen Grenzwertsatz die Standardverteilung. Eine \"ahnliche Rolle spilet bei Punktprozessen der ``Poisson Prozess''. Beim homogenen Poisson Prozess treten Ereignisse mit einer festen Rate $\lambda=\text{const.}$ auf und sind unabh\"angig von der Zeit $t$ und -unabh\"angig von den Zeitpunkten fr\"uherer Ereignisse. Die -Wahrscheinlichkeit zu irgendeiner Zeit ein Ereigniss in einem kleinen -Zeitfenster der Breite $\Delta t$ zu bekommen ist +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 \; . \] +Beim inhomogenen Poisson Prozess h\"angt die Rate $\lambda$ von der +Zeit ab: $\lambda = \lambda(t)$. \begin{figure}[t] \includegraphics[width=1\textwidth]{poissonraster100hz} - \caption{\label{hompoissonfig}Rasterplot von Spikes eine homogenen - Poisson Prozesse mit $\lambda=100$\,Hz.} + \caption{\label{hompoissonfig}Rasterplot von Spikes eines homogenen + Poisson Prozesses mit $\lambda=100$\,Hz.} \end{figure} -Beim inhomogenen Poisson Prozess h\"angt die Rate $\lambda$ von der -Zeit ab: $\lambda = \lambda(t)$. - \begin{figure}[t] \includegraphics[width=0.45\textwidth]{poissonisihexp20hz}\hfill \includegraphics[width=0.45\textwidth]{poissonisihexp100hz} @@ -144,16 +173,17 @@ Zeit ab: $\lambda = \lambda(t)$. Der homogne Poissonprozess hat folgende Eigenschaften: \begin{itemize} -\item Die Intervalle $T$ sind exponentiell verteilt: $p(T) = \lambda e^{-\lambda T}$ . +\item Die Intervalle $T$ sind exponentiell verteilt: $p(T) = \lambda e^{-\lambda T}$ (\figref{hompoissonisihfig}). \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$ . -\item Die seriellen Korrelationen $\rho_k =0$ for $k>0$, da das +\item Die seriellen Korrelationen $\rho_k =0$ f\"ur $k>0$, da das Auftreten der Ereignisse unabh\"angig von der Vorgeschichte ist. Ein solcher Prozess wird auch Erneuerungsprozess genannt (``renewal process''). \item Die Anzahl der Ereignisse $k$ innerhalb eines Fensters der L\"ange W ist Poissonverteilt: \[ P(k) = \frac{(\lambda W)^ke^{\lambda W}}{k!} \] + (\figref{hompoissoncountfig}) \item Der Fano Faktor ist immer $F=1$ . \end{itemize} diff --git a/statistics/code/gaussianpdf.out b/statistics/code/gaussianpdf.out new file mode 100644 index 0000000..55abadf --- /dev/null +++ b/statistics/code/gaussianpdf.out @@ -0,0 +1,5 @@ +>> gaussianpdf +The integral between 1 and 2 is 0.137 +The probability of getting a number between 1 and 2 is 0.138 +The integral between -infinity and +infinity is 1 +I.e. the probability to get any number is 1 diff --git a/statistics/lecture/statistics.tex b/statistics/lecture/statistics.tex index 678fcee..1738bc2 100644 --- a/statistics/lecture/statistics.tex +++ b/statistics/lecture/statistics.tex @@ -48,14 +48,14 @@ und die andere H\"alfte nicht kleiner als der Median ist.} \end{definition} -\begin{exercise}[mymedian.m] +\begin{exercise}{mymedian.m}{} \tr{Write a function that computes the median of a vector.} {Schreibe eine Funktion, die den Median eines Vektors zur\"uckgibt.} \end{exercise} \matlab{} stellt die Funktion \code{median()} zur Berechnung des Medians bereit. -\begin{exercise}[checkmymedian.m] +\begin{exercise}{checkmymedian.m}{} \tr{Write a script that tests whether your median function really returns a median above which are the same number of data than below. In particular the script should test data vectors of @@ -76,7 +76,7 @@ % Das mittlere Quartil entspricht dem Median. % \end{definition} -% \begin{exercise}[quartiles.m] +% \begin{exercise}{quartiles.m}{} % \tr{Write a function that computes the first, second, and third quartile of a vector.} % {Schreibe eine Funktion, die das erste, zweite und dritte Quartil als Vektor zur\"uckgibt.} % \end{exercise} @@ -89,12 +89,12 @@ Die Klassen unterteilen den Wertebereich meist in angrenzende und gleich gro{\ss}e Intervalle. Histogramme k\"onnen verwendet werden, um die Wahrscheinlichkeitsverteilung der Messwerte abzusch\"atzen. -\begin{exercise}[rollthedie.m] +\begin{exercise}{rollthedie.m}{} \tr{Write a function that simulates rolling a die $n$ times.} {Schreibe eine Funktion, die das $n$-malige W\"urfeln mit einem W\"urfel simuliert.} \end{exercise} -\begin{exercise}[diehistograms.m] +\begin{exercise}{diehistograms.m}{} \tr{Plot histograms from rolling the die 20, 100, 1000 times. Use the plain hist(x) function, force 6 bins via hist( x, 6 ), and set meaningfull bins positions.} {Plotte Histogramme von 20, 100, und @@ -128,7 +128,7 @@ des Auftretens der Gr\"o{\ss}e $x_i$ in der $i$-ten Klasse an Meistens haben wir es jedoch mit reellen Messgr\"o{\ss}en zu tun. -\begin{exercise}[gaussianbins.m] +\begin{exercise}{gaussianbins.m}{} \tr{Draw 100 random data from a Gaussian distribution and plot histograms with different bin sizes of the data.} {Ziehe 100 normalverteilte Zufallszahlen und erzeuge Histogramme mit @@ -172,7 +172,7 @@ spricht von einer Wahrscheinlichkeitsdichte. einer Wahrscheinlichkeitsdichtefunktion.} \end{figure} -\begin{exercise}[gaussianpdf.m] +\begin{exercise}{gaussianpdf.m}{gaussianpdf.out} \tr{Plot the Gaussian probability density}{Plotte die Gauss'sche Wahrscheinlichkeitsdichte } \[ p_g(x) = \frac{1}{\sqrt{2\pi\sigma^2}}e^{-\frac{(x-\mu)^2}{2\sigma^2}}\] \tr{What does it mean?}{Was bedeutet die folgende Wahrscheinlichkeit?} @@ -182,7 +182,7 @@ spricht von einer Wahrscheinlichkeitsdichte. \tr{Why?}{Warum?} \end{exercise} -\begin{exercise}[boxwhisker.m] +\begin{exercise}{boxwhisker.m}{} \tr{Generate eine $40 \times 10$ matrix of random numbers and illustrate their distribution in a box-whicker plot (\code{boxplot()} function). How to interpret the plot?}