% filename: inference_example_3.m

% This script plots the likelihood function for Example 3 in the Introduction to Inference presentation.  In the example, five randomly sampled MLB baseball players' heights in inches were 73, 76, 74, 74, and 73.  Assuming them to have come from a normal distribution with mean mu and standard deviation sigma, the goal is to give plausible values for mu and sigma.

% **********

% create the data vector and vectors of candidates for mu and sigma:


data = [73, 76, 74, 74, 73];

mu = [70:0.1:80];
sigma = [0:0.1:5];

% create a matrix of likelihoods, one for each combination of mu and sigma in their vectors.  We accomplish that here by first creating an appropriately-sized matrix of zeros, then filling in the likelihoods one by one using embedded for loops.

L = zeros(length(mu), length(sigma));

for i = 1:length(mu)
  for j = 1:length(sigma)
    L(i,j) = prod(normpdf(data, mu(i), sigma(j)));
  end
end

% make a surface plot of the likelihood function.  Note that we're plotting L' instead of L.  That's a side effect of two inconsistencies with matrices: when indexing them, you always name the row first and then the column; but when you want to plot them, you (might often) want to plot the rows as "y's" and the columns as "x's".

figure
surf(mu, sigma, L')
xlabel('\mu')
ylabel('\sigma')

% make a contour plot of the likelihood function with level curves drawn at 10%, 5%, and 1% of the maximum likelihood.  Also plot the MLE as a red circle ('ro').  

figure
maxL = max(max(L));   % (the max function applied only once to a matrix returns a vector of column maxima, so we apply it twice to find the overall matrix maximum)
levels = maxL*[.1 .05 .01];
contour(mu, sigma, L', levels)
xlabel('\mu')
ylabel('\sigma')
title('MLB player heights, \mu and \sigma')
[index_mu, index_sigma] = find(L == maxL);
hold on
plot(mu(index_mu), sigma(index_sigma),'ro')
hold off