%Biasvardemo.m
%Bias variance Normalization demo 
W1=rand(1,5); 
B1=rand(5,1);
W2=rand(5,1);
B2=rand(1,1);
%DATA1 
INPUT=rand(4000,1)*2-1; 
H=INPUT*W1+repmat(B1',4000,1);
H=1./(1+exp(H));
OUTPUT=H*W2+repmat(B2',4000,1); 
OUTPUT=1./(1+exp(OUTPUT));
for i=1:1:40
    INPUTTRAIN{i}=INPUT((i-1)*100+1:1:(i-1)*100+100,:);
    OUTPUTTRAIN{i}=OUTPUT((i-1)*100+1:1:(i-1)*100+100,:)+ randn(100,1)*sqrt(0.2); 
end
INPUTTEST=rand(40,1)*2-1;
H=INPUTTEST*W1+repmat(B1',40,1); 
H=1./(1+exp(H));
OUTPUT=H*W2+repmat(B2',40,1); 
OUTPUTTEST=1./(1+exp(OUTPUT));
sigma=1.2;
SSECOL=[];
VAR=[];
BIAS=[];
k=0;
REG=[0 0.1:0.1:0.5 1 2 5 10 15 20:20:100];
for i=1:1:length(REG)
    k=k+1;
    for j=1:1:40
        [temp,MSE,W]=performregression(INPUTTRAIN{j}, OUTPUTTRAIN{j},INPUTTEST,OUTPUTTEST,REG(i),sigma); 
        OUTPUTVALR{j}=temp'; 
        figure(k) 
        subplot(2,1,1) 
        plot(OUTPUTVALR{j},'b')
        hold on 
    end
    subplot(2,1,2) 
    MEANESTIMATE=mean(cell2mat(OUTPUTVALR'));
    hold on 
    plot(MEANESTIMATE,'r') 
    hold on
    plot(OUTPUTTEST,'k')
    BIAS=[BIAS sum((OUTPUTTEST'-MEANESTIMATE).^2)/length(OUTPUTTEST)]; 
    VAR=[VAR mean(var(cell2mat(OUTPUTVALR')))]; 
end
figure 
subplot(2,2,1) 
stem(REG,log(BIAS),'r') 
subplot(2,2,2) 
stem(REG,log(VAR),'g')
subplot(2,2,3)
stem(REG,log(BIAS+VAR),'b')
[P1,Q1]=min(BIAS+VAR);
MSECOL=[];
for i=1:1:length(REG) 
    [OUTPUTTESTR,MSE,W]=performregression(INPUTTRAIN{1}, OUTPUTTRAIN{1},INPUTTEST,OUTPUTTEST,REG(i),sigma);
    MSECOL=[MSECOL MSE];
end
subplot(2,2,4) 
stem(REG,log(MSECOL),'k')