clear
n=27;


rs = [eps,0.1:0.05:0.5];
mc = 10000;


for ir =1:length(rs)
  rr=rs(ir);
  f = randn(mc,n);
  fstd = std(f,[],2);
  fm = mean(f,2);
  f = (f - fm(:,ones(1,n)))./fstd(:,ones(1,n));
  o = 20+10*(f+sqrt(1-rr^2)*randn(mc,n)./rr);
    
  r = correlation(o,f);
  r(r<0.)=eps;
  sf = sqrt((1-r.^2)./r.^2);
    
  b13=norminv(1/3)*sqrt(sf.^2+1);
  b23=norminv(2/3)*sqrt(sf.^2+1);
    
  Pb = normcdf(b13(:,ones(1,n)),f,sf(:,ones(1,n)));
  Pa = 1-normcdf(b23(:,ones(1,n)),f,sf(:,ones(1,n)));

  Ob = o < norminv(1/3,mean(o(:)),std(o(:)));
  Oa = o > norminv(2/3,mean(o(:)),std(o(:)));
  rps = mean((Pa-Oa).^2+(Pb-Ob).^2,2);
  rps_c = mean((1/3-Oa).^2+(1/3-Ob).^2,2);
  rpss0(:,ir) = 1 - rps./rps_c;
  
  r=rr*ones(mc,1);
  sf = sqrt((1-r.^2)./r.^2);
    
  b13=norminv(1/3)*sqrt(sf.^2+1);
  b23=norminv(2/3)*sqrt(sf.^2+1);
    
  Pb = normcdf(b13(:,ones(1,n)),f,sf(:,ones(1,n)));
  Pa = 1-normcdf(b23(:,ones(1,n)),f,sf(:,ones(1,n)));

  rps = mean((Pa-Oa).^2+(Pb-Ob).^2,2);
  rps_c = mean((1/3-Oa).^2+(1/3-Ob).^2,2);
  rpss00(:,ir) = 1 - rps./rps_c;
  
  for ii=1:n
    ind = setdiff([1:n],ii);
    r = correlation(o(:,ind),f(:,ind));
    fm = mean(f(:,ind),2);
    fstd = std(f(:,ind),[],2);
    r(r<0.0)=eps;
    sf = sqrt((1-r.^2)./r.^2);
    b13=norminv(1/3)*sqrt(sf.^2+1);
    b23=norminv(2/3)*sqrt(sf.^2+1);
    Pb(:,ii) = normcdf(b13,(f(:,ii)-fm)./fstd,sf);
    Pa(:,ii) = 1-normcdf(b23,(f(:,ii)-fm)./fstd,sf);
  end
    
  rps = mean((Pa-Oa).^2+(Pb-Ob).^2,2);
  rpss1(:,ir) = 1 - rps./rps_c;

  for ii=1:n
    ind = setdiff([1:n],ii);
    r = correlation(o(:,ind),f(:,ind));
    fm = mean(f(:,ind),2);
    fstd = std(f(:,ind),[],2);
    r(r<0.3)=eps;
    sf = sqrt((1-r.^2)./r.^2);
    b13=norminv(1/3)*sqrt(sf.^2+1);
    b23=norminv(2/3)*sqrt(sf.^2+1);
    Pb(:,ii) = normcdf(b13,(f(:,ii)-fm)./fstd,sf);
    Pa(:,ii) = 1-normcdf(b23,(f(:,ii)-fm)./fstd,sf);
  end

  rps = mean((Pa-Oa).^2+(Pb-Ob).^2,2);
  rpss1a(:,ir) = 1 - rps./rps_c;
    
  for ii=1:n
    ind = setdiff([1:n],ii);
    r = correlation(o(:,ind),f(:,ind));
    %r = 0.8*r;
    fm = mean(f(:,ind),2);
    fstd = std(f(:,ind),[],2);
    r(r<0.0)=eps;
    sf = sqrt((1-r.^2)./r.^2)*sqrt(2.00);
    b13=norminv(1/3)*sqrt(sf.^2+1);
    b23=norminv(2/3)*sqrt(sf.^2+1);
    Pb(:,ii) = normcdf(b13,(f(:,ii)-fm)./fstd,sf);
    Pa(:,ii) = 1-normcdf(b23,(f(:,ii)-fm)./fstd,sf);
  end
  rps = mean((Pa-Oa).^2+(Pb-Ob).^2,2);
  rpss2(:,ir) = 1 - rps./rps_c;
    
  for ii=1:n
    ind = setdiff([1:n],ii);
    om = mean(o(:,ind),2);
    ostd = std(o(:,ind),[],2);
    fm = mean(f(:,ind),2);
    a = mean((f(:,ind)-fm(:,ones(1,nnz(ind)))).*(o(:,ind)-om(:,ones(1,nnz(ind)))),2)./std(f(:,ind),[],2);
    fr = a(:,ones(1,n)).*(f-fm(:,ones(1,n)))+om(:,ones(1,n));
    sf = sqrt(mean((o(:,ind)-fr(:,ind)).^2,2));
    b13=norminv(1/3,om,ostd);
    b23=norminv(2/3,om,ostd);
    Pb(:,ii) = normcdf(b13,fr(:,ii),sf);
    Pa(:,ii) = 1-normcdf(b23,fr(:,ii),sf);
  end
  rps = mean((Pa-Oa).^2+(Pb-Ob).^2,2);
  rpssR(:,ir) = 1 - rps./rps_c;
  
  
  
  plot(rs(1:ir),[mean(rpss00);mean(rpss0);mean(rpss1);mean(rpss1a);mean(rpss2);mean(rpssR)]','-x')
  grid
  title(['n=',int2str(n)])
  xlabel('correlation')
  ylabel('RPSS')
  legend('population','sample+th=0','cv-1+th=0','cv-1+th=0.3','cv-1 N/S x200%','regression residual','location','nw')
  drawnow
end
axis('tight')


  
  