clf
a=1;
n=50;
x = randn(1,n);
y = a*x + randn(1,n);

b= regress(y',[x',ones(n,1)]);
xhat=(min(x)-1):0.1:(max(x)+1);
yhat = [xhat',ones(length(xhat),1)]*b;

plot(x,y,'o')
xlabel('x')
ylabel('y')
hold on
h=plot(xhat,yhat,'k--');
set(h,'linewidth',1.5)
hold off

print -dpdf stat1
print -depsc stat1
pause

y = a*x.^2 + 0.2*randn(1,n);

b= regress(y',[x',ones(n,1)]);
xhat=(min(x)-0.5):0.1:(max(x)+0.5);
yhat = [xhat',ones(length(xhat),1)]*b;

plot(x,y,'o')
xlabel('x')
ylabel('y')
hold on
h=plot(xhat,[yhat,a*xhat'.^2],'k--');
set(h(1),'linewidth',1.5)
set(h(2),'color',[1,1,1]*0.5,'linestyle','-')
hold off

print -dpdf stat2
print -depsc stat2
pause
%------------------------------------------
x = randn(1,n);
y = exp(a*x + 0.2*randn(1,n));

b= regress(y',[x',ones(n,1)]);
xhat=(min(x)-0.5):0.1:(max(x)+0.5);
yhat = [xhat',ones(length(xhat),1)]*b;

b=glmfit(x,y,'poisson','log');
yhat2 = glmval(b,xhat,'log');

plot(x,y,'o')
xlabel('x')
ylabel('y')
hold on
h=plot(xhat,[yhat,yhat2],'k--');
set(h(1),'linewidth',1.5)
set(h(2),'color',[1,1,1]*0.5,'linestyle','-')
hold off
legend('data','linear regression','GLM','location','northwest')
print -dpdf stat2a
print -depsc stat2a
pause



subplot(1,2,1)
x = randn(1,n);
y = a*x + 0.1*randn(1,n);

b= regress(y',[x',ones(n,1)]);
xhat=(min(x)-1):0.1:(max(x)+1);
yhat = [xhat',ones(length(xhat),1)]*b;

plot(x,y,'o')
xlabel('x')
ylabel('y')
hold on
h=plot(xhat,yhat,'k--');
set(h,'linewidth',1.5)
hold off
axis('square')
title(['good SSE = ',num2str(sum((y'-[x',ones(n,1)]*b).^2),2)])

subplot(1,2,2)
x = randn(1,n);
y = a*x + 2*randn(1,n);

b= regress(y',[x',ones(n,1)]);
xhat=(min(x)-1):0.1:(max(x)+1);
yhat = [xhat',ones(length(xhat),1)]*b;

plot(x,y,'o')
xlabel('x')
ylabel('y')
hold on
h=plot(xhat,yhat,'k--');
set(h,'linewidth',1.5)
hold off
axis('square')
title(['bad SSE = ',num2str(sum((y'-[x',ones(n,1)]*b).^2),3)])

print -dpdf stat3
print -depsc stat3
pause

%------------------------------------
subplot(1,2,1)
x = randn(1,n);
y = a*x + 0.1*randn(1,n);
x=10*x;
y=10*y;

b= regress(y',[x',ones(n,1)]);
xhat=(min(x)-1):0.1:(max(x)+1);
yhat = [xhat',ones(length(xhat),1)]*b;

plot(x,y,'o')
xlabel('x')
ylabel('y')
hold on
h=plot(xhat,yhat,'k--');
set(h,'linewidth',1.5)
hold off
axis('square')
title(['good SSE = ',num2str(sum((y'-[x',ones(n,1)]*b).^2),2)])

subplot(1,2,2)
x = randn(1,n);
y = a*x + 2*randn(1,n);
x=x/10;
y=y/10;

b= regress(y',[x',ones(n,1)]);
xhat=(min(x)-.1):0.1:(max(x)+.1);
yhat = [xhat',ones(length(xhat),1)]*b;

plot(x,y,'o')
xlabel('x')
ylabel('y')
hold on
h=plot(xhat,yhat,'k--');
set(h,'linewidth',1.5)
hold off
axis('square')
title(['bad SSE = ',num2str(sum((y'-[x',ones(n,1)]*b).^2),3)])

print -dpdf stat4
print -depsc stat4
pause
%------------------------------------
subplot(1,2,1)
x = randn(1,n);
y = a*x + 0.1*randn(1,n);
x=10*x;
y=10*y;

