teaching/scientificComputing
Archived
3
0
Fork 0
Code Issues Pull Requests Releases Wiki Activity

Added matlab code to mle chapter

Browse Source
This commit is contained in:
Jan Benda 2015-10-25 11:25:56 +01:00
parent 93089b4be2
commit effc38f96f
5 changed files with 462 additions and 286 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

51
statistics/code/mlemeanstd.m Normal file
View File

@ -0,0 +1,51 @@
% draw random numbers:
n = 500;
mu = 3.0;
sigma =2.0;
x = randn(n,1)*sigma+mu;
fprintf(' mean of the data is %.2f\n', mean(x))
fprintf('standard deviation of the data is %.2f\n', std(x))
% mean as parameter:
pmus = 2.0:0.01:4.0;
% matrix with the probabilities for each x and pmus:
lms = zeros(length(x), length(pmus));
for i=1:length(pmus)
pmu = pmus(i);
p = exp(-0.5*((x-pmu)/sigma).^2.0)/sqrt(2.0*pi)/sigma;
lms(:,i) = p;
end
lm = prod(lms, 1); % likelihood
loglm = sum(log(lms), 1); % log likelihood
% plot likelihood of mean:
subplot(2, 2, 1);
plot(pmus, lm );
xlabel('mean')
ylabel('likelihood')
subplot(2, 2, 2);
plot(pmus, loglm );
xlabel('mean')
ylabel('log likelihood')
% standard deviation as parameter:
psigs = 1.0:0.01:3.0;
% matrix with the probabilities for each x and psigs:
lms = zeros(length(x), length(psigs));
for i=1:length(psigs)
psig = psigs(i);
p = exp(-0.5*((x-mu)/psig).^2.0)/sqrt(2.0*pi)/psig;
lms(:,i) = p;
end
lm = prod(lms, 1); % likelihood
loglm = sum(log(lms), 1); % log likelihood
% plot likelihood of standard deviation:
subplot(2, 2, 3);
plot(psigs, lm );
xlabel('standard deviation')
ylabel('likelihood')
subplot(2, 2, 4);
plot(psigs, loglm);
xlabel('standard deviation')
ylabel('log likelihood')

27
statistics/code/mlepdffit.m Normal file
View File

@ -0,0 +1,27 @@
% plot gamma pdfs:
xx = 0.0:0.1:10.0;
shapes = [ 1.0, 2.0, 3.0, 5.0];
cc = jet(length(shapes) );
for i=1:length(shapes)
yy = gampdf(xx, shapes(i), 1.0);
plot(xx, yy, '-', 'linewidth', 3, 'color', cc(i,:), ...
'DisplayName', sprintf('s=%.0f', shapes(i)) );
hold on;
end
% generate gamma distributed random numbers:
n = 50;
x = gamrnd(3.0, 1.0, n, 1);
% histogram:
[h,b] = hist(x, 15);
h = h/sum(h)/(b(2)-b(1));
bar(b, h, 1.0, 'DisplayName', 'data');
% maximum likelihood estimate:
p = mle(x, 'distribution', 'gamma');
yy = gampdf(xx, p(1), p(2));
plot(xx, yy, '-k', 'linewidth', 5, 'DisplayName', 'mle' );
hold off;
legend('show');

29
statistics/code/mlepropfit.m Normal file
View File

@ -0,0 +1,29 @@
m = 2.0; % slope
sigma = 1.0; % standard deviation
n = 100; % number of data pairs
% data pairs:
x = 5.0*rand(n, 1);
y = m*x + sigma*randn(n, 1);
% fit:
slope = mleslope(x, y);
fprintf('slopes:\n');
fprintf('original = %.2f\n', m);
fprintf(' fit = %.2f\n', slope);
% lines:
xx = 0.0:0.1:5.0; % x-axis values
yorg = m*xx;
yfit = slope*xx;
% plot:
plot(xx, yorg, '-r', 'linewidth', 5);
hold on;
plot(xx, yfit, '-g', 'linewidth', 2);
plot(x, y, 'ob');
hold off;
legend('data', 'original', 'fit', 'Location', 'NorthWest');
legend('boxoff')
xlabel('x');
ylabel('y');

6
statistics/code/mleslope.m Normal file
View File

@ -0,0 +1,6 @@
function slope = mleslope(x, y )
% Compute the maximum likelihood estimate of the slope
% of a line through the origin
% given the data pairs in the vectors x and y.
slope = sum(x.*y)/sum(x.*x);
end

635
statistics/lecture/descriptivestatistics.tex
View File

