[likelihood] solved standard deviation problem

likelihood/exercises/likelihood-1.tex

@@ -37,10 +37,11 @@ of the standard deviation.
 \end{parts}
 \begin{solution}
   \lstinputlisting{mlestd.m}
+  \lstinputlisting[language={}]{mlestd.out}
   \includegraphics[width=1\textwidth]{mlestd}\\
 
   The more data the smaller the product of the probabilities ($\approx
-  p^n$ with $0 \le p < 1$) and the smaller the sum of the logarithms
+  p^n$ if $0 \le p < 1$) and the smaller the sum of the logarithms
   of the probabilities ($\approx n\log p$, note that $\log p < 0$).
 
   The product eventually gets smaller than the precision of the

likelihood/exercises/mlestd.m

@@ -6,16 +6,19 @@ for k = 1:length(ns)
     fprintf('\nn=%d\n', n);
     % 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))
+    datamean = mean(x);
+    fprintf('                                mean %.3f\n', datamean)
+    % we assume the mean to be given, therefore dof=n:
+    fprintf('                  standard deviation %.3f\n', std(x, 1))
 
     % standard deviation as parameter:
-    psigs = 1.0:0.01:3.0;
+    psigs = 1.0:0.001: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;
+        % same mean as for the standard deviation of the data:
+        p = exp(-0.5*((x-datamean)/psig).^2.0)/sqrt(2.0*pi)/psig;
         lms(:,i) = p;
     end
     lm = prod(lms, 1);              % likelihood
@@ -24,8 +27,9 @@ for k = 1:length(ns)
     if length(maxlm) > 1
         maxlm = maxlm(1);
     end
-    fprintf('      likelihood: standard deviation of the data is %.2f\n', maxlm)
-    fprintf('  log-likelihood: standard deviation of the data is %.2f\n', psigs(loglm==max(loglm)))
+    maxloglm = psigs(loglm==max(loglm));
+    fprintf('      likelihood: standard deviation %.3f\n', maxlm)
+    fprintf('  log-likelihood: standard deviation %.3f\n', maxloglm)
 
     % plot likelihood of standard deviation:
     subplot(2, 2, 2*k-1);