b= regress(y',[x',ones(n,1)]);
xhat=(min(x)-1):0.1:(max(x)+1);
yhat = [xhat',ones(length(xhat),1)]*b;

plot(x,y,'o')
xlabel('x')
ylabel('y')
hold on
h=plot(xhat,yhat,'k--');
set(h,'linewidth',1.5)
hold off
axis('square')
title(['good r^2 = ',num2str(corr2(x,y)^2,3)])

subplot(1,2,2)
x = randn(1,n);
y = a*x + 2*randn(1,n);
x=x/10;
y=y/10;

b= regress(y',[x',ones(n,1)]);
xhat=(min(x)-.1):0.1:(max(x)+.1);
yhat = [xhat',ones(length(xhat),1)]*b;

plot(x,y,'o')
xlabel('x')
ylabel('y')
hold on
h=plot(xhat,yhat,'k--');
set(h,'linewidth',1.5)
hold off
axis('square')
title(['bad r^2 = ',num2str(corr2(x,y)^2,3)])

print -dpdf stat5
print -depsc stat5
pause

%------------------------------------

subplot(1,2,1)
x = randn(1,n);
y = 0.1*x + randn(1,n);

b= regress(y',[x',ones(n,1)]);
xhat=(min(x)-1):0.1:(max(x)+1);
yhat = [xhat',ones(length(xhat),1)]*b;

plot(x,y,'o')
xlabel('x')
ylabel('y')
hold on
h=plot(xhat,yhat,'k--');
set(h,'linewidth',1.5)
hold off
axis('square')
title(['bad r^2 = ',num2str(corr2(x,y)^2,3)])

subplot(1,2,2)
%x = randn(1,n);
%y = a*x + randn(1,n);
x(end+1)=15;y(end+1)=15;

b= regress(y',[x',ones(n+1,1)]);
xhat=(min(x)-.1):0.1:(max(x)+.1);
yhat = [xhat',ones(length(xhat),1)]*b;

h=plot(x,y,'o');
xlabel('x')
ylabel('y')
hold on
plot(x(end),y(end),'xm')
h=plot(xhat,yhat,'k--');
set(h,'linewidth',1.5)
hold off
axis('square')
title(['good? r^2 = ',num2str(corr2(x,y)^2,3)])

print -dpdf stat6
print -depsc stat6
pause
%-----------------------------------------------
m=20
n=m+1;

xb =randn(n,2*(n));
yb =randn(1,2*(n));

x = xb(:,1:(n));
y = yb(:,1:(n));

b= regress(y',[x',ones(n,1)]);
yhat = [x',ones(n,1)]*b;
y1=yhat;

subplot(1,2,1)
plot(yhat,y,'o')
axis('square')
xlabel('yhat')
ylabel('y')
title(['in-sample r^2 = ',num2str(corr2(yhat,y')^2,3)])

x = xb(:,(n+1):end);
y = yb(:,(n+1):end);
yhat = [x',ones(n,1)]*b;
y2=yhat;

subplot(1,2,2)
plot(yhat,y,'o')
axis('square')
xlabel('yhat')
ylabel('y')
title(['independent r^2 = ',num2str(corr2(yhat,y')^2,3)])

print -dpdf stat7
print -depsc stat7
pause

clf
plot([yb;[y1',y2']]','-x')
set(gca,'xlim',[0,2*n+1])
xx=get(gca,'xlim');
yy=get(gca,'ylim');
h=patch([xx(1),n+0.5,n+0.5,xx(1)],[yy([1,1]),yy([2,2])],'b');
set(h,'FaceAlpha',0.2)
h=patch([n+0.5,xx([2,2]),n+0.5],[yy([1,1]),yy([2,2])],'r');
set(h,'FaceAlpha',0.2)

print -dpdf stat7a
print -depsc stat7a
pause

%-----------------------------------------------
m=20
n=40

xb =randn(m,2*(n));
yb =randn(1,2*(n));

x = xb(:,1:(n));
y = yb(:,1:(n));

b= regress(y',[x',ones(n,1)]);
yhat = [x',ones(n,1)]*b;

subplot(1,2,1)
plot(yhat,y,'o')
axis('square')
xlabel('yhat')
ylabel('y')
title(['in-sample r^2 = ',num2str(corr2(yhat,y')^2,3)])

x = xb(:,(n+1):end);
y = yb(:,(n+1):end);
yhat = [x',ones(n,1)]*b;

subplot(1,2,2)
plot(yhat,y,'o')
axis('square')
xlabel('yhat')
ylabel('y')
title(['independent r^2 = ',num2str(corr2(yhat,y')^2,3)])

print -dpdf stat8
print -depsc stat8
pause

%-----------------------------------------------
m=20
n=80

xb =randn(m,2*(n));
yb =randn(1,2*(n));

x = xb(:,1:(n));
y = yb(:,1:(n));

b= regress(y',[x',ones(n,1)]);
yhat = [x',ones(n,1)]*b;

subplot(1,2,1)
plot(yhat,y,'o')
axis('square')
xlabel('yhat')
ylabel('y')
title(['in-sample r^2 = ',num2str(corr2(yhat,y')^2,3)])

x = xb(:,(n+1):end);
y = yb(:,(n+1):end);
yhat = [x',ones(n,1)]*b;

subplot(1,2,2)
plot(yhat,y,'o')
axis('square')
xlabel('yhat')
ylabel('y')
title(['independent r^2 = ',num2str(corr2(yhat,y')^2,3)])

print -dpdf stat9
print -depsc stat9
pause

%-----------------------------------------------
m=20
n=160

xb =randn(m,2*(n));
yb =randn(1,2*(n));

x = xb(:,1:(n));
y = yb(:,1:(n));

b= regress(y',[x',ones(n,1)]);
yhat = [x',ones(n,1)]*b;

subplot(1,2,1)
plot(yhat,y,'o')
axis('square')
xlabel('yhat')
ylabel('y')
title(['in-sample r^2 = ',num2str(corr2(yhat,y')^2,3)])

x = xb(:,(n+1):end);
y = yb(:,(n+1):end);
yhat = [x',ones(n,1)]*b;

subplot(1,2,2)
plot(yhat,y,'o')
axis('square')
xlabel('yhat')
ylabel('y')
title(['independent r^2 = ',num2str(corr2(yhat,y')^2,3)])

print -dpdf stat10
print -depsc stat10
pause

%-----------------------------------------------
m=10
n=50
mc=10000;

xb =randn(mc,2*(n));
yb =randn(1,2*(n));

c=zeros(mc,1);
for ii=1:mc
  c(ii) = corr2(xb(ii,:),yb);
end
[dum,id] = sort(-abs(c));
xb=xb(id(1:20),:);


x = xb(:,1:(n));
y = yb(:,1:(n));

b= regress(y',[x',ones(n,1)]);
yhat = [x',ones(n,1)]*b;

subplot(1,2,1)
plot(yhat,y,'o')
axis('square')
xlabel('yhat')
ylabel('y')
title(['in-sample r^2 = ',num2str(corr2(yhat,y')^2,3)])

x = xb(:,(n+1):end);
y = yb(:,(n+1):end);
yhat = [x',ones(n,1)]*b;

subplot(1,2,2)
plot(yhat,y,'o')
axis('square')
xlabel('yhat')
ylabel('y')
title(['independent r^2 = ',num2str(corr2(yhat,y')^2,3)])

print -dpdf stat11
print -depsc stat11
pause

%-----------------------------------------------
m=10
n=50

xb =randn(mc,2*(n));
yb =randn(1,2*(n));

c=zeros(mc,1);
for ii=1:mc
  c(ii) = corr2(xb(ii,1:n),yb(1:n));
end
[dum,id] = sort(-abs(c));
xb=xb(id(1:20),:);


x = xb(:,1:(n));
y = yb(:,1:(n));

b= regress(y',[x',ones(n,1)]);
yhat = [x',ones(n,1)]*b;

subplot(1,2,1)
plot(yhat,y,'o')
axis('square')
xlabel('yhat')
ylabel('y')
title(['in-sample r^2 = ',num2str(corr2(yhat,y')^2,3)])

x = xb(:,(n+1):end);
y = yb(:,(n+1):end);
yhat = [x',ones(n,1)]*b;

subplot(1,2,2)
plot(yhat,y,'o')
axis('square')
xlabel('yhat')
ylabel('y')
title(['independent r^2 = ',num2str(corr2(yhat,y')^2,3)])

print -dpdf stat12
print -depsc stat12
pause

%-----------------------------------------------
sb

%print -dpdf stat13
%print -depsc stat13
