Reorganized the folders and started a common script for the lectures.
This commit is contained in:
29
likelihood/code/mlemean.m
Normal file
29
likelihood/code/mlemean.m
Normal file
@@ -0,0 +1,29 @@
|
||||
% draw random numbers:
|
||||
n = 100;
|
||||
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(1, 2, 1);
|
||||
plot(pmus, lm );
|
||||
xlabel('mean')
|
||||
ylabel('likelihood')
|
||||
subplot(1, 2, 2);
|
||||
plot(pmus, loglm );
|
||||
xlabel('mean')
|
||||
ylabel('log likelihood')
|
||||
22
likelihood/lecture/Makefile
Normal file
22
likelihood/lecture/Makefile
Normal file
@@ -0,0 +1,22 @@
|
||||
BASENAME=likelihood
|
||||
PYFILES=$(wildcard *.py)
|
||||
PYPDFFILES=$(PYFILES:.py=.pdf)
|
||||
|
||||
pdf : $(BASENAME)-chapter.pdf $(PYPDFFILES)
|
||||
|
||||
$(BASENAME)-chapter.pdf : $(BASENAME)-chapter.tex $(BASENAME).tex
|
||||
pdflatex -interaction=scrollmode $< | tee /dev/stderr | fgrep -q "Rerun to get cross-references right" && pdflatex -interaction=scrollmode $< || true
|
||||
|
||||
$(PYPDFFILES) : %.pdf : %.py
|
||||
python $<
|
||||
|
||||
clean :
|
||||
rm -f *~ $(BASENAME)-chapter.aux $(BASENAME)-chapter.log $(BASENAME)-chapter.out $(BASENAME).aux $(BASENAME).log
|
||||
|
||||
cleanall : clean
|
||||
rm -f $(BASENAME)-chapter.pdf
|
||||
|
||||
watch :
|
||||
while true; do ! make -q pdf && make pdf; sleep 0.5; done
|
||||
|
||||
|
||||
225
likelihood/lecture/likelihood-chapter.tex
Normal file
225
likelihood/lecture/likelihood-chapter.tex
Normal file
@@ -0,0 +1,225 @@
|
||||
\documentclass[12pt]{report}
|
||||
|
||||
%%%%% title %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
||||
\title{\tr{Introduction to Scientific Computing}{Einf\"uhrung in die wissenschaftliche Datenverarbeitung}}
|
||||
\author{Jan Benda\\Abteilung Neuroethologie\\[2ex]\includegraphics[width=0.3\textwidth]{UT_WBMW_Rot_RGB}}
|
||||
\date{WS 15/16}
|
||||
|
||||
%%%% language %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
||||
% \newcommand{\tr}[2]{#1} % en
|
||||
% \usepackage[english]{babel}
|
||||
\newcommand{\tr}[2]{#2} % de
|
||||
\usepackage[german]{babel}
|
||||
|
||||
%%%%% packages %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
||||
\usepackage{pslatex} % nice font for pdf file
|
||||
\usepackage[breaklinks=true,bookmarks=true,bookmarksopen=true,pdfpagemode=UseNone,pdfstartview=FitH,colorlinks=true,citecolor=blue]{hyperref}
|
||||
|
||||
%%%% layout %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
||||
\usepackage[left=25mm,right=25mm,top=20mm,bottom=30mm]{geometry}
|
||||
\setcounter{tocdepth}{1}
|
||||
|
||||
%%%%% section style %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
||||
\usepackage[sf,bf,it,big,clearempty]{titlesec}
|
||||
\setcounter{secnumdepth}{1}
|
||||
|
||||
|
||||
%%%%% units %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
||||
\usepackage[mediumspace,mediumqspace,Gray]{SIunits} % \ohm, \micro
|
||||
|
||||
|
||||
%%%%% figures %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
||||
\usepackage{graphicx}
|
||||
\usepackage{xcolor}
|
||||
\pagecolor{white}
|
||||
|
||||
\newcommand{\ruler}{\par\noindent\setlength{\unitlength}{1mm}\begin{picture}(0,6)%
|
||||
\put(0,4){\line(1,0){170}}%
|
||||
\multiput(0,2)(10,0){18}{\line(0,1){4}}%
|
||||
\multiput(0,3)(1,0){170}{\line(0,1){2}}%
|
||||
\put(0,0){\makebox(0,0){{\tiny 0}}}%
|
||||
\put(10,0){\makebox(0,0){{\tiny 1}}}%
|
||||
\put(20,0){\makebox(0,0){{\tiny 2}}}%
|
||||
\put(30,0){\makebox(0,0){{\tiny 3}}}%
|
||||
\put(40,0){\makebox(0,0){{\tiny 4}}}%
|
||||
\put(50,0){\makebox(0,0){{\tiny 5}}}%
|
||||
\put(60,0){\makebox(0,0){{\tiny 6}}}%
|
||||
\put(70,0){\makebox(0,0){{\tiny 7}}}%
|
||||
\put(80,0){\makebox(0,0){{\tiny 8}}}%
|
||||
\put(90,0){\makebox(0,0){{\tiny 9}}}%
|
||||
\put(100,0){\makebox(0,0){{\tiny 10}}}%
|
||||
\put(110,0){\makebox(0,0){{\tiny 11}}}%
|
||||
\put(120,0){\makebox(0,0){{\tiny 12}}}%
|
||||
\put(130,0){\makebox(0,0){{\tiny 13}}}%
|
||||
\put(140,0){\makebox(0,0){{\tiny 14}}}%
|
||||
\put(150,0){\makebox(0,0){{\tiny 15}}}%
|
||||
\put(160,0){\makebox(0,0){{\tiny 16}}}%
|
||||
\put(170,0){\makebox(0,0){{\tiny 17}}}%
|
||||
\end{picture}\par}
|
||||
|
||||
% figures:
|
||||
\setlength{\fboxsep}{0pt}
|
||||
\newcommand{\texpicture}[1]{{\sffamily\footnotesize\input{#1.tex}}}
|
||||
%\newcommand{\texpicture}[1]{\fbox{\sffamily\footnotesize\input{#1.tex}}}
|
||||
%\newcommand{\texpicture}[1]{\setlength{\fboxsep}{2mm}\fbox{#1}}
|
||||
%\newcommand{\texpicture}[1]{}
|
||||
\newcommand{\figlabel}[1]{\textsf{\textbf{\large \uppercase{#1}}}}
|
||||
|
||||
% maximum number of floats:
|
||||
\setcounter{topnumber}{2}
|
||||
\setcounter{bottomnumber}{0}
|
||||
\setcounter{totalnumber}{2}
|
||||
|
||||
% float placement fractions:
|
||||
\renewcommand{\textfraction}{0.2}
|
||||
\renewcommand{\topfraction}{0.8}
|
||||
\renewcommand{\bottomfraction}{0.0}
|
||||
\renewcommand{\floatpagefraction}{0.5}
|
||||
|
||||
% spacing for floats:
|
||||
\setlength{\floatsep}{12pt plus 2pt minus 2pt}
|
||||
\setlength{\textfloatsep}{20pt plus 4pt minus 2pt}
|
||||
\setlength{\intextsep}{12pt plus 2pt minus 2pt}
|
||||
|
||||
% spacing for a floating page:
|
||||
\makeatletter
|
||||
\setlength{\@fptop}{0pt}
|
||||
\setlength{\@fpsep}{8pt plus 2.0fil}
|
||||
\setlength{\@fpbot}{0pt plus 1.0fil}
|
||||
\makeatother
|
||||
|
||||
% rules for floats:
|
||||
\newcommand{\topfigrule}{\vspace*{10pt}{\hrule height0.4pt}\vspace*{-10.4pt}}
|
||||
\newcommand{\bottomfigrule}{\vspace*{-10.4pt}{\hrule height0.4pt}\vspace*{10pt}}
|
||||
|
||||
% captions:
|
||||
\usepackage[format=plain,singlelinecheck=off,labelfont=bf,font={small,sf}]{caption}
|
||||
|
||||
% put caption on separate float:
|
||||
\newcommand{\breakfloat}{\end{figure}\begin{figure}[t]}
|
||||
|
||||
% references to panels of a figure within the caption:
|
||||
\newcommand{\figitem}[1]{\textsf{\bfseries\uppercase{#1}}}
|
||||
% references to figures:
|
||||
\newcommand{\panel}[1]{\textsf{\uppercase{#1}}}
|
||||
\newcommand{\fref}[1]{\textup{\ref{#1}}}
|
||||
\newcommand{\subfref}[2]{\textup{\ref{#1}}\,\panel{#2}}
|
||||
% references to figures in normal text:
|
||||
\newcommand{\fig}{Fig.}
|
||||
\newcommand{\Fig}{Figure}
|
||||
\newcommand{\figs}{Figs.}
|
||||
\newcommand{\Figs}{Figures}
|
||||
\newcommand{\figref}[1]{\fig~\fref{#1}}
|
||||
\newcommand{\Figref}[1]{\Fig~\fref{#1}}
|
||||
\newcommand{\figsref}[1]{\figs~\fref{#1}}
|
||||
\newcommand{\Figsref}[1]{\Figs~\fref{#1}}
|
||||
\newcommand{\subfigref}[2]{\fig~\subfref{#1}{#2}}
|
||||
\newcommand{\Subfigref}[2]{\Fig~\subfref{#1}{#2}}
|
||||
\newcommand{\subfigsref}[2]{\figs~\subfref{#1}{#2}}
|
||||
\newcommand{\Subfigsref}[2]{\Figs~\subfref{#1}{#2}}
|
||||
% references to figures within bracketed text:
|
||||
\newcommand{\figb}{Fig.}
|
||||
\newcommand{\figsb}{Figs.}
|
||||
\newcommand{\figrefb}[1]{\figb~\fref{#1}}
|
||||
\newcommand{\figsrefb}[1]{\figsb~\fref{#1}}
|
||||
\newcommand{\subfigrefb}[2]{\figb~\subfref{#1}{#2}}
|
||||
\newcommand{\subfigsrefb}[2]{\figsb~\subfref{#1}{#2}}
|
||||
|
||||
% references to tables:
|
||||
\newcommand{\tref}[1]{\textup{\ref{#1}}}
|
||||
% references to tables in normal text:
|
||||
\newcommand{\tab}{Tab.}
|
||||
\newcommand{\Tab}{Table}
|
||||
\newcommand{\tabs}{Tabs.}
|
||||
\newcommand{\Tabs}{Tables}
|
||||
\newcommand{\tabref}[1]{\tab~\tref{#1}}
|
||||
\newcommand{\Tabref}[1]{\Tab~\tref{#1}}
|
||||
\newcommand{\tabsref}[1]{\tabs~\tref{#1}}
|
||||
\newcommand{\Tabsref}[1]{\Tabs~\tref{#1}}
|
||||
% references to tables within bracketed text:
|
||||
\newcommand{\tabb}{Tab.}
|
||||
\newcommand{\tabsb}{Tab.}
|
||||
\newcommand{\tabrefb}[1]{\tabb~\tref{#1}}
|
||||
\newcommand{\tabsrefb}[1]{\tabsb~\tref{#1}}
|
||||
|
||||
|
||||
%%%%% equation references %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
||||
%\newcommand{\eqref}[1]{(\ref{#1})}
|
||||
\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}}
|
||||
\newcommand{\Eqnsref}[1]{\Eqns~\eqref{#1}}
|
||||
|
||||
|
||||
%%%%% listings %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
||||
\usepackage{listings}
|
||||
\lstset{
|
||||
inputpath=../code,
|
||||
basicstyle=\ttfamily\footnotesize,
|
||||
numbers=left,
|
||||
showstringspaces=false,
|
||||
language=Matlab,
|
||||
commentstyle=\itshape\color{darkgray},
|
||||
keywordstyle=\color{blue},
|
||||
stringstyle=\color{green},
|
||||
backgroundcolor=\color{blue!10},
|
||||
breaklines=true,
|
||||
breakautoindent=true,
|
||||
columns=flexible,
|
||||
frame=single,
|
||||
caption={\protect\filename@parse{\lstname}\protect\filename@base},
|
||||
captionpos=t,
|
||||
xleftmargin=1em,
|
||||
xrightmargin=1em,
|
||||
aboveskip=10pt
|
||||
}
|
||||
|
||||
%%%%% math stuff: %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
||||
\usepackage{amsmath}
|
||||
\usepackage{bm}
|
||||
\usepackage{dsfont}
|
||||
\newcommand{\naZ}{\mathds{N}}
|
||||
\newcommand{\gaZ}{\mathds{Z}}
|
||||
\newcommand{\raZ}{\mathds{Q}}
|
||||
\newcommand{\reZ}{\mathds{R}}
|
||||
\newcommand{\reZp}{\mathds{R^+}}
|
||||
\newcommand{\reZpN}{\mathds{R^+_0}}
|
||||
\newcommand{\koZ}{\mathds{C}}
|
||||
|
||||
|
||||
%%%%% structure: %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
||||
\usepackage{ifthen}
|
||||
|
||||
\newcommand{\code}[1]{\texttt{#1}}
|
||||
|
||||
\newcommand{\source}[1]{
|
||||
\begin{flushright}
|
||||
\color{gray}\scriptsize \url{#1}
|
||||
\end{flushright}
|
||||
}
|
||||
|
||||
\newenvironment{definition}[1][]{\medskip\noindent\textbf{Definition}\ifthenelse{\equal{#1}{}}{}{ #1}:\newline}%
|
||||
{\medskip}
|
||||
|
||||
\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}:}\newline \newcommand{\exercisesource}{#1}}%
|
||||
{\ifthenelse{\equal{\exercisesource}{}}{}{\ifthenelse{\value{theexercise}>\value{maxexercise}}{}{\medskip\lstinputlisting{\exercisesource}}}\medskip\stepcounter{theexercise}}
|
||||
|
||||
\graphicspath{{figures/}}
|
||||
|
||||
|
||||
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
||||
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
||||
\begin{document}
|
||||
|
||||
\include{likelihood}
|
||||
|
||||
\end{document}
|
||||
212
likelihood/lecture/likelihood.tex
Normal file
212
likelihood/lecture/likelihood.tex
Normal file
@@ -0,0 +1,212 @@
|
||||
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
||||
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
||||
\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}[mlemean.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 den Mittelwert als Parameter. Vergleiche
|
||||
die Position der Maxima mit den aus den Daten berechneten
|
||||
Mittelwerte.
|
||||
\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.
|
||||
|
||||
|
||||
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
||||
\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}
|
||||
99
likelihood/lecture/mlemean.py
Normal file
99
likelihood/lecture/mlemean.py
Normal file
@@ -0,0 +1,99 @@
|
||||
import numpy as np
|
||||
import matplotlib.pyplot as plt
|
||||
|
||||
plt.xkcd()
|
||||
fig = plt.figure( figsize=(6,5) )
|
||||
|
||||
# the data:
|
||||
n = 40
|
||||
rng = np.random.RandomState(54637281)
|
||||
sigma = 0.5
|
||||
rmu = 2.0
|
||||
xd = rng.randn(n)*sigma+rmu
|
||||
# and possible pdfs:
|
||||
x = np.arange( 0.0, 4.0, 0.01 )
|
||||
mus = [1.5, 2.0, 2.5]
|
||||
g=np.zeros((len(x), len(mus)))
|
||||
for k, mu in enumerate(mus) :
|
||||
g[:,k] = np.exp(-0.5*((x-mu)/sigma)**2.0)/np.sqrt(2.0*np.pi)/sigma
|
||||
# plot it:
|
||||
ax = fig.add_subplot( 2, 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_xlim(0.5, 3.5)
|
||||
ax.set_ylim(-0.02, 0.85)
|
||||
ax.set_xticks( np.arange(0, 5))
|
||||
ax.set_yticks( np.arange(0, 0.9, 0.2))
|
||||
ax.set_xlabel('x')
|
||||
ax.set_ylabel('Probability density')
|
||||
s = 1
|
||||
for mu in mus :
|
||||
r = 5.0*rng.rand()+2.0
|
||||
cs = 'angle3,angleA={:.0f},angleB={:.0f}'.format(90+s*r, 90-s*r)
|
||||
s *= -1
|
||||
ax.annotate('', xy=(mu, 0.02), xycoords='data',
|
||||
xytext=(mu, 0.75), textcoords='data',
|
||||
arrowprops=dict(arrowstyle="->", relpos=(0.5,0.5),
|
||||
connectionstyle=cs), zorder=1 )
|
||||
if mu > rmu :
|
||||
ax.text(mu-0.1, 0.04, '?', zorder=1, ha='right')
|
||||
else :
|
||||
ax.text(mu+0.1, 0.04, '?', zorder=1)
|
||||
for k in xrange(len(mus)) :
|
||||
ax.plot(x, g[:,k], zorder=5)
|
||||
ax.scatter(xd, 0.05*rng.rand(len(xd))+0.2, s=30, zorder=10)
|
||||
|
||||
# likelihood:
|
||||
thetas=np.arange(1.5, 2.6, 0.01)
|
||||
ps=np.zeros((len(xd),len(thetas)));
|
||||
for i, theta in enumerate(thetas) :
|
||||
ps[:,i]=np.exp(-0.5*((xd-theta)/sigma)**2.0)/np.sqrt(2.0*np.pi)/sigma
|
||||
p=np.prod(ps,axis=0)
|
||||
# plot it:
|
||||
ax = fig.add_subplot( 2, 2, 3 )
|
||||
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(r'Parameter $\theta$')
|
||||
ax.set_ylabel('Likelihood')
|
||||
ax.set_xticks( np.arange(1.6, 2.5, 0.4))
|
||||
ax.annotate('Maximum',
|
||||
xy=(2.0, 5.5e-11), xycoords='data',
|
||||
xytext=(1.0, 1.1), textcoords='axes fraction', ha='right',
|
||||
arrowprops=dict(arrowstyle="->", relpos=(0.0,0.5),
|
||||
connectionstyle="angle3,angleA=10,angleB=70") )
|
||||
ax.annotate('',
|
||||
xy=(2.0, 0), xycoords='data',
|
||||
xytext=(2.0, 5e-11), textcoords='data',
|
||||
arrowprops=dict(arrowstyle="->", relpos=(0.5,0.5),
|
||||
connectionstyle="angle3,angleA=90,angleB=80") )
|
||||
ax.plot(thetas,p)
|
||||
|
||||
ax = fig.add_subplot( 2, 2, 4 )
|
||||
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(r'Parameter $\theta$')
|
||||
ax.set_ylabel('Log-Likelihood')
|
||||
ax.set_ylim(-50,-20)
|
||||
ax.set_xticks( np.arange(1.6, 2.5, 0.4))
|
||||
ax.set_yticks( np.arange(-50, -19, 10.0))
|
||||
ax.annotate('Maximum',
|
||||
xy=(2.0, -23), xycoords='data',
|
||||
xytext=(1.0, 1.1), textcoords='axes fraction', ha='right',
|
||||
arrowprops=dict(arrowstyle="->", relpos=(0.0,0.5),
|
||||
connectionstyle="angle3,angleA=10,angleB=70") )
|
||||
ax.annotate('',
|
||||
xy=(2.0, -50), xycoords='data',
|
||||
xytext=(2.0, -26), textcoords='data',
|
||||
arrowprops=dict(arrowstyle="->", relpos=(0.5,0.5),
|
||||
connectionstyle="angle3,angleA=80,angleB=100") )
|
||||
ax.plot(thetas,np.log(p))
|
||||
|
||||
plt.tight_layout();
|
||||
plt.savefig('mlemean.pdf')
|
||||
#plt.show();
|
||||
70
likelihood/lecture/mlepdf.py
Normal file
70
likelihood/lecture/mlepdf.py
Normal file
@@ -0,0 +1,70 @@
|
||||
import numpy as np
|
||||
import scipy.stats as st
|
||||
import scipy.optimize as opt
|
||||
import matplotlib.pyplot as plt
|
||||
|
||||
plt.xkcd()
|
||||
fig = plt.figure( figsize=(6,3) )
|
||||
|
||||
# the data:
|
||||
n = 100
|
||||
shape = 2.0
|
||||
scale = 1.0
|
||||
rng = np.random.RandomState(4637281)
|
||||
xd = rng.gamma(shape, scale, n)
|
||||
|
||||
# true pdf:
|
||||
xx = np.arange(0.0, 10.1, 0.01)
|
||||
rv = st.gamma(shape)
|
||||
yy = rv.pdf(xx)
|
||||
|
||||
# mle fit:
|
||||
a = st.gamma.fit(xd, 5.0)
|
||||
yf = st.gamma.pdf(xx, *a)
|
||||
|
||||
# plot it:
|
||||
ax = fig.add_subplot( 1, 2, 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_xlim(0, 10.0)
|
||||
ax.set_ylim(0.0, 0.42)
|
||||
ax.set_xticks( np.arange(0, 11, 2))
|
||||
ax.set_yticks( np.arange(0, 0.42, 0.1))
|
||||
ax.set_xlabel('x')
|
||||
ax.set_ylabel('Probability density')
|
||||
ax.plot(xx, yy, '-', lw=5, color='#ff0000', label='pdf')
|
||||
ax.plot(xx, yf, '-', lw=2, color='#ffcc00', label='mle')
|
||||
ax.scatter(xd, 0.025*rng.rand(len(xd))+0.05, s=30, zorder=10)
|
||||
ax.legend(loc='upper right', frameon=False)
|
||||
|
||||
# histogram:
|
||||
h,b = np.histogram(xd, np.arange(0, 8.5, 1), density=True)
|
||||
|
||||
# fit histogram:
|
||||
def gammapdf(x, n, l, s) :
|
||||
return st.gamma.pdf(x, n, l, s)
|
||||
popt, pcov = opt.curve_fit(gammapdf, b[:-1]+0.5*(b[1]-b[0]), h)
|
||||
yc = st.gamma.pdf(xx, *popt)
|
||||
|
||||
# plot it:
|
||||
ax = fig.add_subplot( 1, 2, 2 )
|
||||
ax.spines['right'].set_visible(False)
|
||||
ax.spines['top'].set_visible(False)
|
||||
ax.yaxis.set_ticks_position('left')
|
||||
ax.xaxis.set_ticks_position('bottom')
|
||||
ax.set_xlim(0, 10.0)
|
||||
ax.set_xticks( np.arange(0, 11, 2))
|
||||
ax.set_xlabel('x')
|
||||
ax.set_ylim(0.0, 0.42)
|
||||
ax.set_yticks( np.arange(0, 0.42, 0.1))
|
||||
ax.set_ylabel('Probability density')
|
||||
ax.plot(xx, yy, '-', lw=5, color='#ff0000', label='pdf')
|
||||
ax.plot(xx, yc, '-', lw=2, color='#ffcc00', label='fit')
|
||||
ax.bar(b[:-1], h, np.diff(b))
|
||||
ax.legend(loc='upper right', frameon=False)
|
||||
|
||||
plt.tight_layout();
|
||||
plt.savefig('mlepdf.pdf')
|
||||
#plt.show();
|
||||
49
likelihood/lecture/mlepropline.py
Normal file
49
likelihood/lecture/mlepropline.py
Normal file
@@ -0,0 +1,49 @@
|
||||
import numpy as np
|
||||
import matplotlib.pyplot as plt
|
||||
|
||||
plt.xkcd()
|
||||
fig = plt.figure( figsize=(6,4) )
|
||||
|
||||
# the line:
|
||||
slope = 2.0
|
||||
xx = np.arange(0.0, 4.1, 0.1)
|
||||
yy = slope*xx
|
||||
# the data:
|
||||
n = 80
|
||||
rng = np.random.RandomState(218)
|
||||
sigma = 1.5
|
||||
x = 4.0*rng.rand(n)
|
||||
y = slope*x+rng.randn(n)*sigma
|
||||
# fit:
|
||||
slopef = np.sum(x*y)/np.sum(x*x)
|
||||
yf = slopef*xx
|
||||
|
||||
# plot it:
|
||||
ax = fig.add_subplot( 1, 1, 1 )
|
||||
ax.spines['left'].set_position('zero')
|
||||
ax.spines['bottom'].set_position('zero')
|
||||
ax.spines['right'].set_visible(False)
|
||||
ax.spines['top'].set_visible(False)
|
||||
ax.get_xaxis().set_tick_params(direction='inout', length=10, width=2)
|
||||
ax.get_yaxis().set_tick_params(direction='inout', length=10, width=2)
|
||||
ax.yaxis.set_ticks_position('left')
|
||||
ax.xaxis.set_ticks_position('bottom')
|
||||
ax.set_xticks(np.arange(0.0, 4.1))
|
||||
ax.set_xlim(0.0, 4.2)
|
||||
#ax.set_ylim(-1, 5)
|
||||
#ax.set_xticks( np.arange(0, 5))
|
||||
#ax.set_yticks( np.arange(0, 0.9, 0.2))
|
||||
ax.set_xlabel('x')
|
||||
ax.set_ylabel('y')
|
||||
#ax.annotate('', xy=(mu, 0.02), xycoords='data',
|
||||
# xytext=(mu, 0.75), textcoords='data',
|
||||
# arrowprops=dict(arrowstyle="->", relpos=(0.5,0.5),
|
||||
# connectionstyle=cs), zorder=1 )
|
||||
ax.scatter(x, y, label='data', s=50, zorder=10)
|
||||
ax.plot(xx, yy, 'r', lw=6.0, color='#ff0000', label='original', zorder=5)
|
||||
ax.plot(xx, yf, '--', lw=2.0, color='#ffcc00', label='fit', zorder=7)
|
||||
ax.legend(loc='upper left', frameon=False)
|
||||
|
||||
plt.tight_layout();
|
||||
plt.savefig('mlepropline.pdf')
|
||||
#plt.show();
|
||||
Reference in New Issue
Block a user