%
%    height_example.m
%
%  Note: Code requires kde.m and histnorm.m, which can be downloaded fro MATLAB Central.

  clear all
  close all

%
% Input number of students n, true mean mu and standard deviation sigma, and compute heights
% from normal distribution N(mu,sigma^2). 
%
  
  n = input('n = ');
  mu = 67;              % 5'7'' 
  sigma = 2.5;          % 2 sigma = 5"

  heights = mu*ones(1,n) + sigma*randn(1,n);
  students = 1:1:n;

%
% Plot the heights and 2 sigma intervals.
%

  figure(1) 
  hold on
  plot(students,heights,'bx','linewidth',3);
  plot(students,mu*ones(1,n),'k-',students,(mu+2*sigma)*ones(1,n),'r--',students,(mu-2*sigma)*ones(1,n),'r--','linewidth',3);
  hold off
  axis([0 n 60 74])
  set(gca,'Fontsize',[20]);
  legend('Height','2\sigma Intervals','Location','Northeast')
  xlabel('Student')
  ylabel('Height')

%
% Plot histogram and kernel density estimate (kde) for x_i and compare to the original normal distribution
% N(mu,sigma^2).
%

  dx = (mu+3*sigma - mu+3*sigma)/1000;
  x = mu-(3*sigma):dx:mu+(3*sigma);           % Construct grid with +- 3 sigma
  pdf = normpdf(x,mu,sigma);                  % Compute N(mu,sigma^2) on grid x

  [bandwidth_heights,density_heights,heights_mesh,cdf_heights]=kde(heights);

  nbins = 15;
  figure(2)
  h = histnorm(heights,nbins,'plot');         % Compute and plot the histogram with nbins
  % histobj=histogram(heights,nbins,'normalization','pdf'); 
  hold on
  plot(x,pdf,'linewidth',3)
  plot(heights_mesh,density_heights,'r--','linewidth',3)
  hold off
  set(gca,'Fontsize',[20]);
  legend('Histogram','N(\mu,\sigma^2)','KDE','Location','Northeast')

%
% Compute the sample mean and standard deviation estimates and plot the sampling distribution for xbar.
%

  xbar = (1/n)*sum(heights)
  S2 = (1/(n-1))*sum((heights - xbar).^2);
  S = sqrt(S2)

  dx = (mu+(3*sigma/sqrt(n)) - mu+(3*sigma/sqrt(n)))/1000;
  x = mu-(3*sigma/sqrt(n)):dx:mu+(3*sigma/sqrt(n));           % Construct grid with +- 3 sigma/sqrt(n)
  pdf = normpdf(x,mu,sigma/sqrt(n));

  figure(3)
  plot(x,pdf,'linewidth',3);

%
% Compute 95% confidence limits and compare with 2*sigma/sqrt(n).
%

  alpha = 0.05;
                                  % Online T-value calculator: tcrit2 = 2.01669219 when n = 44; 
  tcrit2 = tinv(1-alpha/2,n-1)    % MATLAB t-inverse cummulative function 
  value= tcrit2*S/sqrt(n)         % Defines interval when sigma is not known

  2*sigma/sqrt(n)                 % Defines interval when sigma is known

%
%
%