%
%  likelyhood ratio laptop assignment
%
%  m_s and s_s are the mean and standard deviations for sick
%  m_n and s_b are the mean and standard deviations for normal
%  K_n and K_s are the number of bins to use in generating the pdfs
%  N_pdf is the number of sample points used to generate the data for each
%  N_pts is the total number of patients to test
%  p is the probability of being normal
%
function laptop12(m_s,s_s,K_s,m_n,s_n,K_n,p,N_pdf,N_pts)
%
%  generate data with the appropriate means and variances
%
  D_s = m_s+s_s*randn(1,N_pdf);
  D_n = m_n+s_n*randn(1,N_pdf);
%
%  plot histograms of the data
%
  figure; 
  orient tall
  subplot(2,1,1); hist(D_n); xlabel('z');  title('Normals');
  subplot(2,1,2); hist(D_s); xlabel('z'); title('Sick');
%
%  Get the density for sick
%
  mu_s = mean(D_s);
  sigma_s = std(D_s);
%
%  create and anonymous function for pdf sick
%
  pdf_s = @(z) exp(-(z-mu_s).^2/(2*sigma_s^2))/sqrt(2*pi*sigma_s^2);
%  
  a_s = min(D_s);
  b_s = max(D_s);
%   dz = ??
%   ni = ??
%   pdf = ??
  z = linspace(a_s,b_s,1000);
%  
%  The following sum should equal 1
%
  sum(pdf*(dz(2)-dz(1)))
%
%  plot the estimated densities
%
  figure;
  orient tall
  subplot(3,1,1); plot(dz,pdf,'x',z,pdf_s(z),'-');
  xlabel('z'); ylabel('f_S(z)'); grid; 
  title([' \mu = ' num2str(mu_s,3) ', \sigma = ' num2str(sigma_s,3)]); 
%
%  Get the density for normal
%
  mu_n = mean(D_n);
  sigma_n = std(D_n);
%
%  create and anonymous function for pdf normal
%
  pdf_n = @(z) exp(-(z-mu_n).^2/(2*sigma_n^2))/sqrt(2*pi*sigma_n^2);
%  
  a_n = min(D_n);
  b_n = max(D_n);
%   dz = ??
%   ni = ??
%   pdf = ??
  z = linspace(a_n,b_n,1000);
%  
%  The following sum should equal 1
%
  sum(pdf*(dz(2)-dz(1)))
%
%  plot the estimated densities
%
  subplot(3,1,2); plot(dz,pdf,'x',z,pdf_n(z),'-');
  xlabel('z'); ylabel('f_N(z)'); grid; 
  title([' \mu = ' num2str(mu_n,3) ', \sigma = ' num2str(sigma_n,3)]); 
%
%  Now determine the threshold eta
%
% css and cnn are the costs of making the correct choices
% cns is the cost of choosing normal when sick is true
% csn is the cost of choosing sick when normal is true
%
  cns = 1; csn = 1; css = 0; cnn = 0;
  eta = (1-p)*(cns-css)/(p*(csn-cnn));
%
%  assuming Gausian densities, we can find the threshold
% locations
%  
  Pn = -[1 -2*mu_n mu_n^2]/(2*sigma_n^2);
  Ps =  [1 -2*mu_s mu_s^2]/(2*sigma_s^2);
  Peta = [0 0 -log(sigma_n*eta/sigma_s)];
  P = Pn+Ps+Peta
  disp('Thresholds at :');
  r = roots(P)
%
%  Now process some new patients
%
%  get the normal patients
%
  N_n = N_pts*p;
  z_n = m_n+s_n*randn(1,N_n);
%
%  compare them to the test
%
%   lambda = ??
%   n_correct = sum(??);
  n_false = N_n - n_correct;
%
%  get the sick patients
%
  N_s = N_pts - N_n;
  z_s = m_s + s_s*randn(1,N_s);
%
%  compare them to the test
%
%   lambda = ??
%   s_correct = sum(??);
  s_false = N_s - s_correct;
%
%  prepare to plot the two densities on one graph
%
  a = min(a_s,a_n); b = max(b_s,b_n);
  z = linspace(a,b,1000);
%
  subplot(3,1,3); plot(z,pdf_s(z),'-',z,pdf_n(z),'--',r(1)*[1 1],[0 0.25],'r-',r(2)*[1 1],[0 0.25],'r-'); grid;
  xlabel('z');
  title([' p_n = ' num2str(p) ', N_n = ' num2str(N_n) ' (' num2str(n_correct/N_n,3) ', '  num2str(n_false/N_n,3) ')' ...
                                ',  N_s = ' num2str(N_s) ' (' num2str(s_correct/N_s,3) ', '  num2str(s_false/N_s,3) ')' ]); 
  legend('f_S','f_N');
  
