clear
n=20
mc = 100;

xb =randn(mc,3*(n));
yb = xb(1,:) + randn(1,(3*n));

nmodels = nchoosek(1:mc,2);

clear c1 c2 c3
for ii=1:size(nmodels,1);
  ind = nmodels(ii,:);
  x = xb(ind,1:(n));
  y = yb(:,1:(n));
  b= regress(y',[x',ones(n,1)]);
  yhat = [x',ones(n,1)]*b;
  c1(ii) = corr2(yhat,y');
  
  x = xb(ind,(n+1):(2*n));
  y = yb(:,(n+1):(2*n));
  yhat = [x',ones(n,1)]*b;
  c2(ii) = corr2(yhat,y');

  x = xb(ind,(2*n+1):(3*n));
  y = yb(:,(2*n+1):(3*n));
  yhat = [x',ones(n,1)]*b;
  c3(ii) = corr2(yhat,y');
  [ii/size(nmodels,1),  nnz((c1>0.6)&(c2>0.6)),max(c2)]
end


 
ii=find(c2==max(c2))

%for ii=find((c1>0.7)&(c2>0.7)&(c3<0.5)) 
  ind = nmodels(ii,:);
  x = xb(ind,1:(n));
  y = yb(:,1:(n));
  b= regress(y',[x',ones(n,1)]);
  yhat = [x',ones(n,1)]*b;

  subplot(1,3,1)
  plot(yhat,y,'o')
  title(['data set 1 r^2 = ',num2str(corr2(yhat,y')^2,3)])
  axis('square')
  axis('tight')
  xlabel('yhat')
  ylabel('y')
  
  
  
  x = xb(ind,(n+1):(2*n));
  y = yb(:,(n+1):(2*n));
  yhat = [x',ones(n,1)]*b;

  subplot(1,3,2)
  plot(yhat,y,'o')
  title(['data set 2 r^2 = ',num2str(corr2(yhat,y')^2,3)])
  axis('square')
  axis('tight')
  xlabel('yhat')
  ylabel('y')
  
  x = xb(ind,(2*n+1):(3*n));
  y = yb(:,(2*n+1):(3*n));
  yhat = [x',ones(n,1)]*b;

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

%pause
%end

n=100000;
x =randn(2,3*(n));
y = x(1,:) + randn(1,(3*n));
yhat = [x',ones(3*n,1)]*b;
corr2(yhat,y')^2