diff --git a/likelihood/exercises/exercises01.tex b/likelihood/exercises/exercises01.tex index ba01ed4..4b3dae7 100644 --- a/likelihood/exercises/exercises01.tex +++ b/likelihood/exercises/exercises01.tex @@ -15,7 +15,7 @@ \else \newcommand{\stitle}{} \fi -\header{{\bfseries\large Exercise 12\stitle}}{{\bfseries\large Maximum Likelihood}}{{\bfseries\large January 7th, 2019}} +\header{{\bfseries\large Exercise 12\stitle}}{{\bfseries\large Maximum likelihood}}{{\bfseries\large January 7th, 2019}} \firstpagefooter{Prof. Dr. Jan Benda}{Phone: 29 74573}{Email: jan.benda@uni-tuebingen.de} \runningfooter{}{\thepage}{} @@ -93,14 +93,14 @@ jan.benda@uni-tuebingen.de} Let's compute the likelihood and the log-likelihood for the estimation of the standard deviation. \begin{parts} - \part Draw $n=50$ normaly distributed random numbers with mean - $\mu=3$ and standard deviation $\sigma=2$. + \part Draw $n=50$ random numbers from a normal distribution with + mean $\mu=3$ and standard deviation $\sigma=2$. \part Plot the likelihood (computed as the product of probabilities) and the log-likelihood (sum of the logarithms of the probabilities) - using the standard deviation as the parameter we want to estimate - from the data. Compare the position of the maxima with the standard - deviation that you can compute from the data. + as a function of the standard deviation. Compare the position of the + maxima with the standard deviation that you compute directly from + the data. \part Increase $n$ to 1000. What happens to the likelihood, what happens to the log-likelihood? Why? @@ -111,75 +111,86 @@ of the standard deviation. \end{solution} %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -\question \qt{Maximum-Likelihood-Sch\"atzer einer Ursprungsgeraden} -In der Vorlesung haben wir folgende Formel f\"ur die Maximum-Likelihood -Absch\"atzung der Steigung $\theta$ einer Ursprungsgeraden durch $n$ Datenpunkte $(x_i|y_i)$ mit Standardabweichung $\sigma_i$ hergeleitet: -\[\theta = \frac{\sum_{i=1}^n \frac{x_iy_i}{\sigma_i^2}}{ \sum_{i=1}^n +\question \qt{Maximum-likelihood estimator of a line through the origin} +In the lecture we derived the following equation for an +maximum-likelihood estimate of the slope $\theta$ of a straight line +through the origin fitted to $n$ pairs of data values $(x_i|y_i)$ with +standard deviation $\sigma_i$: +\[\theta = \frac{\sum_{i=1}^n \frac{x_i y_i}{\sigma_i^2}}{ \sum_{i=1}^n \frac{x_i^2}{\sigma_i^2}} \] \begin{parts} - \part \label{mleslopefunc} Schreibe eine Funktion, die in einem $x$ und einem - $y$ Vektor die Datenpaare \"uberreicht bekommt und die Steigung der - Ursprungsgeraden, die die Likelihood maximiert, zur\"uckgibt - ($\sigma=\text{const}$). - - \part - Schreibe ein Skript, das Datenpaare erzeugt, die um eine - Ursprungsgerade mit vorgegebener Steigung streuen. Berechne mit der - Funktion aus \pref{mleslopefunc} die Steigung aus den Daten, - vergleiche mit der wahren Steigung, und plotte die urspr\"ungliche - sowie die gefittete Gerade zusammen mit den Daten. - - \part - Ver\"andere die Anzahl der Datenpunkte, die Steigung, sowie die - Streuung der Daten um die Gerade. + \part \label{mleslopefunc} Write a function that takes two vectors + $x$ and $y$ containing the data pairs and returns the slope, + computed according to this equation. For simplicity we assume + $\sigma_i=\sigma_j=\sigma$ for all $1 \le i \le n$ and $1 \le j \le + n$. How does this simplify the equation for the slope? + \begin{solution} + \lstinputlisting{mleslope.m} + \end{solution} + + \part Write a script that generates data pairs that scatter around a + line through the origin with a given slope. Use the function from + \pref{mleslopefunc} to compute the slope from the generated data. + Compare the computed slope with the true slope that has been used to + generate the data. Plot the data togehther with the line from which + the data were generated and the maximum-likelihood fit. + \begin{solution} + \lstinputlisting{mlepropfit.m} + \includegraphics[width=1\textwidth]{mlepropfit} + \end{solution} + + \part \label{mleslopecomp} Vary the number of data pairs, the slope, + as well as the variance of the data points around the true + line. Under which conditions is the maximum-likelihood estimation of + the slope closer to the true slope? + + \part To answer \pref{mleslopecomp} more precisely, generate for + each condition let's say 1000 data sets and plot a histogram of the + estimated slopes. How does the histogram, its mean and standard + deviation relate to the true slope? \end{parts} \begin{solution} - \lstinputlisting{mleslope.m} - \lstinputlisting{mlepropfit.m} - \includegraphics[width=1\textwidth]{mlepropfit} + \lstinputlisting{mlepropest.m} + \includegraphics[width=1\textwidth]{mlepropest} \end{solution} - +\continue %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -\question \qt{Maximum-Likelihood-Sch\"atzer einer Wahrscheinlichkeitsdichtefunktion} -Verschiedene Wahrscheinlichkeitsdichtefunktionen haben Parameter, die -nicht so einfach wie der Mittelwert und die Standardabweichung einer -Normalverteilung direkt aus den Daten berechnet werden k\"onnen. Solche Parameter -m\"ussen dann aus den Daten mit der Maximum-Likelihood-Methode gefittet werden. - -Um dies zu veranschaulichen ziehen wir uns diesmal nicht normalverteilte Zufallszahlen, sondern Zufallszahlen aus der Gamma-Verteilung. +\question \qt{Maximum-likelihood-estimation of a probability-density function} +Many probability-density functions have parameters that cannot be +computed directly from the data, like, for example, the mean of +normally-distributed data. Such parameter need to be estimated by +means of the maximum-likelihood from the data. + +Let us demonstrate this approach by means of data that are drawn from a +gamma distribution, \begin{parts} - \part - Finde heraus welche \code{matlab} Funktion die - Wahrscheinlichkeitsdichtefunktion (probability density function) der - Gamma-Verteilung berechnet. - - \part - 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. - - \part - Finde heraus mit welcher Funktion Gammaverteilte Zufallszahlen in - \code{matlab} gezogen werden k\"onnen. Erzeuge mit dieser Funktion - 50 Zufallszahlen mit einem der oben geplotteten ``shape'' Parameter. - - \part - Berechne und plotte ein normiertes Histogramm dieser Zufallszahlen. - - \part - Finde heraus mit welcher \code{matlab}-Funktion eine beliebige - Verteilung (``distribution'') an die Zufallszahlen nach der - Maximum-Likelihood Methode gefittet werden kann. Wie wird diese - Funktion benutzt, um die Gammaverteilung an die Daten zu fitten? - - \part - Bestimme mit dieser Funktion die Parameter der Gammaverteilung aus - den Zufallszahlen. - - \part - Plotte anschlie{\ss}end die Gammaverteilung mit den gefitteten - Parametern. + \part Find out which \code{matlab} function computes the + probability-density function of the gamma distribution. + + \part \label{gammaplot} Use this function to plot the + probability-density function of the gamma distribution for various + values of the (positive) ``shape'' parameter. Wet set the ``scale'' + parameter to one. + + \part Find out which \code{matlab} function generates random numbers + that are distributed according to a gamma distribution. Generate + with this function 50 random numbers using one of the values of the + ``shape'' parameter used in \pref{gammaplot}. + + \part Compute and plot a properly normalized histogram of these + random numbers. + + \part Find out which \code{matlab} function fit a distribution to a + vector of random numbers according to the maximum-likelihood method. + How do you need to use this function in order to fit a gamma + distribution to the data? + + \part Estimate with this function the parameter of the gamma + distribution used to generate the data. + + \part Finally, plot the fitted gamma distribution on top of the + normalized histogram of the data. \end{parts} \begin{solution} \lstinputlisting{mlepdffit.m} diff --git a/likelihood/exercises/mlepropest.m b/likelihood/exercises/mlepropest.m new file mode 100644 index 0000000..b77171c --- /dev/null +++ b/likelihood/exercises/mlepropest.m @@ -0,0 +1,25 @@ +m = 2.0; % slope +sigmas = [0.1, 1.0]; % standard deviations +ns = [100, 1000]; % number of data pairs +trials = 1000; % number of data sets + +for i = 1:length(sigmas) + sigma = sigmas(i); + for j = 1:length(ns) + n = ns(j); + slopes = zeros(trials, 1); + for k=1:trials + % data pairs: + x = 5.0*rand(n, 1); + y = m*x + sigma*randn(n, 1); + % fit: + slopes(k) = mleslope(x, y); + end + subplot(2, 2, 2*(i-1)+j); + bins = [1.9:0.005:2.1]; + hist(slopes, bins); + title(sprintf('sigma=%g, n=%d', sigma, n)); + end +end + +savefigpdf(gcf, 'mlepropest.pdf', 12, 7); diff --git a/likelihood/exercises/mlestd.m b/likelihood/exercises/mlestd.m index de37a56..149eaf7 100644 --- a/likelihood/exercises/mlestd.m +++ b/likelihood/exercises/mlestd.m @@ -1,30 +1,35 @@ -% draw random numbers: -n = 50; 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)) +ns = [50, 1000]; +for k = 1:length(ns) + n = ns(k); + % draw random numbers: + 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)) -% 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 + % 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(1, 2, 1); -plot(psigs, lm ); -xlabel('standard deviation') -ylabel('likelihood') -subplot(1, 2, 2); -plot(psigs, loglm); -xlabel('standard deviation') -ylabel('log likelihood') + % plot likelihood of standard deviation: + subplot(2, 2, 2*k-1); + plot(psigs, lm ); + title(sprintf('likelihood n=%d', n)); + xlabel('standard deviation') + ylabel('likelihood') + subplot(2, 2, 2*k); + plot(psigs, loglm); + title(sprintf('log-likelihood n=%d', n)); + xlabel('standard deviation') + ylabel('log likelihood') +end savefigpdf(gcf, 'mlestd.pdf', 15, 5);