%predictivedistribution.m
%Prediction distribution 
v=0.5; 
x=-pi:0.1:pi; 
betaactual=100; 
t=sin(x)+randn(1,length(x))*(1/sqrt(betaactual));
TRAINx=x(1:2:length(x)); 
TRAINt=t(1: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 
u=1; 
bn=[5 6 7 8 9 10 15 20 25 30]; 
for attempt=1:1:10; 
    TPREDICT=[]; 
    mNCOL=[];
    MEAN=linspace(-pi,pi,bn(attempt));
    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';
for i=1:1:100 
    x=-pi:0.1:pi; 
    betaactual=100; 
    t=sin(x)+randn(1,length(x))*(1/sqrt(betaactual));
    TESTx=x(2:2:length(x)); 
    TESTt=t(2:2:length(x));
    %Construncting the matrix PHItest 
    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
    t_predict=PHItest*mN+randn(length(TESTx),1)*(1/ sqrt(beta));
    TPREDICT=[TPREDICT t_predict];
end
figure(1)
subplot(5,2,u) 
plot(TPREDICT)
figure(2)
subplot(5,2,u)
plot(var(TPREDICT'))
s1=[];
SN=inv(alpha*eye(size(PHI,2))+beta*PHI'*PHI);
for j=1:1:size(PHItest,1)
    s1=[s1 (1/beta)+PHItest(j,:)*SN*PHItest(j,:)'];
end
hold on 
plot(s1,'r') 
u=u+1;
end