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

Updated maximum likelihood chapter.

Browse Source 
Common latex header file.
This commit is contained in:
Jan Benda
2015-10-31 20:45:32 +01:00
parent 91f1f3663e
commit 40a5343fa8
20 changed files with 429 additions and 1452 deletions
Download Patch File Download Diff File Expand all files Collapse all files

12
likelihood/code/mlegammafit.m Normal file
View File

@@ -0,0 +1,12 @@
% generate gamma distributed random numbers:
n = 50;
shape = 2.0;
scale = 1.0;
x = gamrnd(shape, scale, n, 1);
% maximum likelihood estimate:
p = mle(x, 'distribution', 'gamma');
% report results:
fprintf('shape=%.2f\n', p(1));
fprintf('scale=%.2f\n', p(2));

13
likelihood/code/mlemean.m
View File

@@ -3,8 +3,8 @@ 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))
fprintf(' mean of the data is %.2f\n', mean(x))
% mean as parameter:
pmus = 2.0:0.01:4.0;
@@ -18,12 +18,19 @@ end
lm = prod(lms, 1); % likelihood
loglm = sum(log(lms), 1); % log likelihood
% position of maxima:
maxlm = max(lm); % height of the maximum
pmumaxlm = pmus(lm==maxlm); % pmu where lm is at maximum
pmumaxloglm = pmus(loglm==max(loglm)); % same for loglm in one line
fprintf(' maximum of likelihood is at %.2f\n', pmumaxlm)
fprintf(' maximum of log-likelihood is at %.2f\n', pmumaxloglm)
% plot likelihood of mean:
subplot(1, 2, 1);
plot(pmus, lm );
plot(pmus, lm);
xlabel('mean')
ylabel('likelihood')
subplot(1, 2, 2);
plot(pmus, loglm );
plot(pmus, loglm);
xlabel('mean')
ylabel('log likelihood')

213
likelihood/lecture/likelihood-chapter.tex
View File

@@ -1,217 +1,8 @@
\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}
\input{../../header}
%%%% 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{red!60!black},
keywordstyle=\color{blue!50!black},
stringstyle=\color{green!50!black},
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}}
\lstset{inputpath=../code}
\graphicspath{{figures/}}

160
likelihood/lecture/likelihood.tex
View File

