%Bayesregression.m
%Bayes regression 
%Training set 
v=0.5; 
x=-pi:0.1:pi; 
beta=10; 
t=sin(x)+randn(1,length(x))*(1/sqrt(beta));
TRAINx=x(1:2:length(x)); 
TRAINt=t(1:2:length(x));
TESTx=x(2:2:length(x)); 
TESTt=t(2:2:length(x)); 
%Estimate alpha and beta 
%Let the number of basis function (Gaussian) be 20 
%and mean of the basis function is uniformly spaced between %-pi to pi 
MEAN=linspace(-pi,pi,20); 
%Construncting the matrix PHI 
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
%Initialize alpha and beta 
alpha=10; 
beta=1; 
Alphacol=[alpha]; 
Betacol=[beta]; 
for iteration=1:1:50 
    [E,D]=eig(beta*PHI'*PHI);
    A=(alpha*eye(length(MEAN)+1)+beta*PHI'*PHI);
    mN=beta*inv(A)*PHI'*TRAINt';
    gamma=0 
    for i=1:1:size(D,2)
        gamma=gamma+D(i,i)/(alpha+D(i,i));
    end
    alpha=gamma/(mN'*mN);
    s=0;
    for i=1:1:length(TRAINx)
        s=s+sum((TRAINt(i)-mN'*PHI(i,:)').^2)
    end
    beta=inv((1/(length(TRAINx)-gamma))*s);
    Alphacol=[Alphacol alpha];
    Betacol=[Betacol beta];
end
%Final mN 
A=(alpha*eye(length(MEAN)+1)+beta*PHI'*PHI); 
mN=beta*inv(A)*PHI'*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*mN;
y_actual=sin(x(2:2:length(x)));
plot(y_predict)
hold on
plot(y_actual)
plot(TESTt)
figure
subplot(2,1,1)
plot(Alphacol)
title('Convergence of alpha')
subplot(2,1,2)
plot(Betacol)
title('Convergence of beta')
