Merge branch 'master' of raven.am28.uni-tuebingen.de:scientificComputing
Conflicts: header.tex
This commit is contained in:
commit
328e06b2c0
2
Makefile
2
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 :
|
||||
|
||||
5
bootstrap/code/bootstrapsem.out
Normal file
5
bootstrap/code/bootstrapsem.out
Normal file
@ -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 -
|
||||
38
bootstrap/code/correlationsignificance.m
Normal file
38
bootstrap/code/correlationsignificance.m
Normal file
@ -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;
|
||||
4
bootstrap/code/correlationsignificance.out
Normal file
4
bootstrap/code/correlationsignificance.out
Normal file
@ -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
|
||||
@ -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}
|
||||
|
||||
|
||||
53
bootstrap/lecture/bootstrapsem.py
Normal file
53
bootstrap/lecture/bootstrapsem.py
Normal file
@ -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();
|
||||
60
bootstrap/lecture/permutecorrelation.py
Normal file
60
bootstrap/lecture/permutecorrelation.py
Normal file
@ -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();
|
||||
33
header.tex
33
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}}
|
||||
|
||||
3
likelihood/code/mlegammafit.out
Normal file
3
likelihood/code/mlegammafit.out
Normal file
@ -0,0 +1,3 @@
|
||||
>> mlegammafit
|
||||
shape=2.10
|
||||
scale=0.95
|
||||
5
likelihood/code/mlemean.out
Normal file
5
likelihood/code/mlemean.out
Normal file
@ -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
|
||||
@ -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
|
||||
|
||||
@ -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')
|
||||
|
||||
@ -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}
|
||||
|
||||
|
||||
5
statistics/code/gaussianpdf.out
Normal file
5
statistics/code/gaussianpdf.out
Normal file
@ -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
|
||||
@ -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?}
|
||||
|
||||
Reference in New Issue
Block a user