%gaussianprocess.m
%Gaussian process 
W1=rand(10,5);
B1=rand(5,1); 
W2=rand(5,3);
B2=rand(3,1);
INPUT=rand(1000,10)*2-1;
H=INPUT*W1+repmat(B1',1000,1);
H=1./(1+exp(H)) 
OUTPUT=H*W2+repmat(B2',1000,1);
OUTPUT=1./(1+exp(OUTPUT));
INPUTTRAIN=INPUT(1:1:400,:);
VALIDTRAIN=INPUT(401:1:700,:)
INPUTTEST=INPUT(701:1:1000,:);
OUTPUTTRAIN=OUTPUT(1:1:400,:)+randn(400,3)*0.2;
OUTPUTVALID=OUTPUT(401:1:700,:)+randn(300,3)*0.2;
OUTPUTTEST=OUTPUT(701:1:1000,:);
SSECOL=[]; 
var=[0.1:0.5:10 10:10:100 100:100:2000];
for iter=1:1:50
%Formulation of kernel matrix 
for i=1:1:length(INPUTTRAIN) 
    for j=1:1:length(INPUTTRAIN) 
        K(i,j)=gbf(INPUTTRAIN(i,:),INPUTTRAIN(j,:), var(iter)); 
    end
end
%Formulation of C matrix 
C=K+0.2*eye(400); %Prediction is obtained as follows 
TARGET=[]; 
for i=1:1:length(VALIDTRAIN) 
    temp=[]; 
    for j=1:1:length(INPUTTRAIN) 
        temp=[temp gbf(VALIDTRAIN(i,:),INPUTTRAIN(j,:), var(iter))]; 
    end
    TARGET=[TARGET (temp*inv(C)*OUTPUTTRAIN)']; 
end
SSE=sum(sum((TARGET-OUTPUTVALID').^2));
SSECOL=[SSECOL SSE];
end
[P,Q]=min(SSECOL);
TARGET=[];
for i=1:1:length(INPUTTEST) 
    temp=[]; 
    for j=1:1:length(INPUTTRAIN)
        temp=[temp gbf(INPUTTEST(i,:),INPUTTRAIN(j,:),Q(1))]; 
    end
    TARGET=[TARGET (temp*inv(C)*OUTPUTTRAIN)'];
end
figure(1)
subplot(3,1,1)
plot(TARGET(1,:))
hold on 
plot(OUTPUTTEST(:,1),'r')
subplot(3,1,2)
plot(TARGET(2,:))
hold on 
plot(OUTPUTTEST(:,2),'r')
subplot(3,1,3)
plot(TARGET(3,:)) 
hold on 
plot(OUTPUTTEST(:,3),'r')