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;
    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
    maxlm = psigs(lm==max(lm));
    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)))

    % 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);