%RVMforclassification.m
close all 
clear all
%RVM sparsity for classification %Multi-class problem
m1=[1 0 0]';
m2=[0 1 0]'; 
m1=rand(3,1)+m1;
m2=rand(3,1)+m2;
N1=10;
N2=10;
POSTRAIN=[];
POSTEST=[];
POSVAL=[];
temp1=rand(length(m1),1);
C1=temp1*temp1'*0.2;
temp2=rand(length(m2),1);
C2=temp2*temp2'*0.2;
[X1,C1est]=genrandn(m1,C1,N1);
[X2,C2est]=genrandn(m2,C2,N2); 
DATA=[X1 X2]; 
figure(1)
plot3(X1(1,:),X1(2,:),X1(3,:),'r*')
hold on 
plot3(X2(1,:),X2(2,:),X2(3,:),'b*')
[V,M]=kmeans(DATA',4);
B1=[]; B2=[]; B3=[]; B4=[]; 
sigma=0.5; 
for i=1:1:size(DATA,2)
    B1=[B1 gbf(DATA(:,i)',M(1,:),sigma)];
    B2=[B2 gbf(DATA(:,i)',M(2,:),sigma)];
    B3=[B3 gbf(DATA(:,i)',M(3,:),sigma)]; 
    B4=[B4 gbf(DATA(:,i)',M(4,:),sigma)]; 
end
%Formulation of matrix 
MAT=[ones(20,1) B1' B2' B3' B4'];
TRAIN=[MAT(1:1:5,:); MAT(11:1:15,:)];
TEST=[MAT(6:1:10,:); MAT(16:1:20,:)];
%Initialize A matrix 
TRAININDEX=[ones(1,5) zeros(1,5);]
TESTINDEX=[ones(1,5) zeros(1,5);] 
%Initialize weights 
%Computation of ynk (nth vector in the kth class) 
%Initializing W and alpha 
%W=0.01*(rand(5,1)*2-1);
%alpha=100*rand(1,5);;
%save INITIALIZE W alpha
load INITIALIZE
%load INITIALIZE 
ERROR=[]; 
for iteration=1:1:10 
    %Soft max computation for the individual classses 
    P1=(1./(1+exp(-TRAIN*W))); 
    P2=1-P1; 
    %Error function computation
P=[P1 P2] 
error=-1*sum(sum([ones(1,5) zeros(1,5);zeros(1,5) ones(1,5)]'.*log(P))); 
ERROR=[ERROR error];
%Updating the weights are obtained as follows. 
%Initialize alphai’s matrix 
A=diag(alpha); 
W=(inv(A))*TRAIN'*(TRAININDEX'-P1);
B=diag([P1.*(1-P1)]);
C=pinv(TRAIN'*B*TRAIN+A);
%Updata W 
for i=1:1:5 
    gamma(i)=1-alpha(i)*C(i,i); 
end
alpha=(gamma)./(W'.^2); 
end
%Checking with test data %Soft max computation for the individual classses 
P1=(1./(1+exp(-TEST*W))); 
P2=1-P1; 
%Error function computation 
P=[P1 P2];
[u,v]=min(P');
REF=[ones(1,5) ones(1,5)*2];
POS=length((find(REF-v)==0))*10 
figure(2)
plot(P(:,1),'r')
hold on
plot(P(:,2),'b')
figure(3)
plot(ERROR)