function output=experiment3()

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%   date:2006-5-23
%%%   content:"bias vers variance relationship"
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clear all
clc

n=100;
m=100;

X=zeros(n,m);
Y=zeros(n,m);
TY=zeros(n,m);

E=zeros(n,4);
bias=zeros(n,4);
var=zeros(n,4);

%%gfr = random('Normal',0,0.1,[2,4]);%%mu=0;sigma=0.1
%%unir = unifrnd(-1,1,[1,10]);%%Random numbers from the continuous uniform distribution
X=unifrnd(-1,1,[n,m]);
Y=X.^2+random('Normal',0,0.1,[n,m]);
TY=X.^2;

for i=1:1:n
x=X(i,:);y=Y(i,:);ty=TY(i,:);
x=x';y=y';ty=ty';

[E(i,1),bias(i,1),var(i,1)]=fitfun1(x,y,ty);
[E(i,2),bias(i,2),var(i,2)]=fitfun2(x,y,ty);
[E(i,3),bias(i,3),var(i,3)]=fitfun3(x,y,ty);
[E(i,4),bias(i,4),var(i,4)]=fitfun4(x,y,ty);
end

for t=1:1:4
    EE=E(:,t);    
    minE=min(EE);MaxE=max(EE);
    temp=minE:(MaxE-minE)*0.1:MaxE;
    histE=zeros(size(temp));
    for i=2:1:length(temp)
        index=find(E>=temp(i-1)&E<temp(i));
        histE(i)=length(index);
    end
    subplot(3,2,t);plot(histE);xlabel(['g' num2str(t) '(x)']);ylabel('MSE');
end
s=1:1:4;
averbias=zeros(1,4);
avervar=zeros(1,4);

averbias(s)=mean(E(:,s));%mean(bias(:,s),1);
avervar(s)=std(var(:,s),0,1);

subplot(3,2,5);h=stem(averbias);xlabel('g1(x)--g4(x)');ylabel('average bias');set(h,'MarkerFaceColor','blue');
subplot(3,2,6);h=stem(avervar);xlabel('g1(x)--g4(x)');ylabel('average variance');set(h,'MarkerFaceColor','red','Marker','square');

ss=0;
%%%%%%%%%%%%%%%%%%%%%%%
function [err,bias,var]=fitfun1(x,y,ty)
z=0.5.*ones(length(y),1);

err=mean((z-ty).^2);
%bias=(mean(z-ty))^2;
var=mean((z-mean(z)).^2);
bias=err-var;
%err=bias+var;
%%%%%%%%%%%%%%%%%%%%%%%%%

function [err,bias,var]=fitfun2(x,y,ty)

z=1.0.*ones(length(y),1);

err=mean((z-ty).^2);
%bias=(mean(z-ty))^2;
var=mean((z-mean(z)).^2);
bias=err-var;
%err=bias+var;
%%%%%%%%%%%%%%%%%%%%%%%%%%

function [err,bias,var]=fitfun3(x,y,ty)
A=zeros(length(x),2);
A(:,1)=1;
A(:,2)=x;
p=A\y;
z1=A*p;

err=mean((z1-ty).^2,1);
%bias=(mean(z1-ty))^2;
var=mean((z1-mean(z1)).^2,1);
bias=err-var;
%err=bias+var;
%%%%%%%%%%%%%%%%%%%%%%%

function [err,bias,var]=fitfun4(x,y,ty)
B=zeros(length(x),4);
B(:,1)=1;
B(:,2)=x;
B(:,3)=x.^2;
B(:,4)=x.^3;
p=B\y;
z2=B*p;

err=mean((z2-ty).^2,1);
%bias=(mean(z2-ty))^2;
var=mean((z2-mean(z2)).^2);
bias=err-var;
%err=bias+var;
%%%%%%%%%%%%%%%%%%%%%%%