rand('state',0); randn('state',0);

% create some data for testing
X=[ 0.0680    0.1279    0.2270    0.4919    0.5611    0.6624    0.7016    0.8422    0.9297    0.9873    0.9988]';
Y=[ 0.3260    0.2880    0.2675    0.4868    0.6155    0.6096    0.7588    0.8465    0.7939    0.8553    0.9751]';
% and more finely spaced points for visualizing the regression function
Xplt = linspace(-.25,1.25,100)'; ax=[-.25 1.25,-.25,1.25];

% generate basic linear features and fit parameters
X1=[1+0*X, X]; th1 = Y'/X1';  Xplt1=[1+0*Xplt, Xplt];

% generate polynomial features and fit parameters
Xp{1}=[1+0*X]; thp{1}=Y'/Xp{1}'; XpP{1}=[1+0*Xplt]; 
for p=2:11, Xp{p}=[Xp{p-1},X.^(p-1)]; thp{p}=Y'/Xp{p}'; XpP{p}=[XpP{p-1},Xplt.^(p-1)]; end;




% Plot a scatter plot of the data
figure(1); plot(X,Y,'r.','markersize',40); axis(ax); set(gca,'fontsize',24); axis(ax); axis square;
print('-dpng','-r72','linRegScatter.png');

% Plot error residual
figure(1); plot(Xplt,XpP{1}*thp{1}','k-',[X X]',[Y Xp{1}*thp{1}']','b-',X,Y,'r.','markersize',40,'linewidth',3); 
axis(ax); set(gca,'fontsize',24); axis(ax); axis square;
print('-dpng','-r72','linRegResid0.png');
figure(1); plot(Xplt,XpP{2}*thp{2}','k-',[X X]',[Y Xp{2}*thp{2}']','b-',X,Y,'r.','markersize',40,'linewidth',3); 
axis(ax); set(gca,'fontsize',24); axis(ax); axis square;
print('-dpng','-r72','linRegResid1.png');


% Plot cost surface in theta0,1
theta0=-1:.01:3; theta1=-1:.01:3; Jth=0*theta0'*theta1;
for i=1:length(theta0), for j=1:length(theta1), Jth(i,j)=J([theta0(i),theta1(j)],X1,Y); end; end;
figure(2); contour(theta0,theta1,Jth',logspace(log10(.05),log10(50),10));
print('-dpng','-r72','linRegCostFn.png');

myview = [
    0.2504   -0.9681         0    0.3589
    0.5414    0.1400    0.8290   -0.7552
    0.8026    0.2076   -0.5592    8.4347
         0         0         0    1.0000];
myview = [
0.9026    0.4306    0.0000   -0.6666
-0.2408    0.5048    0.8290   -0.5465
-0.3569    0.7482   -0.5593    8.7442
0         0         0    1.0000 ];



sub=1:10:length(theta1);
figure(2); imagesc(theta1,theta0,(Jth').^.3); axis('xy'); axis square; %colormap(bone(256));
print -dpng -r72 linRegCostA.png
figure(3); [C,H2]=contour(theta1,theta0,Jth',[.05 .25 1 5 25 50],'r-'); axis square;
print -dpng -r72 linRegCostD.png
figure(4); surf(theta1(sub),theta0(sub),Jth(sub,sub)',(Jth(sub,sub)').^.3); view(myview); %print -djpeg cost3dA.jpg
print -dpng -r72 linRegCostB.png
figure(5); surf(theta1(sub),theta0(sub),Jth(sub,sub)',(Jth(sub,sub)').^.3); view(myview); hold on;
[C,H3]=contour3(theta1,theta0,Jth',[.05 .25 1 5 25 50],'r-');
set(H2,'linewidth',3); set(H3,'linewidth',4); %print -djpeg cost3dB.jpg
print -dpng -r72 linRegCostC.png

% Compute gradient at each point in theta0,1
%for i=1:length(theta0), for j=1:length(theta1), dJth(i,j,:)=dJ([theta0(i),theta1(j)],X1,Y); end; end;
%figure(3); quiver(theta0(1:20:end),theta1(1:20:end),.1*dJth(1:20:end,1:20:end,1),.1*dJth(1:20:end,1:20:end,2));

% Plot gradient sequence in theta space, (x,y) space
theta=[.2 -.5]; alpha=.05; plotIters = [0 1 10 30 90];
for iter=0:100,
	if (any(plotIters==iter))
		figure(3); plot(Xplt,Xplt1*theta','k-',[X X]',[Y X1*theta']','b-',X,Y,'r.','markersize',40,'linewidth',3); axis(ax); axis square;
		print('-dpng','-r72',['linRegGradA' sprintf('%02d',iter) '.png');
		figure(4); clf; contour(theta0,theta1,Jth',logspace(log10(.05),log10(50),10)); 
			hold on; plot(theta(1),theta(2),'b.','markersize',50); axis square;
		print('-dpng','-r72',['linRegGradB' sprintf('%02d',iter) '.png');
	end;
	theta=theta+alpha*dJ(theta,X1,Y);
end;


% Plot 2D regression picture

% Plot 2D regression when features are x, x^2

% Plot various regression fits
figure(3); hold off;
for p=[1,2,4,8,11],
plot(X,Y,'r.',Xplt,XpP{p}*thp{p}','k-','markersize',40,'linewidth',3); set(gca,'fontsize',24);
axis(ax); axis square;
print('-dpng','-r72',['linRegFit' sprintf('%02d',p-1) '.png']);
pause(1); 
end;


% Validation data

Xt = [  0.0019   0.1763    0.3457    0.4103    0.6154    0.7382    0.7919    0.8936    0.9218    0.9855 ]';
%Yt = [1+0*Xt Xt Xt.^2 Xt.^3]*thp{4}' + .05*randn(10,1);
Yt = [     0.3521    0.2848    0.3853    0.3836    0.6139    0.7237    0.7884    0.9374    0.8431    0.8230 ]';
for p=[1 2 4 8 11],
  plot(Xplt,XpP{p}*thp{p}','k-',X,Y,'r.',Xt,Yt,'g.','markersize',40,'linewidth',3); set(gca,'fontsize',24); axis(ax); axis square;
  print('-dpng','-r72',['linRegVal' sprintf('%02d',p-1) '.png']);
end;

% Training/test error versus complexity
Xtp=[1+0*Xt]; val=zeros(1,11); tra=zeros(1,11);
for p=1:11,
	tra(p) = sum( (Y-Xp{p}*thp{p}').^2);
	val(p) = sum( (Yt-Xtp*thp{p}').^2 );
	Xtp = [Xtp Xt.^p];
end;
subplot(2,1,1); plot(1:11,tra,'r-',1:11,val,'g-','linewidth',3); axis([0 12 -.1 .6]); set(gca,'fontsize',20);
text(9.5,.5,'Overfitting','fontsize',16); text(.5,.5,'Underfitting','fontsize',16);
print -dpng -r72 linRegValidateFull.png
!convert -trim linRegValidateFull.png  linRegValidate.png


% generate polynomial features and fit parameters
Xp{1}=[1+0*X]; thr{1}=Y'*Xp{1}*inv(Xp{1}'*Xp{1}+1); XpP{1}=[1+0*Xplt]; 
for p=2:11, Xp{p}=[Xp{p-1},X.^(p-1)]; thr{p}=Y'*Xp{p}*inv(Xp{p}'*Xp{p}+1*eye(p)); XpP{p}=[XpP{p-1},Xplt.^(p-1)]; end;
for p=[1 2 4 8 11],
  plot(Xplt,XpP{p}*thr{p}','k-',X,Y,'r.',Xt,Yt,'g.','markersize',40,'linewidth',3); set(gca,'fontsize',24); axis(ax); axis square;
	pause(2);
  print('-dpng','-r72',['linRegReg' sprintf('%02d',p-1) '.png']);
end;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%ax = [-10,30,-10,30];

X4 = [ones(m,1),X,X.^2,X.^3];  X4plt =[ones(size(Xplt)),Xplt,Xplt.^2,Xplt.^3];
y  = X4*thetaTrue' + .1*randn(m,1);
figure; plot(X,y,'r.','markersize',40); 


figure; plot(Xplt,X4plt*thetaTrue','k-','markersize',40); 


rand(10);
X1 = [ones(m,1),X]; X1plt=[ones(size(Xplt)),Xplt];
y=X1*theta' + 3*randn(m,1);
Xt = 20*rand(m,d);
yt=[1+0*Xt,Xt]*theta' + 3*randn(m,1);
Xt2 = 20*rand(1000,d);
yt2=[1+0*Xt2,Xt2]*theta' + 3*randn(1000,1);

figure(1); H1=plot(X(:,1),y,'r.','markersize',40); set(gca,'fontsize',14);
print -dpng -r40 linRegScatter.png

theta0=-40:.5:40; theta1=-1:.05:3;
Cost=zeros(length(theta0),length(theta1));
for i=1:length(theta0),
  for j=1:length(theta1),
    Cost(i,j)=mean((y-X1*[theta0(i),theta1(j)]').^2);
%    Cost(i,j) = J([theta0(i),theta1(j)],X1,y);
end; end;
myview = [
    0.2504   -0.9681         0    0.3589
    0.5414    0.1400    0.8290   -0.7552
    0.8026    0.2076   -0.5592    8.4347
         0         0         0    1.0000];
figure(2); imagesc(theta1,theta0,sqrt(Cost)); axis('xy'); %colormap(bone(256));
print -dpng -r40 linRegCostA.png
figure(3); [C,H2]=contour(theta1,theta0,Cost,[10,30,80,320,1280],'r-');
print -dpng -r40 linRegCostD.png
figure(4); surf(theta1,theta0,Cost,sqrt(Cost)); view(myview); %print -djpeg cost3dA.jpg
print -dpng -r40 linRegCostB.png
figure(5); surf(theta1,theta0,Cost,sqrt(Cost)); view(myview); hold on;
[C,H3]=contour3(theta1,theta0,Cost,[10,30,80,320,1280],'r-');
set(H2,'linewidth',3); set(H3,'linewidth',4); %print -djpeg cost3dB.jpg
print -dpng -r40 linRegCostC.png

thetah{1}=[-20,0]; iter=1;
figure(6); clf; [C,H2]=contour(theta1,theta0,Cost,[10,30,80,320,1280],'r-'); set(H2,'linewidth',3); 
hold on;  plot(thetah{1}(2),thetah{1}(1),'b.','markersize',40);
print('-dpng','-r40','linRegStep1a.png');
figure(7); clf; H1=plot(X(:,1),y,'r.','markersize',40); hold on; plot(Xplt,X1plt*thetah{iter}','k-','linewidth',4);
print('-dpng','-r40','linRegStep1b.png');
alpha=[3e-2,3e-4]; %pause;
for iter=2:50,
  e=(y - X1*thetah{iter-1}')'*X1;% e = repmat(e,[1,d+1]); e = e.*X1; e=sum(e,1);
  thetah{iter} = thetah{iter-1} + alpha.*e;
%  plot([thetah{iter-1}(1),thetah{iter}(1)],[thetah{iter-1}(2),thetah{iter}(2)],'r-');
  if (log2(iter)==fix(log2(iter)))
		th0=[thetah{iter/2}(2) thetah{iter/2}(1)];
		th1=[thetah{iter}(2) thetah{iter}(1)];
    figure(6); 
		clf; [C,H2]=contour(theta1,theta0,Cost,[10,30,80,320,1280],'r-'); set(H2,'linewidth',3); hold on;
		plot([th0(1) th1(1)],[th0(2) th1(2)], 'b.','markersize',40);
    %figure(6); plot(thetah{iter}(2),thetah{iter}(1),'b.','markersize',40);
		hold on; arrow([th0 2],[th0+.5*(th1-th0) 2],'EdgeColor','b','FaceColor','b'); hold off;
		if (iter==2) print('-dpng','-r40','linRegStep2a.png'); end;
		if (iter==4) print('-dpng','-r40','linRegStep4a.png'); end;
		if (iter==8) print('-dpng','-r40','linRegStep8a.png'); end;

    figure(7); %clf; H1=plot(X(:,1),y,'r.','markersize',40); hold on; 
    hold off; plot(X(:,1),y,'r.',Xplt,X1plt*thetah{iter}','k-','markersize',40,'linewidth',4);
		if (iter==2) print('-dpng','-r40','linRegStep2b.png'); end;
		if (iter==4) print('-dpng','-r40','linRegStep4b.png'); end;
		if (iter==8) print('-dpng','-r40','linRegStep8b.png'); end;
    %pause;
  end;
end;

figure(10); clf; ax=[-5,25,-5,25]
H1=plot(X(:,1),y,'r.','markersize',40); hold on; %ax=axis;
axis(ax); pause;
Xp=ones(size(X)); Xpplt=ones(size(Xplt));
for p=1:10,
  Xp = [Xp,X.^p]; Xpplt=[Xpplt,Xplt.^p];
  thetap{p} = (inv(Xp'*Xp)*Xp'*y)';
  clf; plot(X(:,1),y,'r.',Xt(:,1),yt,'g.',Xplt,Xpplt*thetap{p}','b-','linewidth',3,'markersize',40); %axis(ax); %axis off; 
  title(sprintf('Order %d polynomial',p));
  pause;
end;

Xp=[]; Xt2p=[];
for p=0:10,
  Xp = [Xp,X.^p]; Xt2p = [Xt2p,Xt2.^p];
  thetap{p+1} = (inv(Xp'*Xp)*Xp'*y)';
  err1(p+1)= mean((y-Xp*thetap{p+1}').^2);
  err2(p+1)= mean((yt2-Xt2p*thetap{p+1}').^2);  
end;
figure(11); plot(0:10,err1,'r-',0:10,err2,'g-','linewidth',3);
axis([0,5,0,30]);

H1=plot(X(:,1),y,'r.','markersize',40); hold on; %ax=axis;
axis(ax); pause;
Xp=ones(size(X)); Xpplt=ones(size(Xplt));
for p=1:10,
  Xp = [Xp,X.^p]; Xpplt=[Xpplt,Xplt.^p];
  thetap{p} = (inv(Xp'*Xp)*Xp'*y)';
  clf; plot(X(:,1),y,'r.',Xt(:,1),yt,'g.',Xplt,Xpplt*thetap{p}','b-','linewidth',3,'markersize',40); %axis(ax); %axis off; 
  title(sprintf('Order %d polynomial',p));
  pause;
end;


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%5
%Outliers and robust cost functions
X1 = [ones(m,1),X]; X1plt=[ones(size(Xplt)),Xplt];
th = y' / X1'; figure(20); clf; plot(X,y,'r.',Xplt,th*X1plt','k-','markersize',40,'linewidth',3);
Xout = [X;20]; yout=[y;0]; X1out = [ones(m+1,1),Xout];
th = yout' / X1out'; figure(20); clf; plot(Xout,yout,'r.',Xplt,th*X1plt','k-','markersize',40,'linewidth',3);

th = y' / X1'; err=(yout-th*X1out'); hist(err);


thetah{1}=th; iter=1;
figure(7); clf; H1=plot(Xout,yout,'r.','markersize',40); hold on; plot(Xplt,X1plt*thetah{iter}','k-','linewidth',4);
alpha=[3e-2,3e-4]; pause;
for iter=2:500,
  e=(y - X1*thetah{iter-1}'); e = repmat(e,[1,d+1]);
  e = sign(e).*X1; e=sum(e,1);
  thetah{iter} = thetah{iter-1} + alpha.*e;
end;
figure(7); %clf; H1=plot(X(:,1),y,'r.','markersize',40); hold on; 
plot(Xplt,X1plt*thetah{iter}','k-','linewidth',4);

thetah{1}=th; iter=1;
figure(7); clf; H1=plot(Xout,yout,'r.','markersize',40); hold on; plot(Xplt,X1plt*thetah{iter}','k-','linewidth',4);
alpha=[3e-2,3e-4]; pause;
for iter=2:500,
  e=(yout - X1out*thetah{iter-1}'); e = repmat(e,[1,d+1]);
  e = sign(e).*X1out; e=sum(e,1);
  thetah{iter} = thetah{iter-1} + alpha.*e;
end;
figure(7); %clf; H1=plot(X(:,1),y,'r.','markersize',40); hold on; 
plot(Xplt,X1plt*thetah{iter}','k-','linewidth',4);



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%5
%Bias and variance
d=1; m=200;
theta = [1.7681   -0.8526    0.2735   -0.0105];
rand('state',0); randn('state',0);
X = 20*rand(m,d); Xplt = (0:.05:20)';
X1 = [ones(m,1),X,X.^2,X.^3]; X1plt=[ones(size(Xplt)),Xplt,Xplt.^2,Xplt.^3];
y=X1*theta' + 3*randn(m,1);
figure(1); H1=plot(X(:,1),y,'r.',Xplt,X1plt*theta','k-','markersize',40,'linewidth',3);

Xp=[1+0*X,X,X.^2,X.^3];
Xpp=[1+0*Xplt,Xplt,Xplt.^2,Xplt.^3];
thetahat = (inv(Xp'*Xp)*Xp'*y)';



%Bias and variance
close all;
d=1; m=6; 
theta = [1.7681   0.2735 ];
rand('state',2); randn('state',2);
X = 20*rand(m,d); Xplt = (0:.05:20)';
x2 = 20*rand(m,d); x3 = 20*rand(m,d); 
X1 = [ones(m,1),X]; X1plt=[ones(size(Xplt)),Xplt];
X2 = [ones(m,1),x2]; X3 = [ones(m,1),x3]; 
y=X1*theta' + 1*randn(m,1);
y2=X2*theta' + 1*randn(m,1);
y3=X3*theta' + 1*randn(m,1);
ax = [0 20 -2 12];
%figure(1); H1=plot(X(:,1),y,'r.',Xplt,X1plt*theta','k-','markersize',40,'linewidth',3);
figure(1); H1=plot(X(:,1),y,'r.','markersize',40);   axis(ax);axis square; print('-dpng','-r40','linRegVarianceA1.png');
figure(2); H1=plot(x2(:,1),y2,'g.','markersize',40); axis(ax);axis square; print('-dpng','-r40','linRegVarianceA2.png');
figure(3); H1=plot(x3(:,1),y3,'b.','markersize',40); axis(ax);axis square; print('-dpng','-r40','linRegVarianceA3.png');

X1 = [ones(m,1)]; X2 = [ones(m,1)]; X3 = [ones(m,1)]; 
X1plt=[ones(size(Xplt))];
th1 = (X1\y)'; th2=(X2\y2)'; th3=(X3\y3)';
figure(4); H1=plot(X(:,1),y,'r.',x2(:,1),y2,'g.',x3(:,1),y3,'b.','markersize',40,'linewidth',3);
figure(4); hold on; plot(Xplt,X1plt*th1','r-',Xplt,X1plt*th2','g-',Xplt,X1plt*th3','b-','linewidth',3);
axis(ax);axis square; title('Poly Order 0');
print('-dpng','-r40','linRegVarianceB1.png');

X1 = [ones(m,1),X]; X2 = [ones(m,1),x2]; X3 = [ones(m,1),x3]; 
X1plt=[ones(size(Xplt)),Xplt];
th1 = (X1\y)'; th2=(X2\y2)'; th3=(X3\y3)';
figure(5); H1=plot(X(:,1),y,'r.',x2(:,1),y2,'g.',x3(:,1),y3,'b.','markersize',40,'linewidth',3);
figure(5); hold on; plot(Xplt,X1plt*th1','r-',Xplt,X1plt*th2','g-',Xplt,X1plt*th3','b-','linewidth',3);
axis(ax);axis square; title('Poly Order 1');
print('-dpng','-r40','linRegVarianceB2.png');

X1 = [ones(m,1),X,X.^2,X.^3]; X2 = [ones(m,1),x2,x2.^2,x2.^3]; X3 = [ones(m,1),x3,x3.^2,x3.^3]; 
X1plt=[ones(size(Xplt)),Xplt,Xplt.^2,Xplt.^3];
th1 = (X1\y)'; th2=(X2\y2)'; th3=(X3\y3)';
figure(6); H1=plot(X(:,1),y,'r.',x2(:,1),y2,'g.',x3(:,1),y3,'b.','markersize',40,'linewidth',3);
figure(6); hold on; plot(Xplt,X1plt*th1','r-',Xplt,X1plt*th2','g-',Xplt,X1plt*th3','b-','linewidth',3);
axis(ax);axis square; title('Poly Order 3');
print('-dpng','-r40','linRegVarianceB3.png');

close all;

Xp=[1+0*X,X,X.^2,X.^3];

