% IrisBayesRiskUniformErrorRate.m
%
% This code calculates the Bayes risk error rate associated with the statistical decision 
%  approach for the classification problem example using of the Iris data
% used in R. A. Fisher's 1936 paper that developed classification 
% with calculations based only on measurement of sepal length and a uniform 
%  distribution
% 
% Make sure that the file with the Iris data - RawIrisDatamatrix.txt - is
% in the directory where you are running Matlab then opemn Matlab and run this 
%  .m Matlab file that is in the same directory by typing IrisBayesRisk
%
% load the data and put in in the matrix y
%    here the data have rows for each individual flower measured (in cm)
%    sepal length   sepal width   petal length   petal width  species
%   where species = 1 for Iris sets, 2 for Iris versicolor, 3 for Iris virginica
%   and there are 50 measurements for each species
%
load RawIrisDatamatrix.txt;
y=RawIrisDatamatrix;
%
%   set up the distribution of sepal length for each species
%   as a uniform distribution with the mean sepal length 
%   and +/- 1/2 of the range 
%
a1=mean(y(1:50,1))-.5*range(y(1:50,1));
b1=mean(y(1:50,1))+.5*range(y(1:50,1));
a2=mean(y(51:100,1))-.5*range(y(51:100,1));
b2=mean(y(51:100,1))+.5*range(y(51:100,1));
a3=mean(y(101:150,1))-.5*range(y(101:150,1));
b3=mean(y(101:150,1))+.5*range(y(101:150,1));
dens1=1./(b1-a1);
dens2=1./(b2-a2);
dens3=1./(b3-a3);
%
% define the uniform distributions as pr1, pr2, and pr3
%   for the three species 
%   and the posterior probability distributions which
%   are the Bayes risk functions as br1, br2 and br2
%
syms pr1(z)
pr1(z)=piecewise(z<a1,0, a1<=z<=b1, dens1, z>b1, 0);
syms pr2(z)
pr2(z)=piecewise(z<a2,0, a2<=z<=b2, dens2, z>b2, 0);
syms pr3(z)
pr3(z)=piecewise(z<a3,0, a3<=z<=b3, dens3, z>b3, 0);
syms b1(z)
   br1(z)=pr1(z)/(pr1(z)+pr2(z)+pr3(z));
syms b2(z)
   br2(z)=pr2(z)/(pr1(z)+pr2(z)+pr3(z));
syms b3(z)
   br3(z)=pr3(z)/(pr1(z)+pr2(z)+pr3(z));
%
% loop through each class to find error rate where ic counts 
% the number of correct classifications in each loop
%  and icount sums all correct classifications
%
icount=0;
ic=0;
for n=1:50
seplen=y(n,1);
if (br1(seplen) >= br2(seplen)) && (br1(seplen) >= br3(seplen)) 
     ic=ic+1;
   else
   end;
end;
icount=icount+ic;
ic=0;
for n=51:100
seplen=y(n,1);
if (br2(seplen) >= br1(seplen)) && (br2(seplen) >= br3(seplen)) 
     ic=ic+1;
   else
   end;
end;
icount=icount+ic;
ic=0;
for n=101:150
seplen=y(n,1);
if (br3(seplen) >= br1(seplen)) && (br3(seplen) >= br2(seplen)) 
     ic=ic+1;
   else
   end;
end;
icount=icount+ic;
err=(150.-icount)/150.
disp('training set error rate or empirical Bayes Risk is')
err



   