echo on
% bvv.m - bias vs. variance demostration
% g(x) = 0 x < 0 and =1 for x > 0
% d(x) = g(x) + e,  e: N(0,0.1)
% population: 200 x sampled uniformly over [-0.5 0.5]
% test samples: [-0.5:0.01:0.5]
% D: randomly sampled from population
% F(W,x,n) = w1 + w2*x + w3*x^2 + ... + wn*x^(n-1)
%
% copyright (c) 1996-2001 by Yu Hen Hu
% Created:  9/8/96
% Modified: 8/30/2001
% modified: 12/23/2005: use polyval 

echo off
clear;clf;
npopu=200; % population size
maxtrials=10; % max # of trials

% 1. Sampling point
rand('seed',sum(100*clock));  
% population: npopu samples
ax=rand(npopu,1)-0.5;  ax=sort(ax); % ax is from small to large
% test samples: 101 samples uniformly sampled over [-.5 .5]
bx=[-0.5:0.01:0.5]'; lb=length(bx);

% evaluate g(x) over ax and over bx

disp('Enter your choice')
disp('1 - g(x) is a step function ')
disp('2 - g(x) is a 2nd order polynomial (default)')
gtype=input('Choose type of g(x): ');
if isempty(gtype), gtype=2; end
if gtype==1,
   % g(x) is a step function
   da=[ax >= 0]+randn(npopu,1)*0.1;  gb=[bx >=0];
elseif gtype==2,
   coef=[1 -2 -3]';
   da=polyval(coef,ax)+randn(npopu,1)*0.1; 
   gb = polyval(coef,bx);
end

dsize=input(['size of sample set D (< ' int2str(npopu) ...
   ', default = ' int2str(round(npopu*0.5)) ') = ']);  % finite sampling set size
if isempty(dsize), dsize=round(npopu*0.5); end
   
done=0;
while done==0,
   
npoly=input(['Enter order of polynomials (<' int2str(dsize) ', default = 3) = ']);
if isempty(npoly), npoly=3; end

for ntrial=1:maxtrials,
   % sample D: a subset of the population
   tmp=randomize([ax da]); D = tmp(1:dsize,:); % first col x, 2nd col d(x)
   wd=polyfit(D(:,1),D(:,2),npoly);    % fitting polynomial
   % wd: npoly+1 x 1,
   fwdx(:,ntrial)=polyval(wd,bx);
   % plot the result for each trial
   figure(1),plot(D(:,1),D(:,2),'g.',bx,gb,'k',bx,fwdx(:,ntrial),'b:')
   legend('d(x)','g(x)','F(W,x)'), 
   title(['trial # ' int2str(ntrial) ', polynomial order = ' int2str(npoly)]);
   drawnow
   pause
end
figure(1),plot(ax,da,'g.',bx,gb,'k',bx,fwdx,'b:')
legend('d(x)','g(x)','F(W,x)'),
title(['Ensemble of F(W,x) over ' int2str(maxtrials) ...
      ' sets of Ds. Polynomial order = ' int2str(npoly)])
axis([-.5 .5 -.5 1.5])

% calculate bias and variance

edfwx=mean(fwdx')'; % lb x 1
bias=abs(gb-edfwx).^2; maxbias=max(bias);
tmp=fwdx-edfwx*ones(1,maxtrials); % lb x maxtrials
varian=sqrt((sum(tmp'.*tmp')/maxtrials)');
maxvar=max(varian);
figure(2),subplot(211),plot(bx,bias),ylabel('bias')
axis([-.5 .5 0 maxbias])
title(['Bias and Variance, Polynomial order =' int2str(npoly) '.'])
subplot(212),plot(bx,varian),ylabel('Std')
axis([-.5 .5 0 maxvar])


done=input('To end program enter 1, to continue, press Enter: ')
if isempty(done), done=0; end
end % while loop