%Probabilistic approach to classification - Logistic approach 
%Using Newton Raphson method 
m1=[1 0 0]';
m2=[0 1 0]';
m3=[0 0 1]';
m1=rand(3,1)+m1;
m2=rand(3,1)+m2;
m3=rand(3,1)+m3;
N1=1000;
N2=1000;
N3=1000;
POSTRAIN=[];
POSTEST=[];
POSVAL=[];
POSTRAINH=[];
POSTESTH=[];
POSVALH=[];
ERROR=[];
ERRORH=[]; 
temp1=rand(length(m1),1); 
C1=temp1*temp1'*0.2;
temp2=rand(length(m2),1);
C2=temp2*temp2'*0.2;
temp3=rand(length(m3),1);
C3=temp3*temp3'*0.2;
[X1,C1est]=genrandn(m1,C1,N1);
[X2,C2est]=genrandn(m2,C2,N2);
[X3,C3est]=genrandn(m3,C3,N3);
DATA=[X1 X2 X3]; 
[V,M]=kmeans(DATA',3);
sigma=1:0.5:5;
sigma=1;
for trial=1:1:length(sigma) 
B1=[];
B2=[];
B3=[];
for i=1:1:size(DATA,2)
    B1=[B1 gbf(DATA(:,i)',M(1,:),sigma(trial))]; 
    B2=[B2 gbf(DATA(:,i)',M(2,:),sigma(trial))];
    B3=[B3 gbf(DATA(:,i)',M(3,:),sigma(trial))];
end
%Formulation of matrix 
MAT=[B1' B2' B3'];
TRAIN=[MAT(1:1:400,:); 
MAT(1001:1:1400,:); 
MAT(2001: 1:2400,:)];
TRAININDEX=zeros(1200,3); 
TRAININDEX(1:400,1)=1; 
TRAININDEX(401:800,2)=1; 
TRAININDEX(801:1200,3)=1; 
VAL=[MAT(401:1:700,:); 
MAT(1401:1:1700,:); 
MAT(2401: 1:2700,:)]; 
VALINDEX=zeros(900,3); 
VALINDEX(1:300,1)=1; 
VALINDEX(301:600,2)=1; 
VALINDEX(601:900,3)=1; 
%Initializing the weights 
W=rand(3,3); 
WH=W; 
%Soft max computation for the individual classses 
P1=exp(TRAIN*W(:,1))./(exp(TRAIN*W(:,1))+exp(TRAIN*W(:,2))+exp(TRAIN*W(:,3)));
P2=exp(TRAIN*W(:,2))./(exp(TRAIN*W(:,1))+exp(TRAIN*W(:,2))+exp(TRAIN*W(:,3)));
P3=exp(TRAIN*W(:,3))./(exp(TRAIN*W(:,1))+exp(TRAIN*W(:,2))+exp(TRAIN*W(:,3)));
%Error function computation 
P=[P1 P2 P3]; 
error=-1*sum(sum(TRAININDEX.*log(P)));
ERROR=[];
for iteration=1:1:100
    PH=P;
    I=diag([1 1 1]);
    temp=0;
    H=[];
    for r=1:1:3 
        for s=1:1:3
            temp=0;
            for n=1:1:1200
                temp=temp+TRAIN(n,:)'*TRAIN(n,:)*PH(n,r)*(I(r,s)-PH(n,s));
            end
            HESSIAN{r,s}=-temp;
        end
    end
    WH(:,1)=WH(:,1)-0.01*pinv(HESSIAN{1,1})*TRAIN'*(P(:,1) -TRAININDEX(:,1));
    WH(:,2)=WH(:,2)-0.01*pinv(HESSIAN{2,2})*TRAIN'*(P(:,2) -TRAININDEX(:,2));
    WH(:,3)=WH(:,3)-0.01*pinv(HESSIAN{3,3})*TRAIN'*(P(:,3) -TRAININDEX(:,3)); 
    %Soft max computation for the individual classses 
    PH1=exp(TRAIN*WH(:,1))./(exp(TRAIN*WH(:,1))+ exp(TRAIN*WH(:,2))+exp(TRAIN*WH(:,3))); 
    PH2=exp(TRAIN*WH(:,2))./(exp(TRAIN*WH(:,1))+exp(TRAIN*WH(:,2))+exp(TRAIN*WH(:,3))); 
    PH3=exp(TRAIN*WH(:,3))./(exp(TRAIN*WH(:,1))+exp(TRAIN*WH(:,2))+exp(TRAIN*WH(:,3)));
    %Error function computation 
    PH=[PH1 PH2 PH3]; 
    error=-1*sum(sum(TRAININDEX.*log(P))); 
    ERRORH=[ERRORH error]; 
    W(:,1)=W(:,1)-0.01*TRAIN'*(P(:,1)-TRAININDEX(:,1));
    W(:,2)=W(:,2)-0.01*TRAIN'*(P(:,2)-TRAININDEX(:,2));
    W(:,3)=W(:,3)-0.01*TRAIN'*(P(:,3)-TRAININDEX(:,3));
    %Soft max computation for the individual classses 
    P1=exp(TRAIN*W(:,1))./(exp(TRAIN*W(:,1))+exp(TRAIN*W(:,2))+exp(TRAIN*W(:,3)));
    P2=exp(TRAIN*W(:,2))./(exp(TRAIN*W(:,1))+exp(TRAIN*W(:,2))+exp(TRAIN*W(:,3)));
    P3=exp(TRAIN*W(:,3))./(exp(TRAIN*W(:,1))+exp(TRAIN*W(:,2))+exp(TRAIN*W(:,3)));
    %Error function computation 
    P=[P1 P2 P3]; 
    error=-1*sum(sum(TRAININDEX.*log(P))); 
    ERROR=[ERROR error]; 
end
figure 
plot(ERRORH,'r') 
hold on 
plot(ERROR,'b')
[V,I]=max(PH');
REF=[ones(1,400) ones(1,400)*2 ones(1,400)*3];
[R,L]=find((I-REF)==0);
POSTRAINH=[POSTRAINH (length(L)/length(I))*100];
PH1=exp(VAL*WH(:,1))./(exp(VAL*WH(:,1))+exp(VAL*WH(:,2))+exp(VAL*WH(:,3))); 
PH2=exp(VAL*WH(:,2))./(exp(VAL*WH(:,1))+exp(VAL*WH(:,2))+exp(VAL*WH(:,3)));
PH3=exp(VAL*WH(:,3))./(exp(VAL*WH(:,1))+exp(VAL*WH(:,2))+exp(VAL*WH(:,3)));
PH=[PH1 PH2 PH3]; 
[V,I]=max(PH'); 
REF=[ones(1,300) ones(1,300)*2 ones(1,300)*3];
[R,L]=find((I-REF)==0);
POSVALH=[POSVALH (length(L)/length(I))*100];
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%% 
[V,I]=max(P');
REF=[ones(1,400) ones(1,400)*2 ones(1,400)*3];
[R,L]=find((I-REF)==0);
POSTRAIN=[POSTRAIN (length(L)/length(I))*100];
P1=exp(VAL*W(:,1))./(exp(VAL*W(:,1))+exp(VAL*W(:,2))+exp(VAL*W(:,3)));
P2=exp(VAL*W(:,2))./(exp(VAL*W(:,1))+exp(VAL*W(:,2))+exp(VAL*W(:,3)));
P3=exp(VAL*W(:,3))./(exp(VAL*W(:,1))+exp(VAL*W(:,2))+exp(VAL*W(:,3)));
P=[P1 P2 P3];
[V,I]=max(P');
REF=[ones(1,300) ones(1,300)*2 ones(1,300)*3];
[R,L]=find((I-REF)==0); 
POSVAL=[POSVAL (length(L)/length(I))*100];
end
[V,I]=max(POSVAL);
%Using testing data 
B1=[]; B2=[]; B3=[]; 
for i=1:1:size(DATA,2)
    B1=[B1 gbf(DATA(:,i)',M(1,:),sigma(I))];
    B2=[B2 gbf(DATA(:,i)',M(2,:),sigma(I))];
    B3=[B3 gbf(DATA(:,i)',M(3,:),sigma(I))];
end
%Formulation of matrix
MAT=[B1' B2' B3'];
TEST=[MAT(701:1:1000,:); MAT(1701:1:2000,:); MAT(2701:1:3000,:)];
TESTINDEX=zeros(1200,3);
TESTINDEX(1:300,1)=1;
TESTINDEX(301:600,2)=1; 
TESTINDEX(601:900,3)=1;
P1=exp(TEST*W(:,1))./(exp(TEST*W(:,1))+exp(TEST*W(:,2))+exp(TEST*W(:,3)));
P2=exp(TEST*W(:,2))./(exp(TEST*W(:,1))+exp(TEST*W(:,2))+exp(TEST*W(:,3))); 
P3=exp(TEST*W(:,3))./(exp(TEST*W(:,1))+exp(TEST*W(:,2))+exp(TEST*W(:,3))); 
P=[P1 P2 P3]; 
[V,I]=max(P'); 
REF=[ones(1,300) ones(1,300)*2 ones(1,300)*3]; 
[R,L]=find((I-REF)==0); 
POSTEST=[length(L)/length(I)]*100; 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
PH1=exp(TEST*WH(:,1))./(exp(TEST*WH(:,1))+exp(TEST*WH(:,2))+exp(TEST*WH(:,3))); 
PH2=exp(TEST*WH(:,2))./(exp(TEST*WH(:,1))+exp(TEST*WH(:,2))+exp(TEST*WH(:,3))); 
PH3=exp(TEST*WH(:,3))./(exp(TEST*WH(:,1))+exp(TEST*WH(:,2))+exp(TEST*WH(:,3))); 
PH=[PH1 PH2 PH3]; 
[V,I]=max(PH'); 
REF=[ones(1,300) ones(1,300)*2 ones(1,300)*3]; 
[R,L]=find((I-REF)==0); 
POSTESTH=[length(L)/length(I)]*100; 
figure 
subplot(2,1,1) 
plot(ERROR) 
title('Convergence using gradient technique') 
subplot(2,1,2) 
plot(ERRORH) 
title('Convergence using Newton Raphson technique')
    
