% draw random numbers:
n = 100;
mu = 3.0;
sigma =2.0;
x = randn(n,1)*sigma+mu;
fprintf('standard deviation of the data is %.2f\n', std(x))
fprintf('              mean of the data is %.2f\n', mean(x))

% mean as parameter:
pmus = 2.0:0.01:4.0;
% matrix with the probabilities for each x and pmus:
lms = zeros(length(x), length(pmus));
for i=1:length(pmus)
    pmu = pmus(i);
    p = exp(-0.5*((x-pmu)/sigma).^2.0)/sqrt(2.0*pi)/sigma;
    lms(:,i) = p;
end
lm = prod(lms, 1);          % likelihood
loglm = sum(log(lms), 1);   % log likelihood

% position of maxima:
maxlm = max(lm);                        % height of the maximum
pmumaxlm = pmus(lm==maxlm);             % pmu where lm is at maximum
pmumaxloglm = pmus(loglm==max(loglm));  % same for loglm in one line
fprintf('      maximum of likelihood is at %.2f\n', pmumaxlm)
fprintf('  maximum of log-likelihood is at %.2f\n', pmumaxloglm)

% plot likelihood of mean:
subplot(1, 2, 1);
plot(pmus, lm);
xlabel('mean')
ylabel('likelihood')
subplot(1, 2, 2);
plot(pmus, loglm);
xlabel('mean')
ylabel('log likelihood')