% IrisBayesRiskTriangle.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
%
%  This version uses a triangular distribution for each of the Sepal length 
%  distributions of the three Iris species
% 
% 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 triangle distribution with the mean sepal length 
%   and +/- 1/2 of the range as the ends of the distributions
%   mi holds the means, ai the left end, bi the right end
%   and hi is the height at the center of the distributions 
%
m1=mean(y(1:50,1));
a1=m1-.5*range(y(1:50,1));
b1=m1+.5*range(y(1:50,1));
m2=mean(y(51:100,1));
a2=m2-.5*range(y(51:100,1));
b2=m2+.5*range(y(51:100,1));
m3=mean(y(101:150,1));
a3=m3-.5*range(y(101:150,1));
b3=m3+.5*range(y(101:150,1));
h1=1./(b1-a1);
h2=1./(b2-a2);
h3=1./(b3-a3);
%
% define the triangle 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<=m1, h1*(z-a1)/(m1-a1), m1<=z<=b1, h1*(b1-z)/(b1-m1), z>b1, 0);
syms pr2(z)
pr2(z)=piecewise(z<a2,0, a2<=z<=m2, h2*(z-a2)/(m2-a2), m2<=z<=b2, h2*(b2-z)/(b2-m2), z>b2, 0);
syms pr3(z)
pr3(z)=piecewise(z<a3,0, a3<=z<=m3, h3*(z-a3)/(m3-a3), m3<=z<=b3, h3*(b3-z)/(b3-m3), z>b3, 0);
syms br1(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



   