@ -145,10 +145,10 @@
%%%%% equation references %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%\newcommand{\eqref}[1]{(\ref{#1})}
\newcommand{\eqn}{Eq.}
\newcommand{\Eqn}{Eq.}
\newcommand{\eqns}{Eqs.}
\newcommand{\Eqns}{Eqs.}
\newcommand{\eqn}{\tr{Eq}{Gl}.}
\newcommand{\Eqn}{\tr{Eq}{Gl}.}
\newcommand{\eqns}{\tr{Eqs}{Gln}.}
\newcommand{\Eqns}{\tr{Eqs}{Gln}.}
\newcommand{\eqnref}[1]{\eqn~\eqref{#1}}
\newcommand{\Eqnref}[1]{\Eqn~\eqref{#1}}
\newcommand{\eqnsref}[1]{\eqns~\eqref{#1}}
@ -205,13 +205,13 @@
\newenvironment{definition}[1][]{\medskip\noindent\textbf{Definition}\ifthenelse{\equal{#1}{}}{}{ #1}:\newline}%
{\medskip}
\newcommand{\showlisting}{yes}
%\newcommand{\showlisting}{no}
\newcounter{maxexercise}
\setcounter{maxexercise}{9} % show listings up to exercise maxexercise
\newcounter{theexercise}
\setcounter{theexercise}{1}
\newenvironment{exercise}[1][]{\medskip\noindent\textbf{\tr{Exercise}{\"Ubung}
\arabic{theexercise}:} \stepcounter{theexercise}\newline \newcommand{\exercisesource}{#1}}%
{\ifthenelse{\equal{\exercisesource}{}}{}{\ifthenelse{\equal{\showlisting}{yes}}{\medskip\lstinputlisting{\exercisesource}}{}}\medskip}
\arabic{theexercise}:}\newline \newcommand{\exercisesource}{#1}}%
{\ifthenelse{\equal{\exercisesource}{}}{}{\ifthenelse{\value{theexercise}>\value{maxexercise}}{}{\medskip\lstinputlisting{\exercisesource}}}\medskip\stepcounter{theexercise}}
\graphicspath{{figures/}}
@ -455,6 +455,347 @@ Korrelationskoeffizienten nahe 0 (\figrefb{correlationfig}).
\end{figure}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\chapter{\tr{Bootstrap Methods}{Bootstrap Methoden}}
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 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]
\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
population?}{Warum k\"onnen wir nur eine Stichprobe der
Grundgesamtheit messen?}}
\end{figure}
\begin{figure}[t]
\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
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
Bootstrap-Stichproben generiert werden (resampling) und so
Eigenschaften der Stichprobenverteilung empirisch bestimmt
werden. Aus Hesterberg et al. 2003, Bootstrap Methods and
Permuation Tests}
\end{figure}
\section{Bootstrap des Standardfehlers}
Beim Bootstrap erzeugen wir durch Resampling neue Stichproben und
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.
Plotte ein Histogramm dieser Mittelwerte, sowie deren Mittelwert und
die Standardabweichung.
Was hat das mit dem Standardfehler zu tun?
\end{exercise}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\chapter{\tr{Maximum likelihood estimation}{Maximum-Likelihood Methode}}
In vielen Situationen wollen wir einen oder mehrere Parameter $\theta$
einer Wahrscheinlichkeitsverteilung sch\"atzen, so dass die Verteilung
die Daten $x_1, x_2, \ldots x_n$ am besten beschreibt. Bei der
Maximum-Likelihood-Methode w\"ahlen wir die Parameter so, dass die
Wahrscheinlichkeit, dass die Daten aus der Verteilung stammen, am
gr\"o{\ss}ten ist.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Maximum Likelihood}
Sei $p(x|\theta)$ (lies ``Wahrscheinlichkeit(sdichte) von $x$ gegeben
$\theta$'') die Wahrscheinlichkeits(dichte)verteilung von $x$ mit dem
Parameter(n) $\theta$. Das k\"onnte die Normalverteilung
\begin{equation}
\label{normpdfmean}
p(x|\theta) = \frac{1}{\sqrt{2\pi \sigma^2}}e^{-\frac{(x-\theta)^2}{2\sigma^2}}
\end{equation}
sein mit
fester Standardverteilung $\sigma$ und dem Mittelwert $\mu$ als
Parameter $\theta$.
Wenn nun den $n$ unabh\"angigen Beobachtungen $x_1, x_2, \ldots x_n$
die Wahrscheinlichkeitsverteilung $p(x|\theta)$ zugrundeliegt, dann
ist die Verbundwahrscheinlichkeit $p(x_1,x_2, \ldots x_n|\theta)$ des
Auftretens der Werte $x_1, x_2, \ldots x_n$ gegeben ein bestimmtes $\theta$
\begin{equation}
p(x_1,x_2, \ldots x_n|\theta) = p(x_1|\theta) \cdot p(x_2|\theta)
\ldots p(x_n|\theta) = \prod_{i=1}^n p(x_i|\theta) \; .
\end{equation}
Andersherum gesehen ist das die Likelihood (deutsch immer noch ``Wahrscheinlichleit'')
den Parameter $\theta$ zu haben, gegeben die Me{\ss}werte $x_1, x_2, \ldots x_n$,
\begin{equation}
{\cal L}(\theta|x_1,x_2, \ldots x_n) = p(x_1,x_2, \ldots x_n|\theta)
\end{equation}
Wir sind nun an dem Wert des Parameters $\theta_{mle}$ interessiert, der die
Likelihood maximiert (``mle'': Maximum-Likelihood Estimate):
\begin{equation}
\theta_{mle} = \text{argmax}_{\theta} {\cal L}(\theta|x_1,x_2, \ldots x_n)
\end{equation}
$\text{argmax}_xf(x)$ bezeichnet den Wert des Arguments $x$ der Funktion $f(x)$, bei
dem $f(x)$ ihr globales Maximum annimmt. Wir suchen also den Wert von $\theta$
bei dem die Likelihood ${\cal L}(\theta)$ ihr Maximum hat.
An der Stelle eines Maximums einer Funktion \"andert sich nichts, wenn
man die Funktionswerte mit einer streng monoton steigenden Funktion
transformiert. Aus gleich ersichtlichen mathematischen Gr\"unden wird meistens
das Maximum der logarithmierten Likelihood (``Log-Likelihood'') gesucht:
\begin{eqnarray}
\theta_{mle} & = & \text{argmax}_{\theta}\; {\cal L}(\theta|x_1,x_2, \ldots x_n) \nonumber \\
& = & \text{argmax}_{\theta}\; \log {\cal L}(\theta|x_1,x_2, \ldots x_n) \nonumber \\
& = & \text{argmax}_{\theta}\; \log \prod_{i=1}^n p(x_i|\theta) \nonumber \\
& = & \text{argmax}_{\theta}\; \sum_{i=1}^n \log p(x_i|\theta) \label{loglikelihood}
\end{eqnarray}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Beispiel: Das arithmetische Mittel}
Wenn die Me{\ss}daten $x_1, x_2, \ldots x_n$ der Normalverteilung \eqnref{normpdfmean}
entstammen, und wir den Mittelwert $\mu$ als einzigen Parameter der Verteilung betrachten,
welcher Wert von $\theta$ maximiert dessen Likelhood?
\begin{figure}[t]
\includegraphics[width=1\textwidth]{mlemean}
\caption{\label{mlemeanfig} Maximum Likelihood Estimation des
Mittelwerts. Oben: Die Daten zusammen mit drei m\"oglichen
Normalverteilungen mit unterschiedlichen Mittelwerten (Pfeile) aus
denen die Daten stammen k\"onnten. Unteln links: Die Likelihood
in Abh\"angigkeit des Mittelwerts als Parameter der
Normalverteilungen. Unten rechts: die entsprechende
Log-Likelihood. An der Position des Maximums bei $\theta=2$
\"andert sich nichts (Pfeil).}
\end{figure}
Die Log-Likelihood \eqnref{loglikelihood} ist
\begin{eqnarray*}
\log {\cal L}(\theta|x_1,x_2, \ldots x_n)
& = & \sum_{i=1}^n \log \frac{1}{\sqrt{2\pi \sigma^2}}e^{-\frac{(x_i-\theta)^2}{2\sigma^2}} \\
& = & \sum_{i=1}^n - \log \sqrt{2\pi \sigma^2} -\frac{(x_i-\theta)^2}{2\sigma^2}
\end{eqnarray*}
Zur Bestimmung des Maximums der Log-Likelihood berechnen wir deren Ableitung
nach dem Parameter $\theta$ und setzen diese gleich Null:
\begin{eqnarray*}
\frac{\text{d}}{\text{d}\theta} \log {\cal L}(\theta|x_1,x_2, \ldots x_n) & = & \sum_{i=1}^n \frac{2(x_i-\theta)}{2\sigma^2} \;\; = \;\; 0 \\
\Leftrightarrow \quad \sum_{i=1}^n x_i - \sum_{i=1}^n x_i \theta & = & 0 \\
\Leftrightarrow \quad n \theta & = & \sum_{i=1}^n x_i \\
\Leftrightarrow \quad \theta & = & \frac{1}{n} \sum_{i=1}^n x_i
\end{eqnarray*}
Der Maximum-Likelihood-Estimator ist das arithmetische Mittel der Daten. D.h.
das arithmetische Mittel maximiert die Wahrscheinlichkeit, dass die Daten aus einer
Normalverteilung mit diesem Mittelwert gezogen worden sind.
\begin{exercise}[mlemeanstd.m]
Ziehe $n=50$ normalverteilte Zufallsvariablen mit einem Mittelwert $\ne 0$
und einer Standardabweichung $\ne 1$.
Plotte die Likelihood (aus dem Produkt der Wahrscheinlichkeiten) und
die Log-Likelihood (aus der Summe der logarithmierten
Wahrscheinlichkeiten) f\"ur (1) den Mittelwert und (2) die
Standardabweichung. Vergleiche die Position der Maxima mit den
aus den Daten berechneten Mittelwerten und Standardabweichungen.
Erh\"ohe $n$ auf 1000. Was passiert mit der Likelihood, was mit der Log-Likelihood?
\end{exercise}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Kurvenfit als Maximum Likelihood Estimation}
Beim Kurvenfit soll eine Funktion $f(x;\theta)$ mit den Parametern
$\theta$ an die Datenpaare $(x_i|y_i)$ durch Anpassung der Parameter
$\theta$ gefittet werden. Wenn wir annehmen, dass die $y_i$ um die
entsprechenden Funktionswerte $f(x_i;\theta)$ mit einer
Standardabweichung $\sigma_i$ normalverteilt streuen, dann lautet die
Log-Likelihood
\begin{eqnarray*}
\log {\cal L}(\theta|x_1,x_2, \ldots x_n)
& = & \sum_{i=1}^n \log \frac{1}{\sqrt{2\pi \sigma_i^2}}e^{-\frac{(y_i-f(x_i;\theta))^2}{2\sigma_i^2}} \\
& = & \sum_{i=1}^n - \log \sqrt{2\pi \sigma_i^2} -\frac{(x_i-f(y_i;\theta))^2}{2\sigma_i^2} \\
\end{eqnarray*}
Der einzige Unterschied zum vorherigen Beispiel ist, dass die
Mittelwerte der Normalverteilungen nun durch die Funktionswerte
gegeben sind.
Der Parameter $\theta$ soll so gew\"ahlt werden, dass die
Log-Likelihood maximal wird. Der erste Term der Summe ist
unabh\"angig von $\theta$ und kann deshalb bei der Suche nach dem
Maximum weggelassen werden.
\begin{eqnarray*}
& = & - \frac{1}{2} \sum_{i=1}^n \left( \frac{y_i-f(x_i;\theta)}{\sigma_i} \right)^2
\end{eqnarray*}
Anstatt nach dem Maximum zu suchen, k\"onnen wir auch das Vorzeichen der Log-Likelihood
umdrehen und nach dem Minimum suchen. Dabei k\"onnen wir auch den Faktor $1/2$ vor der Summe vernachl\"assigen --- auch das \"andert nichts an der Position des Minimums.
\begin{equation}
\theta_{mle} = \text{argmin}_{\theta} \; \sum_{i=1}^n \left( \frac{y_i-f(x_i;\theta)}{\sigma_i} \right)^2 \;\; = \;\; \text{argmin}_{\theta} \; \chi^2
\end{equation}
Die Summer der quadratischen Abst\"ande normiert auf die jeweiligen
Standardabweichungen wird auch mit $\chi^2$ bezeichnet. Der Wert des
Parameters $\theta$ welcher den quadratischen Abstand minimiert ist
also identisch mit der Maximierung der Wahrscheinlichkeit, dass die
Daten tats\"achlich aus der Funktion stammen k\"onnen. Minimierung des
$\chi^2$ ist also ein Maximum-Likelihood Estimate.
\begin{figure}[t]
\includegraphics[width=1\textwidth]{mlepropline}
\caption{\label{mleproplinefig} Maximum Likelihood Estimation der
Steigung einer Ursprungsgeraden.}
\end{figure}
\subsection{Beispiel: einfache Proportionalit\"at}
Als Funktion nehmen wir die Ursprungsgerade
\[ f(x) = \theta x \]
mit Steigung $\theta$. Die $\chi^2$-Summe lautet damit
\[ \chi^2 = \sum_{i=1}^n \left( \frac{y_i-\theta x_i}{\sigma_i} \right)^2 \; . \]
Zur Bestimmung des Minimums berechnen wir wieder die erste Ableitung nach $\theta$
und setzen diese gleich Null:
\begin{eqnarray}
\frac{\text{d}}{\text{d}\theta}\chi^2 & = & \frac{\text{d}}{\text{d}\theta} \sum_{i=1}^n \left( \frac{y_i-\theta x_i}{\sigma_i} \right)^2 \nonumber \\
& = & \sum_{i=1}^n \frac{\text{d}}{\text{d}\theta} \left( \frac{y_i-\theta x_i}{\sigma_i} \right)^2 \nonumber \\
& = & -2 \sum_{i=1}^n \frac{x_i}{\sigma_i} \left( \frac{y_i-\theta x_i}{\sigma_i} \right) \nonumber \\
& = & -2 \sum_{i=1}^n \left( \frac{x_iy_i}{\sigma_i^2} - \theta \frac{x_i^2}{\sigma_i^2} \right) \;\; = \;\; 0 \nonumber \\
\Leftrightarrow \quad \theta \sum_{i=1}^n \frac{x_i^2}{\sigma_i^2} & = & \sum_{i=1}^n \frac{x_iy_i}{\sigma_i^2} \nonumber \\
\Leftrightarrow \quad \theta & = & \frac{\sum_{i=1}^n \frac{x_iy_i}{\sigma_i^2}}{ \sum_{i=1}^n \frac{x_i^2}{\sigma_i^2}} \label{mleslope}
\end{eqnarray}
Damit haben wir nun einen anlytischen Ausdruck f\"ur die Bestimmung
der Steigung $\theta$ des Regressionsgeraden gewonnen. Ein
Gradientenabstieg ist f\"ur das Fitten der Geradensteigung also gar nicht
n\"otig. Das gilt allgemein f\"ur das Fitten von Koeffizienten von
linear kombinierten Basisfunktionen. Parameter die nichtlinear in
einer Funktion enthalten sind k\"onnen aber nicht analytisch aus den
Daten berechnet werden. Da bleibt dann nur auf numerische Verfahren
zur Optimierung der Kostenfunktion, wie z.B. der Gradientenabstieg,
zur\"uckzugreifen.
\begin{exercise}[mleslope.m]
Schreibe eine Funktion, die in einem $x$ und einem $y$ Vektor die
Datenpaare \"uberreicht bekommt und die Steigung der
Ursprungsgeraden \eqnref{mleslope}, die die Likelihood maximiert,
zur\"uckgibt ($\sigma=1$).
\end{exercise}
\begin{exercise}[mlepropfit.m]
Schreibe ein Skript, das Datenpaare erzeugt, die um eine
Ursprungsgerade mit vorgegebener Steigung streuen. Berechne mit der
Funktion die Steigung aus den Daten, vergleiche mit der wahren
Steigung, und plotte die urspr\"ungliche sowie die gefittete Gerade
zusammen mit den Daten.
Ver\"andere die Anzahl der Datenpunkte, die Steigung, sowie die
Streuung der Daten um die Gerade.
\end{exercise}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Fits von Wahrscheinlichkeitsverteilungen}
Zum Abschluss betrachten wir noch den Fall, bei dem wir die Parameter
einer Wahrscheinlichkeitsdichtefunktion (z.B. Mittelwert und
Standardabweichung der Normalverteilung) an ein Datenset fitten wolle.
Ein erster Gedanke k\"onnte sein, die
Wahrscheinlichkeitsdichtefunktion durch Minimierung des quadratischen
Abstands an ein Histogram der Daten zu fitten. Das ist aber aus
folgenden Gr\"unden nicht die Methode der Wahl: (i)
Wahrscheinlichkeitsdichten k\"onnen nur positiv sein. Darum k\"onnen
insbesondere bei kleinen Werten die Daten nicht symmetrisch streuen,
wie es normalverteilte Daten machen sollten. (ii) Die Datenwerte sind
nicht unabh\"angig, da das normierte Histogram sich zu Eins
aufintegriert. Die beiden Annahmen normalverteilte und unabh\"angige Daten
die die Minimierung des quadratischen Abstands zu einem Maximum
Likelihood Estimator machen sind also verletzt. (iii) Das Histgramm
h\"angt von der Wahl der Klassenbreite ab.
Den direkten Weg, eine Wahrscheinlichkeitsdichtefunktion an ein
Datenset zu fitten, haben wir oben schon bei dem Beispiel zur
Absch\"atzung des Mittelwertes einer Normalverteilung gesehen ---
Maximum Likelihood! Wir suchen einfach die Parameter $\theta$ der
gesuchten Wahrscheinlichkeitsdichtefunktion bei der die Log-Likelihood
\eqnref{loglikelihood} maximal wird. Das ist im allgemeinen ein
nichtlinieares Optimierungsproblem, das mit numerischen Verfahren, wie
z.B. dem Gradientenabstieg, gel\"ost wird.
\begin{figure}[t]
\includegraphics[width=1\textwidth]{mlepdf}
\caption{\label{mlepdffig} Maximum Likelihood Estimation einer
Wahrscheinlichkeitsdichtefunktion. Links: die 100 Datenpunkte, die aus der Gammaverteilung
2. Ordnung (rot) gezogen worden sind. Der Maximum-Likelihood-Fit ist orange dargestellt.
Rechts: das normierte Histogramm der Daten zusammen mit der \"uber Minimierung
des quadratischen Abstands zum Histogramm berechneten Fits ist potentiell schlechter.}
\end{figure}
\begin{exercise}[mlepdffit.m]
Zur Abwechslung ziehen wir uns diesmal Zufallszahlen, die nicht
einer Normalverteilung entstammen, sonder aus der Gamma-Verteilung.
Finde heraus welche Funktion die Wahrscheinlichkeitsdichtefunktion
(probability density function) der Gamma-Verteilung in \code{matlab}
berechnet.
Plotte mit Hilfe dieser Funktion die Wahrscheinlichkeitsdichtefunktion
der Gamma-Verteilung f\"ur verschiedene Werte des (positiven) ``shape'' Parameters.
Den ``scale'' Parameter setzen wir auf Eins.
Finde heraus mit welcher Funktion Gamma-verteilte Zufallszahlen in
\code{matlab} gezogen werden k\"onnen. Erzeuge mit dieser Funktion
50 Zufallszahlen mit einem der oben geplotteten ``shape'' Parameter.
Berechne und plotte ein normiertes Histogramm dieser Zufallszahlen.
Finde heraus mit welcher \code{matlab}-Funktion die Gammaverteilung
an die Zufallszahlen nach der Maximum-Likelihood Methode gefittet
werden kann. Bestimme mit dieser Funktion die Parameter der
Gammaverteilung aus den Zufallszahlen. Plotte anschlie{\ss}end
die Gammaverteilung mit den gefitteten Parametern.
\end{exercise}
\end{document}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Statistics}
What is "a statistic"? % dt. Sch\"atzfunktion
\begin{definition}[statistic]
A statistic (singular) is a single measure of some attribute of a
sample (e.g., its arithmetic mean value). It is calculated by
applying a function (statistical algorithm) to the values of the
items of the sample, which are known together as a set of data.
\source{http://en.wikipedia.org/wiki/Statistic}
\end{definition}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Data types}
@ -574,281 +915,3 @@ Korrelationskoeffizienten nahe 0 (\figrefb{correlationfig}).
\end{itemize}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\chapter{\tr{Bootstrap Methods}{Bootstrap Methoden}}
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 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]
\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
population?}{Warum k\"onnen wir nur eine Stichprobe der
Grundgesamtheit messen?}}
\end{figure}
\begin{figure}[t]
\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
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
Bootstrap-Stichproben generiert werden (resampling) und so
Eigenschaften der Stichprobenverteilung empirisch bestimmt
werden. Aus Hesterberg et al. 2003, Bootstrap Methods and
Permuation Tests}
\end{figure}
\section{Bootstrap des Standardfehlers}
Beim Bootstrap erzeugen wir durch Resampling neue Stichproben und
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.
Plotte ein Histogramm dieser Mittelwerte, sowie deren Mittelwert und
die Standardabweichung.
Was hat das mit dem Standardfehler zu tun?
\end{exercise}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\chapter{\tr{Maximum likelihood estimation}{Maximum-Likelihood Methode}}
In vielen Situationen wollen wir einen oder mehrere Parameter $\theta$
einer Wahrscheinlichkeitsverteilung sch\"atzen, so dass die Verteilung
die Daten $x_1, x_2, \ldots x_n$ am besten beschreibt. Bei der
Maximum-Likelihood-Methode w\"ahlen wir die Parameter so, dass die
Wahrscheinlichkeit, dass die Daten aus der Verteilung stammen, am
gr\"o{\ss}ten ist.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Maximum Likelihood}
Sei $p(x|\theta)$ (lies ``Wahrscheinlichkeit(sdichte) von $x$ gegeben
$\theta$'') die Wahrscheinlichkeits(dichte)verteilung von $x$ mit dem
Parameter(n) $\theta$. Das k\"onnte die Normalverteilung
\begin{equation}
\label{normpdfmean}
p(x|\theta) = \frac{1}{\sqrt{2\pi \sigma^2}}e^{-\frac{(x-\theta)^2}{2\sigma^2}}
\end{equation}
sein mit
fester Standardverteilung $\sigma$ und dem Mittelwert $\mu$ als
Parameter $\theta$.
Wenn nun den $n$ unabh\"angigen Beobachtungen $x_1, x_2, \ldots x_n$
die Wahrscheinlichkeitsverteilung $p(x|\theta)$ zugrundeliegt, dann
ist die Verbundwahrscheinlichkeit $p(x_1,x_2, \ldots x_n|\theta)$ des
Auftretens der Werte $x_1, x_2, \ldots x_n$ gegeben ein bestimmtes $\theta$
\[ p(x_1,x_2, \ldots x_n|\theta) = p(x_1|\theta) \cdot p(x_2|\theta)
\ldots p(x_n|\theta) = \prod_{i=1}^n p(x_i|\theta) \; .\]
Andersherum gesehen ist das die Likelihood (deutsch immer noch ``Wahrscheinlichleit'')
den Parameter $\theta$ zu haben, gegeben die Me{\ss}werte $x_1, x_2, \ldots x_n$,
\[ {\cal L}(\theta|x_1,x_2, \ldots x_n) = p(x_1,x_2, \ldots x_n|\theta) \]
Wir sind nun an dem Wert des Parameters $\theta_{mle}$ interessiert, der die
Likelihood maximiert (``mle'': Maximum-Likelihood Estimate):
\[ \theta_{mle} = \text{argmax}_{\theta} {\cal L}(\theta|x_1,x_2,
\ldots x_n) \]
$\text{argmax}_xf(x)$ bezeichnet den Wert des Arguments $x$ der Funktion $f(x)$, bei
dem $f(x)$ ihr globales Maximum annimmt. Wir suchen also den Wert von $\theta$
bei dem die Likelihood ${\cal L}(\theta)$ ihr Maximum hat.
An der Stelle eines Maximums einer Funktion \"andert sich nichts, wenn
man die Funktionswerte mit einer streng monoton steigenden Funktion
transformiert. Aus gleich ersichtlichen mathematischen Gr\"unden wird meistens
das Maximum der logarithmierten Likelihood (``Log-Likelihood'') gesucht:
\begin{eqnarray}
\theta_{mle} & = & \text{argmax}_{\theta}\; {\cal L}(\theta|x_1,x_2, \ldots x_n) \nonumber \\
& = & \text{argmax}_{\theta}\; \log {\cal L}(\theta|x_1,x_2, \ldots x_n) \nonumber \\
& = & \text{argmax}_{\theta}\; \log \prod_{i=1}^n p(x_i|\theta) \nonumber \\
& = & \text{argmax}_{\theta}\; \sum_{i=1}^n \log p(x_i|\theta) \label{loglikelihood}
\end{eqnarray}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Beispiel: Das arithmetische Mittel}
Wenn die Me{\ss}daten $x_1, x_2, \ldots x_n$ der Normalverteilung \eqnref{normpdfmean}
entstammen, und wir den Mittelwert $\mu$ als einzigen Parameter der Verteilung betrachten,
welcher Wert von $\theta$ maximiert dessen Likelhood?
\begin{figure}[t]
\includegraphics[width=1\textwidth]{mlemean}
\caption{\label{mlemeanfig} Maximum Likelihood Estimation des
Mittelwerts. Oben: Die Daten zusammen mit drei m\"oglichen
Normalverteilungen mit unterschiedlichen Mittelwerten (Pfeile) aus
denen die Daten stammen k\"onnten. Unteln links: Die Likelihood
in Abh\"angigkeit des Mittelwerts als Parameter der
Normalverteilungen. Unten rechts: die entsprechende
Log-Likelihood. An der Position des Maximums bei $\theta=2$
\"andert sich nichts (Pfeil).}
\end{figure}
Die Log-Likelihood \eqnref{loglikelihood} ist
\begin{eqnarray*}
\log {\cal L}(\theta|x_1,x_2, \ldots x_n)
& = & \sum_{i=1}^n \log \frac{1}{\sqrt{2\pi \sigma^2}}e^{-\frac{(x_i-\theta)^2}{2\sigma^2}} \\
& = & \sum_{i=1}^n - \log \sqrt{2\pi \sigma^2} -\frac{(x_i-\theta)^2}{2\sigma^2}
\end{eqnarray*}
Zur Bestimmung des Maximums der Log-Likelihood berechnen wir deren Ableitung
nach dem Parameter $\theta$ und setzen diese gleich Null:
\begin{eqnarray*}
\frac{\text{d}}{\text{d}\theta} \log {\cal L}(\theta|x_1,x_2, \ldots x_n) & = & \sum_{i=1}^n \frac{2(x_i-\theta)}{2\sigma^2} \;\; = \;\; 0 \\
\Leftrightarrow \quad \sum_{i=1}^n x_i - \sum_{i=1}^n x_i \theta & = & 0 \\
\Leftrightarrow \quad n \theta & = & \sum_{i=1}^n x_i \\
\Leftrightarrow \quad \theta & = & \frac{1}{n} \sum_{i=1}^n x_i
\end{eqnarray*}
Der Maximum-Likelihood-Estimator ist das arithmetische Mittel der Daten. D.h.
das arithmetische Mittel maximiert die Wahrscheinlichkeit, dass die Daten aus einer
Normalverteilung mit diesem Mittelwert gezogen worden sind.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Kurvenfit als Maximum Likelihood Estimation}
Beim Kurvenfit soll eine Funktion $f(x;\theta)$ mit den Parametern
$\theta$ an die Datenpaare $(x_i|y_i)$ durch Anpassung der Parameter
$\theta$ gefittet werden. Wenn wir annehmen, dass die $y_i$ um die
entsprechenden Funktionswerte $f(x_i;\theta)$ mit einer
Standardabweichung $\sigma_i$ normalverteilt streuen, dann lautet die
Log-Likelihood
\begin{eqnarray*}
\log {\cal L}(\theta|x_1,x_2, \ldots x_n)
& = & \sum_{i=1}^n \log \frac{1}{\sqrt{2\pi \sigma_i^2}}e^{-\frac{(y_i-f(x_i;\theta))^2}{2\sigma_i^2}} \\
& = & \sum_{i=1}^n - \log \sqrt{2\pi \sigma_i^2} -\frac{(x_i-f(y_i;\theta))^2}{2\sigma_i^2} \\
\end{eqnarray*}
Der einzige Unterschied zum vorherigen Beispiel ist, dass die
Mittelwerte der Normalverteilungen nun durch die Funktionswerte
gegeben sind.
Der Parameter $\theta$ soll so gew\"ahlt werden, dass die
Log-Likelihood maximal wird. Der erste Term der Summe ist
unabh\"angig von $\theta$ und kann deshalb bei der Suche nach dem
Maximum weggelassen werden.
\begin{eqnarray*}
& = & - \frac{1}{2} \sum_{i=1}^n \left( \frac{y_i-f(x_i;\theta)}{\sigma_i} \right)^2
\end{eqnarray*}
Anstatt nach dem Maximum zu suchen, k\"onnen wir auch das Vorzeichen der Log-Likelihood
umdrehen und nach dem Minimum suchen. Dabei k\"onnen wir auch den Faktor $1/2$ vor der Summe vernachl\"assigen --- auch das \"andert nichts an der Position des Minimums.
\begin{eqnarray*}
\theta_{mle} & = & \text{argmin}_{\theta} \; \sum_{i=1}^n \left( \frac{y_i-f(x_i;\theta)}{\sigma_i} \right)^2 \;\; = \;\; \text{argmin}_{\theta} \; \chi^2
\end{eqnarray*}
Die Summer der quadratischen Abst\"ande normiert auf die jeweiligen
Standardabweichungen wird auch mit $\chi^2$ bezeichnet. Der Wert des
Parameters $\theta$ welcher den quadratischen Abstand minimiert ist
also identisch mit der Maximierung der Wahrscheinlichkeit, dass die
Daten tats\"achlich aus der Funktion stammen k\"onnen. Minimierung des
$\chi^2$ ist also ein Maximum-Likelihood Estimate.
\begin{figure}[t]
\includegraphics[width=1\textwidth]{mlepropline}
\caption{\label{mleproplinefig} Maximum Likelihood Estimation der
Steigung einer Ursprungsgeraden.}
\end{figure}
\subsection{Beispiel: einfache Proportionalit\"at}
Als Funktion nehmen wir die Ursprungsgerade
\[ f(x) = \theta x \]
mit Steigung $\theta$. Die $\chi^2$-Summe lautet damit
\[ \chi^2 = \sum_{i=1}^n \left( \frac{y_i-\theta x_i}{\sigma_i} \right)^2 \; . \]
Zur Bestimmung des Minimums berechnen wir wieder die erste Ableitung nach $\theta$
und setzen diese gleich Null:
\begin{eqnarray*}
\frac{\text{d}}{\text{d}\theta}\chi^2 & = & \frac{\text{d}}{\text{d}\theta} \sum_{i=1}^n \left( \frac{y_i-\theta x_i}{\sigma_i} \right)^2 \\
& = & \sum_{i=1}^n \frac{\text{d}}{\text{d}\theta} \left( \frac{y_i-\theta x_i}{\sigma_i} \right)^2 \\
& = & -2 \sum_{i=1}^n \frac{x_i}{\sigma_i} \left( \frac{y_i-\theta x_i}{\sigma_i} \right) \\
& = & -2 \sum_{i=1}^n \left( \frac{x_iy_i}{\sigma_i^2} - \theta \frac{x_i^2}{\sigma_i^2} \right) \;\; = \;\; 0 \\
\Leftrightarrow \quad \theta \sum_{i=1}^n \frac{x_i^2}{\sigma_i^2} & = & \sum_{i=1}^n \frac{x_iy_i}{\sigma_i^2} \\
\Leftrightarrow \quad \theta & = & \frac{\sum_{i=1}^n \frac{x_iy_i}{\sigma_i^2}}{ \sum_{i=1}^n \frac{x_i^2}{\sigma_i^2}}
\end{eqnarray*}
Damit haben wir nun einen anlytischen Ausdruck f\"ur die Bestimmung
der Steigung $\theta$ des Regressionsgeraden gewonnen. Ein
Gradientenabstieg ist f\"ur das Fitten der Geradensteigung also gar nicht
n\"otig. Das gilt allgemein f\"ur das fitten von Koeffizienten von
linear kombinierten Basisfunktionen. Parameter die nichtlinear in
einer Funktion enthalten sind k\"onnen aber nicht analytisch aus den
Daten berechnet werden. Da bleibt dann nur auf numerische Verfahren
zur Optimierung der Kostenfunktion, wie z.B. der Gradientenabstieg,
zur\"uckzugreifen.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Fits von Wahrscheinlichkeitsverteilungen}
Zum Abschluss betrachten wir noch den Fall, bei dem wir die Parameter
einer Wahrscheinlichkeitsdichtefunktion (z.B. Mittelwert und
Standardabweichung der Normalverteilung) an ein Datenset fitten wolle.
Ein erster Gedanke k\"onnte sein, die
Wahrscheinlichkeitsdichtefunktion durch Minimierung des quadratischen
Abstands an ein Histogram der Daten zu fitten. Das ist aber aus
folgenden Gr\"unden nicht die Methode der Wahl: (i)
Wahrscheinlichkeitsdichten k\"onnen nur positiv sein. Darum k\"onnen
insbesondere bei kleinen Werten die Daten nicht symmetrisch streuen,
wie es normalverteilte Daten machen sollten. (ii) Die Datenwerte sind
nicht unabh\"angig, da das normierte Histogram sich zu Eins
aufintegriert. Die beiden Annahmen normalverteilte und unabh\"angige Daten
die die Minimierung des quadratischen Abstands zu einem Maximum
Likelihood Estimator machen sind also verletzt.
Den direkten Weg, eine Wahrscheinlichkeitsdichtefunktion an ein
Datenset zu fitten, haben wir oben schon bei dem Beispiel zur
Absch\"atzung des Mittelwertes einer Normalverteilung gesehen ---
Maximum Likelihood! Wir suchen einfach die Parameter $\theta$ der
gesuchten Wahrscheinlichkeitsdichtefunktion bei der die Log-Likelihood
\eqnref{loglikelihood} maximal wird. Das ist im allgemeinen ein
nichtlinieares Optimierungsproblem, das mit numerischen Verfahren, wie
z.B. dem Gradientenabstieg, gel\"ost wird.
\begin{figure}[t]
\includegraphics[width=1\textwidth]{mlepdf}
\caption{\label{mlepdffig} Maximum Likelihood Estimation einer
Wahrscheinlichkeitsdichtefunktion. Links: die 100 Datenpunkte, die aus der Gammaverteilung
2. Ordnung (rot) gezogen worden sind. Der Maximum-Likelihood-Fit ist orange dargestellt.
Rechts: das normierte Histogramm der Daten zusammen mit der \"uber Minimierung
des quadratischen Abstands zum Histogramm berechneten Fits.}
\end{figure}
\end{document}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Statistics}
What is "a statistic"? % dt. Sch\"atzfunktion
\begin{definition}[statistic]
A statistic (singular) is a single measure of some attribute of a
sample (e.g., its arithmetic mean value). It is calculated by
applying a function (statistical algorithm) to the values of the
items of the sample, which are known together as a set of data.
\source{http://en.wikipedia.org/wiki/Statistic}
\end{definition}