test_lssvmard_rbf.m

The goal of SPID is to provide the user with tools capable to simulate, preprocess, process and clas

% Generating training data
ndata = 100;
x1 = rand(ndata, 1);
x2 = x1 + 0.05*randn(ndata, 1);
x3 = 0.5 + 0.5*randn(ndata, 1);
x = [x1, x2, x3];
t = sin(2*pi*x1) + 0.1*randn(ndata, 1);

mx1 = mean(x1); sx1 = std(x1); x1norm = (x1 - mx1)/sx1;
mx2 = mean(x2); sx2 = std(x2); x2norm = (x2 - mx2)/sx2;
mx3 = mean(x3); sx3 = std(x3); x3norm = (x3 - mx3)/sx3;
Xtrnorm = [x1norm x2norm x3norm];
Ytrnorm = t;

% Generating test data
ndata=1000;
x1 = rand(ndata, 1); x1norm = (x1 - mx1)/sx1;
x2 = x1 + 0.05*randn(ndata, 1); x2norm = (x2 - mx2)/sx2;
x3 = 0.5 + 0.5*randn(ndata, 1); x3norm = (x3 - mx3)/sx3;
xtest = [x1, x2, x3]; ttest = sin(2*pi*x1) + 0.1*randn(ndata, 1);
Xtstnorm = [x1norm x2norm x3norm]; Ytest = ttest;

t_start=cputime;
N_inputs = size(Xtrnorm, 2);


%sigStart = 1; gamStart = 1;
nIter=N_inputs*20;
knl=1; % RBF kernel

% Initial gam and sig
dispOpt=0;
nlsig=4; knl=1; sigs=[]; % for RBF kernel; nlsig=4 four layers of sig seaching
%nlsig=1; knl=1; sigs=[39.3]; % for RBF kernel with sigma fixed
%nlsig=1; knl=2; sigs=[1:4]; % for polynomial kernel, sig -> the order
[sig_opt, gam_opt, L3_opt, neffPar, mu_opt, zeta_opt]=...
         optSigGamLssvmBay(Xtrnorm, Ytrnorm, sigs, nlsig, dispOpt,knl), 

     
Osig1=sig_opt; gam1=gam_opt;
%Osig1=0; gam1=1;

% Initialize the relevance vector Obarpil (proportional to 1/sig)
if Osig1<=1
    Obarpi1=0.999*ones(1,N_inputs);
else
    Obarpi1=1/Osig1*ones(1,N_inputs);
end

   
Obunc1=log(Obarpi1./(1-Obarpi1));

%Obarpi1=Obarpi2; gam1=gam2; 
%Obarpi1(1,1)=Osig1;
%for i=2:N_inputs,
%    Obarpi1(1,i)=Obarpi1(1,i-1)*0.8;
%end,

global fixGamma,
fixGamma=0, % Fix gamma during the iteration or not. 0- no, 1-yes

% unconstrained optimization
%[Obarpi2, Omu2, Ozeta2, OL32] = optim_ardL3r(Xtrnorm, Ytrnorm, Obarpi1 , gam1, nIter,knl);

%unconstrained optimization
[Obunc2, Omu2, Ozeta2, OL32] = optim_ardL3r_unc(Xtrnorm, Ytrnorm, Obunc1 , gam1, nIter,knl); 
Obarpi2=1./(1+exp(-Obunc2)); % calculating the relavance vector

et=(cputime-t_start)/60; % elapsed time
gam2=Ozeta2/Omu2, % optimal gamma

Obarpi2,  % optimal relevance 

% gamards: store the optimal results
gamards.theta=Obarpi2; gamards.gam=gam2; gamards.L3=OL32; 
gamards.mu=Omu2; gamards.zeta=Ozeta2; gamards.et=et;

% The new X should be transformed to 
for i = 1:length(gamards.theta),
  Xtrnormnew(:,i) = Xtrnorm(:,i)*gamards.theta(i);
end

% The new X should be transformed  
for i = 1:length(gamards.theta),
  Xtstnormnew(:,i) = Xtstnorm(:,i)*gamards.theta(i);
end

%
% Example for binary classification problem 
%
% Getting the binary classificatiton labels
Ytrn2class = t; Ytrn2class(t>0)=1; Ytrn2class(t<=0)=-1; 
Ytst2class = ttest; Ytst2class(ttest>0)=1; Ytst2class(ttest<=0)=-1;

% Training
lssvcBmod0=lssvcBay_train(Xtrnorm,Ytrn2class,1:size(Xtrnormnew,2),{'kernelType','rbf','maxsteps',-1});
% Test
[Ytrn_class0, Ytrn_probpos0]=lssvcBay_predict(Xtrnorm, lssvcBmod0);
[Ytst_class0, Ytst_probpos0]=lssvcBay_predict(Xtstnorm, lssvcBmod0);

[TN1,TP1,FP1,FN1,auctrn0,setrn0]=roc_tf(Ytrn_probpos0(Ytrn2class<=0),Ytrn_probpos0(Ytrn2class>0));   
[TN2,TP2,FP2,FN2,auctst0,setst0]=roc_tf(Ytst_probpos0(Ytst2class<=0),Ytst_probpos0(Ytst2class>0));   
acctrn0=sum(Ytrn_class0==Ytrn2class)/length(Ytrn2class);
acctst0=sum(Ytst_class0==Ytst2class)/length(Ytst2class);
fprintf('Modelling without weighting (on input variables)\n')
fprintf('Train: ACC=%.2f%% AUC=%.4f (se=%.4f)\n', acctrn0*100, auctrn0, setrn0);  
fprintf('Test: ACC=%.2f%% AUC=%.4f (se=%.4f)\n', acctst0*100, auctst0, setst0);

% Training
lssvcBmod=lssvcBay_train(Xtrnormnew,Ytrn2class,1:size(Xtrnormnew,2),{'kernelType','rbf','maxsteps',-1});
% Test
[Ytrn_class, Ytrn_probpos]=lssvcBay_predict(Xtrnormnew, lssvcBmod);
[Ytst_class, Ytst_probpos]=lssvcBay_predict(Xtstnormnew, lssvcBmod);

[TN1,TP1,FP1,FN1,auctrn,setrn]=roc_tf(Ytrn_probpos(Ytrn2class<=0),Ytrn_probpos(Ytrn2class>0));   
[TN2,TP2,FP2,FN2,auctst,setst]=roc_tf(Ytst_probpos(Ytst2class<=0),Ytst_probpos(Ytst2class>0));   
acctrn=sum(Ytrn_class==Ytrn2class)/length(Ytrn2class);
acctst=sum(Ytst_class==Ytst2class)/length(Ytst2class);
fprintf('Modelling using the ARD weighted inputs!\n')
fprintf('Train: ACC=%.2f%% AUC=%.4f (se=%.4f)\n', acctrn*100, auctrn, setrn);  
fprintf('Test: ACC=%.2f%% AUC=%.4f (se=%.4f)\n', acctst*100, auctst, setst);



return