
disp(' ');
disp('Pearson (rho), Kendall tau, and Spearman (R) correlation coefficients:');
disp('See Section 8.2 in text for Kendall tau and Section 8.5 for Spearman R');

disp(' ');
disp('NOTE: Bootstrap methods do quite well for Kendall`s tau (see below');
disp('  and Section 8.4), but worse for Pearson`s rho and Spearman`s R');
disp('  Perhaps Efron`s Bias-Correct Bootstrap (Section 8.4) will do better?');

disp(' ');
disp('Data: Amount of pre-packed engine lubricant and engine survival times:');
disp('(Source of data unknown.)');
disp(' ');
disp('Lubricant  Survival          Lubricant Survival');

ww2 = [ 1.   29     23     11.   84    429;
        2.   21     36     12.   37    489;
        3.   30     67     13.   94    504;
        4.   45    104     14.   56    925;
        5.   48    147     15.   67    980;
        6.   92    164     16.   79   1254;
        7.   74    247     17.   70   2961;
        8.   61    304     18.   93   5351;
        9.   99    355     19.   76   9915;
       10.   24    408     20.   72  11301  ];

fprintf (' %3d.  %2d  %5d          %3d.   %2d  %5d\n', ww2');  %'

disp(' ');
disp('Initializing random-number generator at');

Startseed=123456;
disp(Startseed);
rand('state',Startseed);

disp('Number of permutations and/or bootstrap replications:');
nsims=10000;
disp(nsims);

% Rewrite data as one nnx3 array
ww = [ ww2(:,1) ww2(:,2) ww2(:,3);  ww2(:,4) ww2(:,5) ww2(:,6)];
nn=size(ww,1);      % Number of rows
ww=sortrows(ww,1);  % Sort on ordinal number

disp('Data sorted into one nnx3 matrix:');
fprintf(' %3d   %3d  %5d\n',ww');  %'

% Fetch Xs and Ys as nnx1 column vectors
%   and wdat as [X Y];
xx=ww(:,2);  yy=ww(:,3);
wdat = [xx yy];

disp(' ');
disp('Pearson rho for X vs Y and Student-t P-value:');
Meanx=sum(xx)/nn;  Meany=sum(yy)/nn;
xxc=xx-Meanx;  yyc=yy-Meany;
sxx=sum(xxc.*xxc);  sxy=sum(xxc.*yyc);  syy=sum(yyc.*yyc);

rho = sxy/sqrt(sxx*syy);
tt = rho*sqrt((nn-2)/(1-rho*rho));
df=nn-2;
Tpval = pvtdist(tt,df);

fprintf ('  Rho=%g   T=%g   P=%6.4f  (df=%d)\n', rho, tt, Tpval, df);

disp(' ');
disp('Permutation test for P-value of Pearson rho:');
disp('(Permuting values of Y for fixed X.)');
Nge=0;  absrho=abs(rho);
for ns=1:nsims;
  permyy = yy(randperm(nn));   % Randomly permute Y-values
  Meany=sum(permyy)/nn;  yyc=permyy-Meany;
  sxx=sum(xxc.*xxc);  sxy=sum(xxc.*yyc);  syy=sum(yyc.*yyc);
  rand_rho = sxy/sqrt(sxx*syy);
  if abs(rand_rho) >= absrho  Nge=Nge+1;  end;
  end;

disp('Number greater-or-equal (for absolute value) and number of permutations:');
fprintf ('    %d    %d\n', Nge,nsims);
phat=Nge/nsims;  pwid=1.960*sqrt(phat*(1-phat)/nsims);
plow=phat-pwid;  phigh=phat+pwid;
disp('95% Conf.Int. for true permutation-P-value bracketed by estimate:');
fprintf ('    %6.4f    %6.4f    %6.4f    (nsims=%d)\n', plow,phat,phigh,nsims);

disp(' ');
disp('Fisher z-transform 95% Conf.Int. for true rho bracketing rho:');
z0=0.50*log((1+rho)/(1-rho));
zwid=1.960/sqrt(nn-3);
z1=z0-zwid;  rh1=(exp(2*z1)-1)/(exp(2*z1)+1);
z2=z0+zwid;  rh2=(exp(2*z2)-1)/(exp(2*z2)+1);
fprintf('    %6.4f    %6.4f    %6.4f\n', rh1, rho, rh2);

disp('Bootstrap 95% Percentile Conf.Int. for rho bracketing Median:');
rboot=zeros(1,nsims);
bias=0;  nle=0;  % For bootstrap bias and bootstrap one-sided P-value
for ns=1:nsims;
  bkk = ceil(nn*rand(nn,1));    % Bootstrap from [1:nn]'  '
  bxx=xx(bkk);  byy=yy(bkk);    % Xs and Ys from bootstrapped pairs
  Meanbx=sum(bxx)/nn;  Meanby=sum(byy)/nn;
  xxbc=bxx-Meanbx;  yybc=byy-Meanby;
  sxxb=sum(xxbc.*xxbc);  sxyb=sum(xxbc.*yybc);  syyb=sum(yybc.*yybc);
  brho = sxyb/sqrt(sxxb*syyb);
  rboot(ns) = brho;
  if brho<=rho bias=bias+1;  end;
  if brho<=0  nle=nle+1;  end;     % Assuming observed rho>=0
  end;
% Sort rboot and find quantiles
rboot=sort(rboot);
bmed=median(rboot);
% This assumes that nsims >= 40 !!
blim1=rboot(floor(nsims*0.025));
blim2=rboot(ceil(nsims*0.975));
fprintf('    %6.4f    %6.4f    %6.4f    Bias: %6.4f\n', blim1,bmed,blim2, bias/nsims);
% This assumes that rho >= 40 !!
fprintf('Bootstrap P-value for H_0:rho=0:  P=%6.4f (2-sided)\n', 2*nle/nsims);
 

% X vs log(Y)?
lyy=log(yy);
disp(' ');
disp('X vs log(Y):');
wwll = [ [1:nn]'  xx  lyy];  %'
fprintf(' %3d   %3d    %5.3f\n',wwll');  %'

disp(' ');
disp('Pearson rho for X vs log(Y) and Student-t P-value:');
Meanx=sum(xx)/nn;  Meany=sum(lyy)/nn;
xxc=xx-Meanx;  yyc=lyy-Meany;
sxx=sum(xxc.*xxc);  sxy=sum(xxc.*yyc);  syy=sum(yyc.*yyc);

rho = sxy/sqrt(sxx*syy);
tt = rho*sqrt((nn-2)/(1-rho*rho));
df=nn-2;
Tpval = pvtdist(tt,df);

fprintf ('  Rho=%g   T=%g   P=%6.4f  (df=%d)\n', rho, tt, Tpval, df);

disp(' ');
disp('Permutation test for P-value of Pearson rho:');
disp('(Permuting values of log(Y) for fixed X.)');
Nge=0;  absrho=abs(rho);
for ns=1:nsims;
  permyy = lyy(randperm(nn));   % Randomly permute Y-values
  Meanx=sum(xx)/nn;  Meany=sum(permyy)/nn;
  xxc=xx-Meanx;  yyc=permyy-Meany;
  sxx=sum(xxc.*xxc);  sxy=sum(xxc.*yyc);  syy=sum(yyc.*yyc);
  rand_rho = sxy/sqrt(sxx*syy);
  if abs(rand_rho) >= absrho  Nge=Nge+1;  end;
  end;
disp('Number greater-or-equal (absolute value) and number of permutations:');
fprintf ('    %d    %d\n', Nge,nsims);
phat=Nge/nsims;  pwid=1.960*sqrt(phat*(1-phat)/nsims);
plow=phat-pwid;  phigh=phat+pwid;
disp('95% Conf.Int. for true permutation-P-value bracketed by estimate:');
fprintf ('    %6.4f    %6.4f    %6.4f    (nsims=%d)\n', plow,phat,phigh,nsims);

disp(' ');
disp('Fisher z-transform 95% Conf.Int. for true rho bracketing rho:');
z0=0.50*log((1+rho)/(1-rho));
zwid=1.960/sqrt(nn-3);
z1=z0-zwid;  rh1=(exp(2*z1)-1)/(exp(2*z1)+1);
z2=z0+zwid;  rh2=(exp(2*z2)-1)/(exp(2*z2)+1);
fprintf('    %6.4f    %6.4f    %6.4f\n', rh1, rho, rh2);

disp('Bootstrap 95% Percentile Conf.Int. for rho bracketing Median:');
rboot=zeros(1,nsims);
bias=0;  nle=0;  % For bootstrap bias and bootstrap one-sided P-value
for ns=1:nsims;
  bkk = ceil(nn*rand(nn,1));     % Bootstrap from [1:nn]'  '
  bxx=xx(bkk);  byy=lyy(bkk);    % Xs and log(Y)s from bootstrapped pairs
  Meanbx=sum(bxx)/nn;  Meanby=sum(byy)/nn;
  xxbc=bxx-Meanbx;  yybc=byy-Meanby;
  sxxb=sum(xxbc.*xxbc);  sxyb=sum(xxbc.*yybc);  syyb=sum(yybc.*yybc);
  brho = sxyb/sqrt(sxxb*syyb);
  rboot(ns) = brho;
  if brho<=rho bias=bias+1;  end;
  if brho<=0  nle=nle+1;  end;     % Assuming observed rho>=0
  end;
% Sort rboot and find quantiles
rboot=sort(rboot);
bmed=median(rboot);
% This assumes that nsims >= 40 !!
blim1=rboot(floor(nsims*0.025));
blim2=rboot(ceil(nsims*0.975));
fprintf('    %6.4f    %6.4f    %6.4f    Bias: %6.4f\n', blim1,bmed,blim2, bias/nsims);
% This assumes that rho >= 40 !!
fprintf('Bootstrap P-value for H_0:rho=0:  P=%6.4f (2-sided)\n', 2*nle/nsims);
 

disp(' ');
disp(' ');
disp('Kendall K and Kendall tau (assuming no ties):');
kstat=0;
for i=1:nn;  for j=i+1:nn;
  qq = (xx(j)-xx(i))*(yy(j)-yy(i));
  if (qq>0) kstat=kstat+1;
  elseif (qq<0) kstat=kstat-1;  end;
  end;  end;
tau=(2*kstat)/(nn*(nn-1));
zval = kstat/sqrt(nn*(nn-1)*(2*nn+5)/18);
Kpval = normpval(zval);
fprintf('  K=%g   tau=%g   Z=%g   P=%6.4f (2-sided)\n', kstat,tau,zval,Kpval);

disp(' ');
disp('Permutation test for P-value of (absolute value of) Kendall tau:');
Nge=0;  abs_kstat=abs(kstat);
for ns=1:nsims;
  permyy = yy(randperm(nn));   % Randomly permute Y values
  randkk=0;
  for i=1:nn;  for j=i+1:nn;
    qq = (xx(j)-xx(i))*(permyy(j)-permyy(i));
    if (qq>0) randkk=randkk+1;
    elseif (qq<0) randkk=randkk-1;  end;
  end;  end;
  if abs(randkk) >= abs_kstat  Nge=Nge+1;  end;
  end;

disp(' ');
disp('Nge and number of permutations:');
fprintf('    Nge=%d    nsims=%d\n', Nge,nsims);

disp('95% Conf.Int. for true two-sided P-value bracketed by estimate:');
phat = Nge/nsims;
pwid = 1.960*sqrt(phat*(1-phat)/nsims);
disp([phat-pwid  phat  phat+pwid]);

disp('Asymptotically distribution-free 95% conf.int. for tau (Section 8.3):');
cc=zeros(1,nn);
for i=1:nn;  csum=0;
  for j=1:nn;
    qq = (xx(j)-xx(i))*(yy(j)-yy(i));
    if (qq>0) csum=csum+1;
    elseif (qq<0) csum=csum-1;  end;
  end;  cc(i)=csum;  end;
cmean=sum(cc)/nn;  ccx=cc-cmean;
c2sum=sum(ccx.*ccx);
cvarin = (2*(nn-2)/((nn-1)*nn*(nn+1)))*c2sum + 1-tau*tau;
cvar = 2*cvarin/(nn*(nn-1));
cwid = 1.960*sqrt(cvar);
disp([tau-cwid  tau  tau+cwid]);

disp('Bootstrap 95% Percentile Conf.Int. for Kendall tau bracketing Median:');
rboot=zeros(1,nsims);  bias=0;
bias=0;  nle=0;  % For bootstrap bias and bootstrap one-sided P-value
for ns=1:nsims;
  bww = wdat(ceil(nn*rand(1,nn)),:);   % Bootstrap rows of wdat
  bxx=bww(:,1);  byy=bww(:,2);         % X and Y columns of bootstrapped wdat
  randkk=0;
  for i=1:nn;  for j=i+1:nn;
    qq = (bxx(j)-bxx(i))*(byy(j)-byy(i));
    if (qq>0) randkk=randkk+1;
    elseif (qq<0) randkk=randkk-1;  end;
  end;  end;
  rboot(ns)=2*randkk/(nn*(nn-1));
  if randkk <= kstat  bias=bias+1;  end;
  if randkk<=0  nle=nle+1;  end;
  end;
rboot=sort(rboot);
bmed=median(rboot);
blim1=rboot(floor(nsims*0.025));
blim2=rboot(ceil(nsims*0.975));
fprintf('    %6.4f    %6.4f    %6.4f    Bias: %6.4f\n', blim1,bmed,blim2, bias/nsims);
fprintf('Bootstrap P-value for H_0:tau=0:  P=%6.4f (2-sided)\n', 2*nle/nsims);

 
disp(' ');
disp(' ');
disp('Within-sample midranks for Spearman R:');
disp('Lubricant, Lubricant midranks, Survival, Survival midranks, Tiegroups:');

[xmidr  xtiegr] = getmidranks(xx);
[ymidr  ytiegr] = getmidranks(yy);
ww = [ww(:,1)  xx  xmidr  yy  ymidr  xtiegr  ytiegr];

fprintf(' %3d    %3d  %4.1f    %5d  %4.1f    %2d  %2d\n', ww');  %'

% Center midranks at average rank
rxc=xmidr - (nn+1)/2;  ryc = ymidr - (nn+1)/2;
% Spearman R correlation coefficient
rr = 12*sum(rxc.*ryc)/((nn-1)*nn*(nn+1));
% Large-sample normal approximation
zz = sqrt(nn-1)*rr;
Rpval = normpval(zz);

disp(' ');
disp(' ');
disp('Spearman R, large-sample Z, large-sample P-value:');
fprintf('    %6.4f    %6.4f    %6.4f (2-sided)\n', rr, zz, Rpval);

disp(' ');
disp('Permutation test for P-value of (absolute value of) Spearman R:');
Nge=0;  AbsRscore = abs(sum(rxc.*ryc));  % Spearman R without normalizaton
for ns=1:nsims;
  permryc = ryc(randperm(nn));   % Randomly permute CENTERED Y midranks
  AbsRandScore = abs(sum(rxc.*permryc));
  if AbsRandScore >= AbsRscore  Nge=Nge+1;  end;
  end;

disp(' ');
disp('Nge and number of permutations:');
fprintf('    Nge=%d    nsims=%d\n', Nge,nsims);

disp('95% Conf.Int. for true two-sided P-value bracketed by estimate:');
phat = Nge/nsims;
pwid = 1.960*sqrt(phat*(1-phat)/nsims);
disp([phat-pwid  phat  phat+pwid]);

disp('Large-sample 95% Conf.Int. for R bracketing R (Section 8.5):');
pwid = 1.960/sqrt(nn-1);
disp([ rr-pwid  rr  rr+pwid ]);

disp('Bootstrap 95% Percentile Conf.Int. for R bracketing Median:');
rboot=zeros(1,nsims);
bias=0;  nle=0;  % For bootstrap bias and bootstrap one-sided P-value
% Scaled Spearman R summands: Osamp is an nnx1 column vector
osamp=12*rxc.*ryc/((nn-1)*nn*(nn+1));
for ns=1:nsims;
  rboot(ns)=sum(osamp(ceil(nn*rand(1,nn))));
  if rboot(ns)<=rr bias=bias+1;  end;
  if rboot(ns)<=0  nle=nle+1;  end;
  end;
rboot=sort(rboot);
bmed=median(rboot);
blim1=rboot(floor(nsims*0.025));
blim2=rboot(ceil(nsims*0.975));
fprintf('    %6.4f    %6.4f    %6.4f    Bias: %6.4f\n', blim1,bmed,blim2, bias/nsims);
fprintf('Bootstrap P-value for H_0:R=0:  P=%6.4f (2-sided)\n', 2*nle/nsims);
 
