%regressiondemo.m
%Data generation with noise 
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:700,:);
INPUTTEST=INPUT(701:1:1000,:);
OUTPUTTRAIN=OUTPUT(1:1:700,:)+randn(700,3)*0.2; 
OUTPUTTEST=OUTPUT(701:1:1000,:);
figure(1)
colormap(copper)
subplot(3,1,1)
imagesc(INPUTTRAIN');
subplot(3,1,2) 
imagesc(OUTPUT(1:1:700,:)')
subplot(3,1,3) 
imagesc(OUTPUTTRAIN')

%Using Gaussian kernel based polynomial regression 
[V,M]=kmeans(INPUTTRAIN,4); 
[V1,I1]=min(sum((INPUTTRAIN-repmat(M(1,:), 700,1)).^2)); 
[V2,I2]=min(sum((INPUTTRAIN-repmat(M(2,:), 700,1)).^2)); 
[V3,I3]=min(sum((INPUTTRAIN-repmat(M(3,:), 700,1)).^2)); 
[V4,I4]=min(sum((INPUTTRAIN-repmat(M(3,:), 700,1)).^2)); 
sigma=mean([V1 V2 V3 V4]); 
B1=[]; B2=[]; B3=[]; B4=[]; 
for i=1:1:size(INPUTTRAIN,1)
    B1=[B1 GBF(INPUTTRAIN(i,:),M(1,:),sigma)]; 
    B2=[B2 GBF(INPUTTRAIN(i,:),M(2,:),sigma)]; 
    B3=[B3 GBF(INPUTTRAIN(i,:),M(3,:),sigma)];
    B4=[B4 GBF(INPUTTRAIN(i,:),M(4,:),sigma)];
end
%Formulation of matrix 
MAT=[ones(700,1) B1' B2' B3' B4'];
W=pinv(MAT)*OUTPUTTRAIN; 
save W W
%TESTING 
B1=[]; B2=[]; B3=[]; B4=[]; 
for i=1:1:size(INPUTTEST,1)
    B1=[B1 GBF(INPUTTEST(i,:),M(1,:),sigma)]; 
    B2=[B2 GBF(INPUTTEST(i,:),M(2,:),sigma)]; 
    B3=[B3 GBF(INPUTTEST(i,:),M(3,:),sigma)]; 
    B4=[B4 GBF(INPUTTEST(i,:),M(4,:),sigma)]; 
end
MAT=[ones(300,1) B1' B2' B3' B4'] 
OUTPUTTESTR=MAT*W;
figure(2) 
subplot(5,2,1) 
plot(INPUTTEST(:,1)) 
subplot(5,2,2) 
plot(INPUTTEST(:,2)) 
subplot(5,2,3)
plot(INPUTTEST(:,3))
subplot(5,2,4)
plot(INPUTTEST(:,4)) 
subplot(5,2,5) 
plot(INPUTTEST(:,5))
subplot(5,2,6) 
plot(INPUTTEST(:,6))
subplot(5,2,7)
plot(INPUTTEST(:,7)) 
subplot(5,2,8)
plot(INPUTTEST(:,8))
subplot(5,2,9)
plot(INPUTTEST(:,9))
subplot(5,2,10) 
plot(INPUTTEST(:,10))
figure(3)
subplot(3,1,1)
plot(OUTPUTTEST(:,1))
hold on
plot(OUTPUTTESTR(:,1))
subplot(3,1,2)
plot(OUTPUTTEST(:,2))
hold on
plot(OUTPUTTESTR(:,2))
subplot(3,1,3) 
plot(OUTPUTTEST(:,3))
hold on 
plot(OUTPUTTESTR(:,3))
figure(4)
colormap(copper) 
subplot(3,1,1)
imagesc(INPUTTEST')
subplot(3,1,2)
imagesc(OUTPUTTEST')
subplot(3,1,3) 
imagesc(OUTPUTTESTR')
%Prior probability density function of INPUT data is %uniform distribution between -1 to 1 %Likelihood function of output given INPUT data %Given the first INPUT data, output is distributed as follows.
figure(5)
colormap(copper) 
OUTPUTNOISE=[];
for i=1:1:1000 
    OUTPUTNOISE=[OUTPUTNOISE; OUTPUT(1,:)+randn(1,3)* 0.2];
end
scatter3(OUTPUTNOISE(:,1),OUTPUTNOISE(:,2), OUTPUTNOISE(:,3),'g*');
hold on 
M1=mean(OUTPUTNOISE);
scatter3(M1(1),M1(2),M1(3),'or');
hold on
scatter3(OUTPUT(1,1),OUTPUT(1,2),OUTPUT(1,3),'+r');
figure(6)
colormap(copper)
OUTPUTNOISE=[];
for i=1:1:1000
    OUTPUTNOISE=[OUTPUTNOISE; OUTPUT(2,:)+randn(1,3)* 0.2];
end
scatter3(OUTPUTNOISE(:,1),OUTPUTNOISE(:,2), OUTPUTNOISE(:,3),'g*');
M2=mean(OUTPUTNOISE); 
hold on scatter3(M2(1),M2(2),M2(3),’or’); 
hold on scatter3(OUTPUT(2,1),OUTPUT(2,2),OUTPUT(2,3),’+r’);
figure(7)
colormap(copper)
OUTPUTNOISE=[]; 
for i=1:1:1000 
    OUTPUTNOISE=[OUTPUTNOISE; OUTPUT(3,:)+randn(1,3)* 0.2]; 
end
scatter3(OUTPUTNOISE(:,1),OUTPUTNOISE(:,2), OUTPUTNOISE(:,3),'g*');
M3=mean(OUTPUTNOISE); hold on scatter3(M3(1),M3(2),M3(3),’or’);
hold on 
scatter3(OUTPUT(3,1),OUTPUT(3,2),OUTPUT(3,3),'+r');