%
% Examples 9.7
% 
% Written by:
% -- 
% John L. Weatherwax                2008-02-20
% 
% email: wax@alum.mit.edu
% 
% Please send comments and especially bug reports to the
% above email address.
% 
%-----

clear all; 
close all; 
clc; 

addpath('../../Code/CSTool');

n = 100;   % <- the size of our random sample
M = 100;   % <- the number of Monte-Carlo estimates 
all_pcc = zeros(1,M); 

for mi=1:M,
  % draw some samples from the proposed densities: 
  % 
  cl = rand(1,n); 
  lab1 = find( cl<=0.6 ); lab2 = setdiff( 1:n, lab1 );
  n1   = length(lab1);      n2 = length(lab2); 
  c1 = normrnd(-1,1,1,n1); 
  c2 = normrnd(+1,1,1,n2); 
  
  ncc=0; 
  
  % Loop first through all of the patterns corresponding to the first class.
  % Here correct classification is obtained if p1_hat*pxgc1 > p2_hat*pxgc2;
  p1_hat = n1/n;      
  muc1   = mean(c1);  % These will be the same for this part ...
  varc1  = var(c1);
  for i = 1:n2
    % Get the test point and the training set
    c2_train = c2;
    % This is the testing point.
    x = c2_train(i);
    % Delete this point from training set ... the result is the training set.
    c2_train(i)=[]; 
    p2_hat = (n2-1)/n; 
    muc2   = mean(c2_train);
    varc2  = var(c2_train);
    pxgc1=normpdf(x,muc2,sqrt(varc2));
    pxgc2=normpdf(x,muc1,sqrt(varc1));
    if p1_hat*pxgc1 > p2_hat*pxgc2 % then we correctly classified it
      ncc = ncc+1;
    end
  end

  % Loop through all of the patterns corresponding to the second class.
  % Here correct classification is obtained if p2_hat*pxgc2 > p1_hat*pxxgver
  p2_hat = n2/n; 
  muc2   = mean(c2);
  var2   = var(c2);
  % These remain the same for the following.
  for i = 1:n1
    % Get the test point and training set
    c1_train = c1;
    % This is the testing point.
    x = c1_train(i);
    % Delete this point from training set ... the result is the training set.
    c1_train(i)=[];
    muc1 = mean(c1_train);
    varc1 = var(c1_train);
    pxgc1 = normpdf(x,muc2,sqrt(varc2));
    pxgc2 = normpdf(x,muc1,sqrt(varc1));
    if pxgc2 > pxgc1 % then we correctly classified it
      ncc = ncc+1;
    end
  end

  pcc = ncc/n;
  %fprintf('pcc = %10.6f\n',pcc); 
  all_pcc(mi) = pcc; 
end

figure; hist( all_pcc ); 
xlabel( 'pcc value' ); ylabel( 'pcc count' ); 
title( sprintf('histogram of the p_{cc} NMC=%d',M) ); 
saveas( gcf, 'prob_9_11_pcc_hist', 'eps' );