clear
nt = 30;
mc=2000
rs = 0:0.05:1;
rms = sqrt(1-rs.^2);

rcv = zeros(mc,length(rs));
rmscv = zeros(mc,length(rs));
for ir = 1:length(rs)
  for im=1:mc
  r = rs(ir);
  x = randn(1,nt);
  y = r*x + sqrt(1-r^2)*randn(1,nt);
  X = [x',ones(nt,1)];
  %X = [ones(nt,1)];
  V = X*inv(X'*X)*X';%HAT matrix
  res = y*(eye(nt)-V);
  res_cv = res./(1-diag(V))';%cv residual w/o cv!!!!!!!!!!!!!! see book
  yp = y - res_cv;

  rcv(im,ir) = corr2(y,yp);

  rmscv(im,ir) = sqrt(mean((y-yp).^2));

  end
  
end

figure(1)

plot(rs,[rs;mean(rcv)])
axis('square')
grid
legend('','cv','location','southeast')
xlabel('population correlation')
title('cv correlation')
filename='bias_c_cv';
export_fig('-pdf',filename)


pause

plot(rs.^2,[1-rms.^2;1-mean(rmscv.^2)])
axis('square')
grid
xlabel('population correlation')
title('cv 1 - normalized error variance')
filename='bias_e_cv';
export_fig('-pdf',filename)

