%regression-RVMsparse.m
%Relevant Vector Machine 
%(sparsity measurment technique) 
%Bayes regression 
%Training set 
clear all 
close all
v=0.5; 
scal=1;
x=-pi:0.1:pi;
betaactual=10;
t=sin(x)+randn(1,length(x))*(1/sqrt(betaactual));
TRAINx=x(1:2:length(x));
TRAINt=t(1:2:length(x))*scal;
TESTx=x(2:2:length(x));
TESTt=t(2:2:length(x))*scal;
%Estimate alpha and beta %Construncting the matrix PHI 
MEAN=linspace(-pi,pi,100);
PHI=[];
for i=1:1:length(TRAINx) 
    temp=[1]; 
    for j=1:1:length(MEAN)
        temp=[temp gbf(TRAINx(i),MEAN(j),v)];
    end
    PHI=[PHI;temp];
end
%The alpha vector and beta is solved using the sparsity analysis as given %below 
beta=1; 
alpha=ones(1,size(MEAN,2)+1)*inf; 
alpha(1)=rand*0.1; MODEL=[1];
Betacol=[beta]; 
for iteration=1:1:10 
    for k=1:1:size(MEAN,2)+1 temp=0; 
        for i=1:1:length(size(MEAN,2)+1) 
            if(i~=k) temp=temp+(1/alpha(i))*PHI(:,i)*PHI(:,i)';
            end
            temp=temp+(1/beta)*eye(size(temp,2));
        end
        s(k)=PHI(:,k)'*inv(temp)*PHI(:,k);
        q(k)=PHI(:,k)'*inv(temp)*TRAINt';
    end
    for k=1:1:size(MEAN,2)+1 
        if(~(isnan(q(k))|(isnan(s(k))))) 
            if(((q(k)^2)>s(k))&(alpha(k)<inf))
                alpha(k)=(s(k)^2)/((q(k)^2)-s(k));
            elseif((q(k)^2)>s(k))&((alpha(k)==inf))
                alpha(k)=(s(k)^2)/((q(k)^2)-s(k)); 
                if(length(find((MODEL-k)==0))==0) 
                    MODEL=[MODEL k]; 
                end
            elseif(((q(k)^2)<=s(k))&((alpha(k)<inf)))
                alpha(k)=inf;
                MODEL=nonzeros(MODEL-k)+k;
                MODEL=MODEL';
            end
        end
    end
    gamma=0;
    A=0;
    A=inv(diag(alpha(MODEL))+beta*PHI(:,MODEL)'*PHI (:,MODEL));
    mN=beta*A*PHI(:,MODEL)'*TRAINt';
    for i=1:1:length(MODEL)
        gamma(MODEL(i))=1-alpha(MODEL(i))*A(i,i);
    end
    s1=sum(([TRAINt- (PHI(:,MODEL)*mN)']).^2); 
    beta=inv((s1/(length(TRAINx)-sum(gamma))));
    Alphacol{iteration}=alpha;
    Betacol=[Betacol beta];
end
A=inv(diag(alpha(MODEL))+beta*PHI(:,MODEL)'*PHI (:,MODEL));
mN=beta*A*PHI(:,MODEL)'*TRAINt';
%Check with test data 
PHItest=[]; 
for i=1:1:length(TESTx) 
    temp=[1]; 
    for j=1:1:length(MEAN) 
        temp=[temp gbf(TESTx(i),MEAN(j),v)]; 
    end
    PHItest=[PHItest;temp];
end
y_predict=PHItest(:,MODEL)*mN;
y_actual=scal*sin(x(2:2:length(x)));
plot(y_predict,'r')
hold on 
plot(y_actual,'b')
plot(TESTt,'k')
figure 
plot(Betacol,'k')
title('Convergence of beta')