@@ -1,13 +1,14 @@
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\chapter{\tr{Maximum likelihood estimation}{Maximum-Likelihood-Sch\"atzer}}
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.
Maximum-Likelihood-Sch\"atzer w\"ahlen die Parameter so, dass die
Wahrscheinlichkeit, dass die Daten aus der Verteilung stammen, am
gr\"o{\ss}ten ist.
die Daten $x_1, x_2, \ldots x_n$ am besten beschreibt.
Maximum-Likelihood-Sch\"atzer (maximum likelihood estimate, mle)
w\"ahlen die Parameter so, dass die Wahrscheinlichkeit, dass die Daten
aus der Verteilung stammen, am gr\"o{\ss}ten ist.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Maximum Likelihood}
@@ -34,9 +35,14 @@ 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}
Hinter dieser Umformung steht eigentlich der Satz von Bayes. Bei der
einfachen Gleichsetzung von ${\cal L}$ mit $p$ fehlen
Normierungsfaktoren, so dass ${\cal L}$ sich nicht auf Eins
aufintegriert, und ${\cal L}$ deshalb keine
Wahrscheinlichkeit(sdichte) ist.
Wir sind nun an dem Wert des Parameters $\theta_{mle}$ interessiert, der die
Likelihood maximiert (``mle'': Maximum-Likelihood Estimate):
Likelihood maximiert (Maximum-Likelihood Estimate ``mle''):
\begin{equation}
\theta_{mle} = \text{argmax}_{\theta} {\cal L}(\theta|x_1,x_2, \ldots x_n)
\end{equation}
@@ -46,8 +52,9 @@ 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:
transformiert. Aus numerischen und 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 \\
@@ -79,19 +86,24 @@ 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}
& = & \sum_{i=1}^n - \log \sqrt{2\pi \sigma^2} -\frac{(x_i-\theta)^2}{2\sigma^2} \; .
\end{eqnarray*}
Der Logarithmus hat die sch\"one Eigenschaft die Exponentialfunktion
der Normalverteilung auszul\"oschen, da der Logarithmus die
Umkehrfunktion der Exponentialfunktion ist ($\log(e^x)=x$).
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
\Leftrightarrow \quad \theta & = & \frac{1}{n} \sum_{i=1}^n x_i \;\; = \;\; \bar x
\end{eqnarray*}
Der Maximum-Likelihood-Sch\"atzer 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.
Der Maximum-Likelihood-Sch\"atzer ist das arithmetische Mittel $\bar
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]
Ziehe $n=50$ normalverteilte Zufallsvariablen mit einem Mittelwert $\ne 0$
@@ -100,8 +112,9 @@ Normalverteilung mit diesem Mittelwert gezogen worden sind.
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
die Position der Maxima mit dem aus den Daten berechneten
Mittelwert.
\newpage
\end{exercise}
@@ -114,7 +127,7 @@ 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)
\log {\cal L}(\theta|(x_1,y_1,\sigma_1), \ldots, (x_n,y_n,\sigma_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*}
@@ -140,9 +153,12 @@ 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 eine Maximum-Likelihood Sch\"atzung. Aber nur, wenn
die Daten normalverteilt um die Funktion streuen! Bei anderen
Verteilungen m\"usste man die Log-Likelihood entsprechend
$\chi^2$ ist also eine Maximum-Likelihood Sch\"atzung.
An der Herleitung sehen wir aber auch, dass die Minimierung des
quadratischen Abstands nur dann eine Maximum-Likelihood Absch\"atzung
ist, wenn die Daten normalverteilt um die Funktion streuen. Bei
anderen Verteilungen m\"usste man die Log-Likelihood entsprechend
\eqnref{loglikelihood} ausrechnen und maximieren.
\begin{figure}[t]
@@ -168,21 +184,32 @@ und setzen diese gleich Null:
\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 im Gegensatz dazu nicht
analytisch aus den Daten berechnet werden. F\"ur diesen Fall bleibt
dann nur auf numerische Verfahren zur Optimierung der Kostenfunktion,
wie z.B. der Gradientenabstieg, zur\"uckzugreifen.
der Steigung $\theta$ des Regressionsgeraden gewonnen
(\figref{mleproplinefig}).
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. Wie z.B. die
die Steigung $m$ und der y-Achsenabschnitt $b$ einer Geradengleichung
\[ y = m \cdot x +b \]
oder allgemeiner die Koeffizienten $a_k$ eines Polynoms
\[ y = \sum_{k=0}^N a_k x^k = a_o + a_1x + a_2x^2 + a_3x^4 + \ldots \]
\matlabfun{polyfit}.
Parameter, die nichtlinear in einer Funktion enthalten sind, k\"onnen
im Gegensatz dazu nicht analytisch aus den Daten berechnet
werden. z.B. die Rate $\lambda$ eines exponentiellen Zerfalls
\[ y = c \cdot e^{\lambda x} \quad , \quad c, \lambda \in \reZ \; . \]
F\"ur diesen Fall bleibt dann nur auf numerische Verfahren zur
Optimierung der Kostenfunktion, wie z.B. der Gradientenabstieg,
zur\"uckzugreifen \matlabfun{lsqcurvefit}.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\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 wollen.
Jetzt betrachten wir noch den Fall, bei dem wir die Parameter einer
Wahrscheinlichkeitsdichtefunktion (z.B. den shape-Parameter einer
Gamma-Verteilung) an ein Datenset fitten wollen.
Ein erster Gedanke k\"onnte sein, die
Wahrscheinlichkeitsdichtefunktion durch Minimierung des quadratischen
@@ -196,16 +223,7 @@ aufintegriert. Die beiden Annahmen normalverteilte und unabh\"angige
Daten, die die Minimierung des quadratischen Abstands
\eqnref{chisqmin} zu einem Maximum-Likelihood Sch\"atzer machen, sind
also verletzt. (iii) Das Histogramm 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.
Klassenbreite ab (\figref{mlepdffig}).
\begin{figure}[t]
\includegraphics[width=1\textwidth]{mlepdf}
@@ -216,3 +234,69 @@ z.B. dem Gradientenabstieg, gel\"ost wird.
normierte Histogramm der Daten zusammen mit dem \"uber Minimierung
des quadratischen Abstands zum Histogramm berechneten Fit.}
\end{figure}
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 \matlabfun{mle}.
\begin{exercise}[mlegammafit.m]
Erzeuge Gammaverteilte Zufallszahlen und benutze Maximum-Likelihood,
um die Parameter der Gammafunktion aus den Daten zu bestimmen.
\newpage
\end{exercise}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Neuronale Kodierung}
In sensorischen Systemen kodieren Populationen von Neuronen mit ihrer
Aktivit\"at Eigenschaften von sensorischen Stimuli. z.B. im visuellen
Kortex V1 die Orientierung eines Balkens. Traditionell wird die
Antwort der Neurone f\"ur verschiedene Stimuli (z.B. verschiedene
Orientierungen des Balkens) gemessen. Die mittlere Antwort der Neurone
als Funktion eines Stimulusparameters ist dann die ``Tuning-curve''
(z.B. Feuerrate als Funktion des Orientierungswinkels).
\begin{figure}[tp]
\includegraphics[width=1\textwidth]{mlecoding}
\caption{\label{mlecodingfig} Maximum Likelihood Sch\"atzung eines
Stimulusparameters aus der Aktivit\"at einer Population von
Neuronen. Oben: Die Tuning-Kurve eines einzelnen Neurons in
Abh\"angigkeit von der Orientierung eines Balkens. Der Stimulus
der die st\"akste Aktivit\"at in diesem Neuron hervorruft ist ein
senkrechter Balken (Pfeil, $\phi_i=90$\,\degree. Die rote Fl\"ache
deutet die Variabilit\"at $p(r)$ der Aktivit\"at $r$ um die
Tuning-Kurve herum an. Mitte: Jedes Neuron in der Population hat
eine andere bevorzugte Orientierung des Stimulus (farbige Linien).
Ein Stimulus einer bestimmten Orientierung aktiviert die Neurone
in spezifischer Weise (Punkte). Unten: Die Log-Likelihood dieser
Aktivit\"aten wir maximal in der N\"ahe der wahren Orientierung
des Stimulus.}
\end{figure}
Das Gehirn ist aber mit dem umgekehrten Problem konfrontiert: gegeben
eine bestimmte Aktivit\"at der Neurone in der Population, was war der
Stimulus (die Orientierung des Balkens)? Eine m\"ogliche Antwort ist
im Sinne von Maximum-Likelihood: es war der Stimulus f\"ur den das
Aktivit\"atsmuster am wahrscheinlichsten ist.
Bleiben wir mit einem Beispiel bei den orientierungssensitiven Zellen
des V1. Das Tuning $\Omega_i(\phi)$ der Zellen $i$ auf ihre bevorzugte
Orientierung $\phi_i$ l\"asst sich gut mit einer van-Mises Funktion
(entspricht der Gaussfunktion auf einer zyklischen x-Achse)
beschreiben (\figref{mlecodingfig}):
\[ \Omega_i(\phi) = c \cdot e^{\cos(2(\phi-\phi_i))} \quad , \quad c
\in \reZ \]
Die Aktivit\"at der Neurone approximieren wir hier mit einer
Normalverteilung um die Tuning-Kurve mit Standardabweichung
$\sigma=\Omega/4$ proportional zu $\Omega$, so dass die
Wahrscheinlichkeit $p_i(r|\phi)$ des $i$-ten Neurons die Aktivit\"at $r$ zu
haben, wenn ein Stimulus mit Orientierung $\phi$ anliegt, gegeben ist durch
\[ p_i(r|\phi) = \frac{1}{\sqrt{2\pi}\Omega_i(\phi)/4} e^{-\frac{1}{2}\left(\frac{r-\Omega_i(\phi)}{\Omega_i(\phi)/4}\right)^2} \; . \]
Die Log-Likelihood der Stimulusorientierung $\phi$ gegeben die
Aktivit\"aten $r_1$, $r_2$, ... $r_n$ ist damit
\[ {\cal L}(\phi|r_1, r_2, \ldots r_n) = \sum_{i=1}^n \log p_i(r_i|\phi) \]

101
likelihood/lecture/mlecoding.py Normal file
View File

@@ -0,0 +1,101 @@
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.cm as cm
plt.xkcd()
fig = plt.figure( figsize=(6,6.8) )
rng = np.random.RandomState(4637281)
lmarg=0.1
rmarg=0.1
ax = fig.add_axes([lmarg, 0.75, 1.0-rmarg, 0.25])
ax.spines['bottom'].set_position('zero')
ax.spines['left'].set_visible(False)
ax.spines['right'].set_visible(False)
ax.spines['top'].set_visible(False)
ax.xaxis.set_ticks_position('bottom')
ax.get_yaxis().set_visible(False)
ax.set_xlim(0.0, np.pi)
ax.set_xticks(np.arange(0.125*np.pi, 1.*np.pi, 0.125*np.pi))
ax.set_xticklabels([])
ax.set_ylim(0.0, 3.5)
ax.text(-0.2, 0.5*3.5, 'Activity', rotation='vertical', va='center')
ax.annotate('Tuning curve',
xy=(0.42*np.pi, 2.5), xycoords='data',
xytext=(0.3*np.pi, 3.2), textcoords='data', ha='right',
arrowprops=dict(arrowstyle="->", relpos=(1.0,0.5),
connectionstyle="angle3,angleA=-10,angleB=110") )
ax.annotate('',
xy=(0.5*np.pi, 0.1), xycoords='data',
xytext=(0.5*np.pi, 2.6), textcoords='data',
arrowprops=dict(arrowstyle="->", relpos=(0.5,0.5),
connectionstyle="angle3,angleA=80,angleB=90") )
ax.text(0.52*np.pi, 0.7, 'preferred\norientation')
ax.plot([0, 0], [0.0, 3.5], 'k', zorder=10, clip_on=False)
xx = np.arange(0.0, 2.0*np.pi, 0.01)
pp = 0.5*np.pi
yy = np.exp(np.cos(2.0*(xx+pp)))
ax.fill_between(xx, yy+0.25*yy, yy-0.25*yy, color=cm.autumn(0.3, 1), alpha=0.5)
ax.plot(xx, yy, color=cm.autumn(0.0, 1))
ax = fig.add_axes([lmarg, 0.34, 1.0-rmarg, 0.38])
ax.spines['bottom'].set_position('zero')
ax.spines['left'].set_visible(False)
ax.spines['right'].set_visible(False)
ax.spines['top'].set_visible(False)
ax.xaxis.set_ticks_position('bottom')
ax.get_yaxis().set_visible(False)
ax.set_xlim(0.0, np.pi)
ax.set_xticks(np.arange(0.125*np.pi, 1.*np.pi, 0.125*np.pi))
ax.set_xticklabels([])
ax.set_ylim(-1.5, 3.0)
ax.text(0.5*np.pi, -1.8, 'Orientation', ha='center')
ax.text(-0.2, 0.5*3.5, 'Activity', rotation='vertical', va='center')
ax.plot([0, 0], [0.0, 3.0], 'k', zorder=10, clip_on=False)
xx = np.arange(0.0, 1.0*np.pi, 0.01)
prefphases = np.arange(0.125*np.pi, 1.*np.pi, 0.125*np.pi)
responses = []
xresponse = 0.475*np.pi
for pp in prefphases :
yy = np.exp(np.cos(2.0*(xx+pp)))
ax.plot(xx, yy, color=cm.autumn(2.0*np.abs(pp/np.pi-0.5), 1))
y = np.exp(np.cos(2.0*(xresponse+pp)))
responses.append(y + rng.randn()*0.25*y)
ax.plot(xresponse, y, '.', markersize=20, color=cm.autumn(2.0*np.abs(pp/np.pi-0.5), 1))
r=0.3
y=-0.8
ax.plot([pp-0.5*r*np.cos(pp), pp+0.5*r*np.cos(pp)], [y-r*np.sin(pp), y+r*np.sin(pp)], 'k', lw=6)
responses = np.array(responses)
ax = fig.add_axes([lmarg, 0.05, 1.0-rmarg, 0.22])
ax.spines['left'].set_visible(False)
ax.spines['right'].set_visible(False)
ax.spines['top'].set_visible(False)
ax.xaxis.set_ticks_position('bottom')
ax.get_yaxis().set_visible(False)
ax.set_xlim(0.0, np.pi)
ax.set_xticks(np.arange(0.125*np.pi, 1.*np.pi, 0.125*np.pi))
ax.set_xticklabels([])
ax.set_ylim(-1600, 0)
ax.set_xlabel('Orientiation')
ax.text(-0.2, -800, 'Log-Likelihood', rotation='vertical', va='center')
ax.plot([0, 0], [-1600, 0], 'k', zorder=10, clip_on=False)
phases = np.linspace(0.0, 1.1*np.pi, 100)
probs = np.zeros((len(responses), len(phases)))
for k, (pp, r) in enumerate(zip(prefphases, responses)) :
y = np.exp(np.cos(2.0*(phases+pp)))
sigma = 0.1*y
probs[k,:] = np.exp(-0.5*((r-y)/sigma)**2.0)/np.sqrt(2.0*np.pi)/sigma
loglikelihood = np.sum(np.log(probs), 0)
maxl = np.max(loglikelihood)
maxp = phases[np.argmax(loglikelihood)]
ax.annotate('',
xy=(maxp, -1600), xycoords='data',
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.plot(phases, loglikelihood, '-b')
plt.savefig('mlecoding.pdf')
#plt.show();