k=1;
VAR=[];
BIASSQUARE=[];
for lambda=linspace(0,10,25)
yestmean=0;    
for iteration=1:1:1000
x=rand(1,158)*(pi/2);
x=sort(x);
y=sin(x)+2*randn(1,length(x));
x1=linspace(0,pi/2,158);
z=sin(x1);
x=x';
M=[x x.^2 x.^3 x.^4 x.^5 x.^6 x.^7 x.^8 x.^9];    
weights=pinv(lambda*diag(ones(1,size(M,2)))+ M'*M)*M'*y';  
yest=M*weights;
figure(1)
subplot(5,5,k)
plot(z,'b')
hold on
plot(yest,'r')
YEST{iteration}=yest;
yestmean=yestmean+yest;
end
yestmean=yestmean/1000;
variance=0;
biassquare=0;
for iteration=1:1:1000
variance=variance+(YEST{iteration}-yestmean).^2;
end
variance=sum(variance)/(1000*158);
VAR=[VAR variance];
figure(2)
biassquare=(yestmean-z').^2;
plot(yestmean,'r')
hold on
plot(z,'b')
biassquare=sum(biassquare)/158;
BIASSQUARE=[BIASSQUARE  biassquare];
k=k+1;
end    
figure
plot(VAR)
title('Variance versus regularization constant')
figure
plot(BIASSQUARE)
title('Bias square versus regularization constant')


