function [Cost,theta0,theta1]=linear_regression
rand('state',0); randn('state',0);
d=1; m=10;
theta = [2, .7];
X = 20*rand(m,d); Xplt = (0:.05:20)';
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)';





function cost = J(theta,X,y)
  yhat = X*theta';
  cost = (y-yhat)*(y-yhat)';

function dJ(theta,X,y)




