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

% 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')
savefigpdf(gcf, 'mlestd.pdf', 15, 5);