optsiggamlssvmbay.m

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

function [sig_opt, gam_opt, L3_opt, neffPar, mu_opt, zeta_opt, sloga_opt, slogb_opt]=...
         optSigGamLssvmBay(Xtrain, Ytrain, Sigma0, nlsig, dispOpt,kenn), 

% function [sig_opt, gam_opt, L3_opt, neffPar, mu_opt, zeta_opt, sloga_opt, slogb_opt]=...
%         optSigGamLssvmBay(Xtrain, Ytrain, Sigma0, nlsig, dispOpt,kenn), 
% dispOpt: >0 display the evolution of sigma and gamma, and warnig message about optimization
%          =0 not display the evoltion of sigma and gamma, but disp the warning message
%          <0 display nothing
% kenn: kernel type: 1- rbf, 2 - linear

N_inputs = size(Xtrain, 2);

if nargin<3,
   Sigma0 = sqrt(N_inputs)*[0.1 1 2 3 5];
elseif isempty(Sigma0),
   Sigma0 = sqrt(N_inputs)*[0.1 1 2 3 5]; 
end,

if nargin<4
   nlsig=1;
end,
if nargin<5,
   dispOpt=0;
end,
if nargin<6,
   kenn=1;
end,

global logKK0,
N=length(Ytrain);
logKK0=zeros(N);
if kenn==1,
	for i=1:N,
	 for j=1:i,
	  logKK0(i,j) = -(Xtrain(i,:)-Xtrain(j,:))*(Xtrain(i,:)-Xtrain(j,:))';
	  logKK0(j,i) = logKK0(i,j);
	 end
	end
elseif kenn==2,
	for i=1:N,
	 for j=1:i,
	  logKK0(i,j) = Xtrain(i,:)*Xtrain(j,:)'+1;
	  logKK0(j,i) = logKK0(i,j); 
	 end
	end
elseif kenn<0,
    logKK0=[]; % additivive kernel  % add by LC Sept25,2004    
end,

numSteps=2; 

	N_inputs = size(Xtrain, 2);
    if dispOpt>0,
       	disp(['#inputs : ', num2str(N_inputs), ' x ', num2str(length(Ytrain))]),
    end,
    
	  Sigma = [];
	  msig=0;
      for lsig=1:nlsig,
        L3=[]; BPlikeli2=[];
        for i = 1:length(Sigma0),
         if Sigma0(i)~=msig,
            if dispOpt>0,
	   	      disp(['Sigma ', num2str(i), ': ', num2str(Sigma0(i))]),
            end,
	    	if i == 1,
      			startgamma = 1;
	      	 else 
   	   		startgamma = gammas(i-1);
		   	end

		   	[Bmu(i), Bzeta(i), Bnumeffpar(i), BPlikeli3(i), BPlikeli2(i), Bsloga(i), Bslogb(i),ew,ed,exitflag] ...
               = optim_gammaBayLSSVMardr(startgamma, Sigma0(i), Xtrain, Ytrain, kenn);
            if exitflag>0,
%                L3(i)= 1/2*BPlikeli2(i)-log(Bsloga(i))-log(Bslogb(i)); % used before Mar 24, 2004, 
%                L3(i)=-log(BPlikeli3(i));  % Revision on Mar 24, 2004, L3->-log(PosteriorLikelihood)
                L3(i)=-BPlikeli3(i);  
            else
                if dispOpt>=0,
                disp(['!!! warning (', num2str(exitflag), '): sig= ', num2str(Sigma0(i))]),
                end
                L3(i)=inf;                % Revision on Mar 24, 2004, L3->-log(PosteriorLikelihood)
%                L3(i)=inf;                 
%                L3(i)= 1/2*BPlikeli2(i)-log(Bsloga(i))-log(Bslogb(i));

             end, 
             gammas(i) = abs(Bzeta(i)/Bmu(i));
         else
            L3(i)=mL3; gammas(i)=mgam; 
            Bmu(i)=mmu; Bzeta(i)=mzeta; Bnumeffpar(i)=mnumeff;
            BPlikeli3(i)=mpl3; BPlikeli2(i)=mpl2;
            Bsloga(i)=msa; Bslogb(i)=msb;
             
         end,
        end, % for i=1:length(Sigma0),
      
    	if lsig==nlsig, 
        	 Sigma=Sigma0;  break, 
	    end,
        
        % store the output for the optimal sigma at this level,
        [mL3, imin]=min(L3);
        msig=Sigma0(imin); mgam=gammas(imin); 
        mmu=Bmu(imin); mzeta=Bzeta(imin); mnumeff=Bnumeffpar(imin);
        mpl3=BPlikeli3(imin); mpl2=BPlikeli2(imin); 
        msa=Bsloga(imin); msb=Bslogb(imin);
       
        % refine the range of sigs
        Sigmar=[Sigma0(imin)];
        if imin==1,
           Sigma0=[Sigma0(1)*0.5, Sigma0];
           imin=imin+1;
        elseif imin==length(Sigma0),
           Sigma0=[Sigma0, Sigma0(length(Sigma0))*1.5];
        end
        stepWidthLeft=(Sigma0(imin)-Sigma0(imin-1))/(numSteps+1);
        for j=1:numSteps,
            Sigmar=[Sigma0(imin)-j*stepWidthLeft, Sigmar];
        end
       
       stepWidthRight=(Sigma0(imin+1)-Sigma0(imin))/(numSteps+1);
       for j=1:numSteps,
           Sigmar=[Sigmar, Sigma0(imin)+j*stepWidthRight];
       end

       if (length(Sigmar)==1), break, end,
       Sigma0=Sigmar;       
       
     end  % for lsig        

     [L3_opt, iopt] = min(L3);
     sig_opt = Sigma(iopt);
     gam_opt = gammas(iopt);
     neffPar = Bnumeffpar(iopt);
     mu_opt=Bmu(iopt);
     zeta_opt=Bzeta(iopt);
     sloga_opt=Bsloga(iopt);
     slogb_opt=Bslogb(iopt);

     if dispOpt>0,
        figure, plot(L3),
    end,

return,