mu = 3.0;
sigma =2.0;
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

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