clear
N=40;
nt=25;
MC = 10;
rs = 0:0.1:1;
threshold = 0.3;
C = zeros(MC,N)+NaN;
for n=1:N
  for mc=1:MC
    x=randn(n,nt);
    y = randn(1,nt);

    %compute correlation between y and pool of predictors
    c = correlation(x,y(ones(n,1),:));
    %choose the ones whose correlation exceeds the threshold
    good = find(abs(c)>threshold);
    n2 = min([nt,length(good)]);
    x2 = x(good,:);

    %PCA of retained predictors
    [u,s,v] = svd(x2,'econ');
    x2 = v';

    %keep PCs whose corelation exceeds the threshold
    c2 = correlation(x2,y(ones(n2,1),:));
    [dum,isort]=sort(-c2);
    x3 = x2(isort,:);
    [b,se,pval,inmodel,stats,nextstep,h]=...
        stepwisefit(x2',y','maxiter',8,'display','off');
    x3=x3(inmodel,:)
    %cross-validated regression
    for ii=1:nt
      ind = setdiff([1:n],ii);
      b= regress(y',[x3;ones(1,nt)]');
      xp(ii) = [x3(:,ii);1]'*b;
    end
   C(mc,n)= corr2(y,xp);
  end
  plot([1:N],[prctile(C,95)],'-x')
  legend('95% sig.','location','southeast')
  xlabel('selection pool size')
  grid
  drawnow
end