% IrisBayesRiskTriangle.m
%
% This code calculates the Bayes risk 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;
%
%   now plot histograms of sepal length for the species 
%       Iris setosa, Iris versicola, Iris virginica which are 
%       plotted in that order
%
figure
hist(y(1:50,1),10);
figure
hist(y(51:100,1),10);
figure
hist(y(101:150,1),10);
%
%   now do a scatterplot of sepal length and sepal width 
%   for each of the 3 species
%
figure
xlim([4 8]);
ylim([1.5 4.5]);
plot(y(1:50,1),y(1:50,2),'+');
hold
plot(y(51:100,1),y(51:100,2),'*');
plot(y(101:150,1),y(101:150,2),'o');
hold off
%
%   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);
%
%  graph the triangle distributions for sepal length for 
% the three species
%
figure
xlim([4 8]);
ylim([0 1]);
plot([a1,m1],[0,h1],'r',[m1,b1],[h1,0],'r');
hold
plot([a2,m2],[0,h2],'k',[m2,b2],[h2,0],'k');
plot([a3,m3],[0,h3],'b',[m3,b3],[h3,0],'b');
hold off
%
% 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));
%   
%  graph the three Bayes risk functions
%
figure
xlim([4 8]);
ylim([0 1]);
fplot(br1,[4,8]);
hold
fplot(br2,[4,8]);
fplot(br3,[4,8]);
%
%  calculate the three Bayes risk functions for particular 
%   values of sepal length and classify an individual with 
%   that sepal length as the species with maximum Bayes risk
%
s=1;
while s>0;
yy=input('sepal length at which to calculate the 3 Bayes risks:    ')
r1=eval(br1(yy));
r2=eval(br2(yy));
r3=eval(br3(yy));
disp('values of Bayes Risk for Iris setosa, Iris versicola, Iris virginica are:')
r1,r2,r3
s=input('Do you want to stop - if so enter 0, if not enter 1')
end



   