mu = 3.0;
sigma =2.0;
ns = [50, 1000];
for k = 1:length(ns)
    n = ns(k);
    fprintf('\nn=%d\n', n);
    % draw random numbers:
    x = randn(n,1)*sigma+mu;
    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.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);
        % 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
    loglm = sum(log(lms), 1);   % log likelihood
    maxlm = psigs(lm==max(lm));
    if length(maxlm) > 1
        maxlm = maxlm(1);
    end
    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);
    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, 10);