/* Compare three types of correlation measures for
     paired data (X_i,Y_i) (like height and weight)
     (Pearson rho, Spearman R, and Kendall tau)
   Also, jackknife and bootstrap confidence intervals for all three */

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
/* Prototypes for functions in  statlib.c :
   Compile by e.g.   gcc -o ProgName Progname.c statlib.c   with
     ProgName.c, statlib.c, and statlib.h in the default directory. */
#include "statlib.h"

#define DTOL  1e-12   /* Minimum floating-point error */

#define MAXN 50       /* Allow for up to this many pairs */

/* The `const' says that the program should not change these values */
/*   without a compiler crash or at least a very strong warning. */

/* Text for this example */
const char *title_data=
 "Prepacked Lubricants and Survival Times for Engines\n"
 "(Source of data unknown.)";

/* Number of pairs (subjects) */
const int nn=20;

const char *xname="Lubricant", *yname="Survival";

const double xx[] =
 {  29,  21,  30,  45,  48,     92,   74,   61,   99,   24,
    84,  37,  94,  56,  67,     79,   70,   93,   76,   72 };
 
const double yy[] =
 {  23,  36,  67, 104, 147,    164,  247,  304,  355,  408,
   429, 489, 504, 925, 980,   1254, 2961, 5351, 9915, 11301 };
 
/* Default starting random-number seed and the number of */
/*   permutations in permutation procedures */
unsigned long wseed=123456;
int nsims=100000;

/* Number of bootstrap replications in bootstrap procedures */
#define NBOOT 10000

int nboot=NBOOT;

/* Prototypes to calculate the three types of correlation coefficients */
void do_pearson(int nn, const double *xx, const double *yy);

/* rx[] and ry[] are ranks within each sample */
void do_spearman(int nn, const double *rx, const double *ry);

/* Kendall correlation coefficient and P-values */
/* Returns the large-sample Z statistic and the permutation P-value CI */
void do_kendall(int nn, const double *xx, const double *yy,
  double *zkendptr, double *p1ptr, double *p2ptr, double *p3ptr);

/* Pearson correlation coefficient revisited */
void do_pearson_revisited(int nn, const double *xx, const double *logyy);

/* Jackknife and Bootstrap for Pearson correlation coefficient */
void do_jackknife_pearson(int nn, const double *xx, const double *yy);
void do_bootstrap_pearson(int nn, const double *xx, const double *yy);

/* Bootstrap for Kendall tau */
void do_bootstrap_kendall(int nn, const double *xx, const double *yy, 
  double zkendstat, double p1kend, double p2kend, double p3kend);

void show_data (int nn, const double *xx, const double *yy, 
    const char *xname, const char *yname);
void show_data_and_midranks (int nn, const double *xx, const double *yy, 
    const double *rx, const double *ry,
    const char *xname, const char *yname);
void show_logdata (int nn, const double *xx, const double *yy, 
    const char *xname, const char *yname);

int main(int argc, char **argv)
 { int i;
   /* Midranks and tiegroups for each sample */
   double rx[MAXN], ry[MAXN], logyy[MAXN];
   double zkendstat, p1kend,p2kend,p3kend;
   /* Helpful for segmentation faults and some other fatal errors */
   /* This function is in statlib.c */
   set_newsignals();

   printf (
   "Three different types of correlation coefficients:\n"
   "  Pearson rho, Spearman R, and Kendall tau\n"
   "Also, jackknife and bootstrap confidence intervals for all three.\n");
   /* Initialize random-number seed */
   /* Enter `(ProgramName) <Number>' to set starting RN seed=<Number> */
   if (argc>=2) wseed=atol(argv[1]);
   wseed=setrand(wseed);   /* If wseed=0, wseed=value from system clock */
   printf ("\nStarting random-number seed: %lu\n",wseed);
   printf ("(Enter a number on the command line "
     "to specify a different seed.)\n");

   /* Describe and display the data */
   printf ("\n%s\n\n", title_data);
   show_data(nn,xx,yy, xname,yname);

   /* Make sure that nn is big enough */
   if (nn<4) { printf ("\nN=%d is too small!\n",nn);  exit(1); }

   /* Calculate and analyze Pearson correlation coefficient: */
   do_pearson(nn,xx,yy);

   /* Find midranks and tiegroup sizes for both xx[] and yy[] */
   /* makeranks() is in statlib.c */
   /* Tie-group information not needed */
   makeranks(rx, xx,nn, NULL,NULL);
   makeranks(ry, yy,nn, NULL,NULL);

   printf ("\nData and within-sample midranks for both samples:\n\n");
   show_data_and_midranks(nn, xx,yy, rx,ry, xname,yname);

   /* Spearman correlation coefficient: */
   do_spearman(nn,rx,ry);

   /* Kendall correlation coefficient: */
   /* Returns the large-sample Kendall Z statistic to avoid recomputing it */
   /*   as well as the confidence interval for the permutation-test P-value */
   do_kendall(nn,xx,yy, &zkendstat, &p1kend, &p2kend, &p3kend);

   /* Pearson correlation coefficient, revisited */
   printf ("\nREVISITING the Pearson correlation coefficient:\n");
   printf ("For Xs and log(Y)s:\n");
   /* logyy[] is log(Y) */
   for (i=0; i<nn; i++)  logyy[i]=log(yy[i]);
   do_pearson_revisited(nn,xx,logyy);
   printf ("\nLog-transformed data:\n");
   show_logdata(nn,xx,logyy, xname,"LogSurv");

   /* Jackknifes and bootstraps? */
   printf ("\nCONFIDENCE INTERVERVALS FOR RHO, R, and TAU?\n");
   printf ("ALTERNATIVE P-VALUES?\n");
   printf ("JACKKNIFES AND BOOTSTRAPS:\n");

   printf ("\nUsing the jackknife:\n");
   printf ("Pearson rho for Xs and Ys:\n");
   printf ("(See Section 5.2 in text for a discussion of jackknifes,\n");
   printf ("   or else the Jackknife Handout on the Math408 Web site.)\n");
   do_jackknife_pearson(nn,xx,yy);
   printf ("\nPearson rho for Xs and log(Y)s:\n");
   do_jackknife_pearson(nn,xx,logyy);

   printf ("\nUsing the bootstrap:\n");
   printf ("Pearson rho for Xs and Ys:\n");
   printf ("(See Section 8.4 in text for a discussion of bootstraps,\n");
   printf ("   or else the Bootstrap Handout on the Math408 Web site.)\n");
   do_bootstrap_pearson(nn,xx,yy);
   printf ("\nPearson rho for Xs and log(Y)s:\n");
   do_bootstrap_pearson(nn,xx,logyy);

   printf ("\nBootstraps for Kendall tau for Xs and Ys:\n");
   do_bootstrap_kendall(nn,xx,yy, zkendstat, p1kend,p2kend,p3kend);

   return 0; }


void pline(int aa)
 { char nbuff[88];
   int i;
   for (i=0; i<aa; i++) nbuff[i]='-';  nbuff[i]=0;
   printf ("%s\n",nbuff); }


void show_data (int nn, const double *xx, const double *yy, 
    const char *xname, const char *yname)
 { int i, nhalf=nn/2;
   if ( (nn%2)!=0) nhalf++;
   /* Show the data in two columns */
   printf ("    %-9s  %-9s     %-9s  %-9s\n", xname,yname, xname,yname);
   pline(4*9+12);
   for (i=0; i<nhalf; i++)
    { int i2=i+nhalf;  /* Index of second column */
      /* The first column */
      printf ("    %2d.  %3g %6g", i+1, xx[i], yy[i]);
      /* If nn is odd, the last line will have only one entry */
      if (i<nhalf-1 || (nn%2)==0)
        printf ("          %2d.  %3g %6g", i2+1, xx[i2], yy[i2]);
      printf("\n"); }
   pline(4*9+12); }


void show_data_and_midranks (int nn, const double *xx, const double *yy, 
    const double *rx, const double *ry,
    const char *xname, const char *yname)
 { int i, nhalf=nn/2;
   if ( (nn%2)!=0) nhalf++;
   /* Show the data in two columns */
   printf ("        %-10s   %-10s        %-10s   %-10s\n", 
      xname,yname, xname,yname);
   pline(4*13+8);
   for (i=0; i<nhalf; i++)
    { int i2=i+nhalf;
      /* The first column */
      printf ("    %2d.  %3g (%2g) %6g (%2g)", 
        i+1, xx[i],rx[i], yy[i],ry[i]);
      /* If nn is odd, the last line will have only one entry */
      if (i<nhalf-1 || (nn%2)==0)
        printf ("      %2d.  %3g (%2g) %6g (%2g)",
          i2+1, xx[i2],rx[i2], yy[i2],ry[i2]);
      printf("\n"); }
   pline(4*13+8); }


void show_logdata (int nn, const double *xx, const double *yy, 
    const char *xname, const char *yname)
 { int i, nhalf=nn/2;
   if ( (nn%2)!=0) nhalf++;
   /* Show the data in two columns */
   printf ("    %-9s  %-9s     %-9s  %-9s\n", xname,yname, xname,yname);
   pline(4*9+12);
   for (i=0; i<nhalf; i++)
    { int i2=i+nhalf;  /* Index of second column */
      /* The first column */
      printf ("    %2d.  %3g %6.2f", i+1, xx[i], yy[i]);
      /* If nn is odd, the last line will have only one entry */
      if (i<nhalf-1 || (nn%2)==0)
        printf ("          %2d.  %3g %6.2f", i2+1, xx[i2], yy[i2]);
      printf("\n"); }
   pline(4*9+12); }




/* Show P-value and 95% CI for true P-value for a permutation test */
/* Done in a subroutine to guarantee the same format */

void showpval (int nbig, int nsims)
 { double pwid, pval=(double)nbig/nsims;
   pwid=1.960*sqrt(pval*(1-pval)/nsims);
   printf ("    P=%d/%s=%g   95%% CI for true P: (%8.6f, %8.6f)\n",
     nbig,commnum(nsims), pval, pval-pwid, pval+pwid); }


/* Analyze the Pearson correlation coefficient */

/* Return a single Pearson correlation coefficient */

double get_pearson_corr(int nn, const double *xx, const double *yy)
 { int i;
   double sx,sy, xmean,ymean, sxx,sxy,syy;
   sx=sy=sxx=sxy=syy=0.0;
   /* Find the means */
   for (i=0; i<nn; i++)
    { sx+=xx[i];  sy+=yy[i]; }
   xmean=sx/nn;  ymean=sy/nn;
   /* Find the two variances and the covariance */
   for (i=0; i<nn; i++)
    { double cx=xx[i]-xmean,  cy=yy[i]-ymean;
      sxx+=cx*cx;  sxy+=cx*cy;  syy+=cy*cy; }
   if (sxx<DTOL || syy<DTOL)
    { printf("Variance of X or Y is zero!\n");  exit(1); }
   return  sxy/sqrt(sxx*syy); }


/* Permutation test for a Pearson-rho-based score carried out */
/*   by fixing the Xs and shuffling the Ys */

void shuffle_pearson_display(double obs_rho, int nn, 
     const double *xx, const double *yy)
 { int i, nbig,ns;
   double ytemp[MAXN], pp,pwid;
   /* Copy Ys to a buffer that is safe to shuffle */
   for (i=0; i<nn; i++) ytemp[i]=yy[i];
   /* Carry out the permutation test using `nsims' simulations */
   for (nbig=ns=0; ns<nsims; ns++)
    { double newrho;
      shuffle_dub(ytemp,nn);
      newrho=get_pearson_corr(nn,xx,ytemp);
      /* Be conservative about machine round-off error */
      if (fabs(newrho)>=fabs(obs_rho)-DTOL) nbig++; }
   /* Display permutation test P-value and 95% CI */
   pp=(double)nbig/nsims;
   pwid=1.960*sqrt(pp*(1-pp)/nsims);
   printf ("  P=%d/%s=%g  95%% CI for true P: (%6.4f, %6.4f, %6.4f)\n",
     nbig,commnum(nsims), pp, pp-pwid,pp,pp+pwid); }



void do_pearson(int nn, const double *xx, const double *yy)
 { double rho,tt;
   printf ("\nPearson correlation coefficient (n=%d):\n", nn);
   rho=get_pearson_corr(nn,xx,yy);
   tt=rho*sqrt((nn-2)/(1-rho*rho));
   printf ("  Rho=%g   T=%g  Classical P=%7.5f (Student's t, df=%d)\n",
     rho, tt, pvttest(nn-2,tt), nn-2);

   /* Do `nsims' permutations */
   printf ("Permutation test based on %s permutations:\n", commnum(nsims));
   printf ("  Each permutation fixes the Xs, permutes the Ys\n");
   shuffle_pearson_display(rho, nn,xx,yy);
        }



void do_spearman(int nn, const double *rx, const double *ry)
 { double rr,zz;
   /* The Spearman correlation coefficient (with tie corrections) */
   /*   is the same as the Pearson correlation for within-sample ranks */
   /* The large-sample approximation is different, however */
   printf ("\nSpearman correlation coefficient (n=%d):\n", nn);
   rr=get_pearson_corr(nn,rx,ry);
   zz=rr*sqrt(nn-1);
   printf ("  R=%g   Z=%g   Large-sample P=%7.5f (two-sided normal)\n",
     rr, zz, normpval(zz));
   /* Permutation test based on permutations of the Y-values, */
   /*   keeping the X-values fixed: */
   printf ("Doing %s permutations fixing RankXs, permuting RankYs:\n", 
      commnum(nsims));
   /* Everything after this point is the same with the within-sample */
   /*   ranks/midranks replacing the original sample */
   shuffle_pearson_display(rr, nn,rx,ry); }



/* The Pearson correlation coefficient revisited */

void do_pearson_revisited(int nn, const double *xx, const double *logyy)
 { double rho,tt;
   rho=get_pearson_corr(nn,xx,logyy);
   tt=rho*sqrt((nn-2)/(1-rho*rho));
   printf ("  Rho=%g   T=%g  Classical P=%7.5f (Student's t, df=%d)\n",
     rho, tt, pvttest(nn-2,tt), nn-2);

   /* Do `nsims' permutations */
   printf ("  Doing %s permutations fixing Xs, permuting log(Y)s:\n", 
     commnum(nsims));
   shuffle_pearson_display(rho, nn,xx,logyy);
      }



int get_kendall_kk(int nn,
   const double *xx, const double *yy, int *ntieptr)
 { int i,j, kk,nties;
   kk=nties=0;
   for (i=0; i<nn; i++)
    { double xi=xx[i], yi=yy[i];
      for (j=i+1; j<nn; j++)
       { double ww;
         if ( (ww = (xx[j]-xi)*(yy[j]-yi))==0.0) nties++;
         else if (ww > 0.0)  kk++;
         else kk--; }}
   if (ntieptr!=NULL) *ntieptr=nties;
   return kk; }


/* The next function is probably the easiest way to handle tie corrections */
/*   for the large-sample approximation of the Kendall K statistic */
/* Note that it fills three integers in an integer list with tie sums */

void get_kendall_tiesums (int *kt, int ntg, const int *tielist)
 { int i;
   kt[0]=kt[1]=kt[2]=0;  /* Initialize at zero */
   for (i=0; i<ntg; i++)
    { int tt=tielist[i];
      kt[0]+=tt*(tt-1);  kt[1]+=tt*(tt-1)*(tt-2);
      kt[2]+=tt*(tt-1)*(2*tt+5); }}


double get_kendall_var (int nn,
  int nxtg, const int *xtielist, int nytg, const int *ytielist)
 { double vark;  int tx[5], ty[5];
   get_kendall_tiesums(tx,nxtg,xtielist);
   get_kendall_tiesums(ty,nytg,ytielist);
   vark = (double)(nn*(nn-1)*(2*nn+5) - tx[2] - ty[2])/18
     + (double)(tx[1]*ty[1])/(9*nn*(nn-1)*(nn-2))
     + (double)(tx[0]*ty[0])/(2*nn*(nn-1));
   return vark; }


/* Returns large-sample Z score and CI for permutation-test P-value */

void do_kendall(int nn, const double *xx, const double *yy,
  double *zkendptr, double *p1ptr, double *p2ptr, double *p3ptr)
 { int i,kk,nties, nbig,ns;
   int nxtg,nytg, xtielist[MAXN], ytielist[MAXN];
   double tau,vark,zz, mranks[MAXN], yytemp[MAXN], pp,pwid;
   printf ("\nKendall correlation coefficient (n=%d):\n", nn);
   kk = get_kendall_kk(nn,xx,yy,&nties);
   tau = (double)(2*kk)/(nn*(nn-1));
   printf ("  Kendall K=%d  Kendall tau=%g  %d tie%s\n", 
     kk, tau, nties, (nties==1) ? "" : "s");
   /* The tie correction for the large-sample approximation for */
   /*   the P-value needs the tie-group structure for each sample */
   /* Information written to mranks[] is not used */
   makeranks(mranks, xx,nn, xtielist,&nxtg);
   makeranks(mranks, yy,nn, ytielist,&nytg);
   /* The variance itself needs only n and the tiegroups */
   vark=get_kendall_var(nn, nxtg,xtielist, nytg,ytielist);
   /* The large-sample Z score */
   zz=kk/sqrt(vark);
   printf ("  Large-sample Z=%6.4f  P=%6.4f (2-sided, tie-corrected)\n",
     zz, normpval(zz));
   /* Permutation test based on permutations of the Y-values, */
   /*   keeping the X-values fixed: */
   /* Return values of K; don't worry about taus in permutation test */
   for (i=0; i<nn; i++) yytemp[i]=yy[i];
   printf ("  Doing %s permutations fixing Xs, permuting Ys:\n", 
      commnum(nsims));
   for (nbig=ns=0; ns<nsims; ns++)
    { int newkk;
      shuffle_dub(yytemp,nn);
      newkk=get_kendall_kk(nn,xx,yytemp,NULL);
      /* get_kendall_kk() returns an integer, so no round-off error */
      /* abs() returns the absolute value of an integer */
      /* abs() is a standard C function */
      /* These are integers, so no round-off error */
      if (abs(newkk)>=abs(kk)) nbig++; }
   /* Display permutation test P-value and 95% CI */
   /* Display permutation test P-value and 95% CI */
   pp=(double)nbig/nsims;
   pwid=1.960*sqrt(pp*(1-pp)/nsims);
   printf ("  P=%d/%s=%g  95%% CI for true P: (%6.4f, %6.4f, %6.4f)\n",
     nbig,commnum(nsims), pp, pp-pwid,pp,pp+pwid);
   if (zkendptr!=NULL) *zkendptr=zz;
   if (p1ptr!=NULL) *p1ptr=pp-pwid;
   if (p2ptr!=NULL) *p2ptr=pp;
   if (p3ptr!=NULL) *p3ptr=pp+pwid; }




void do_jackknife_pearson(int nn, const double *xx, const double *yy)
 { int i,j,r;
   /* Jackknife delete-1 and pseudo-values */
   double rho, del1value[MAXN], psvalue[MAXN];
   double jsum, jmean,jvar,jstd,tt, tcrit,jzwid,jtwid;
   rho=get_pearson_corr(nn,xx,yy);
   printf ("Jackknife of Pearson rho (n=%d, original rho=%6.4f)\n", nn,rho);
   /* Construct the jackknife delete-1 and pseudo-values */
   for (i=0; i<nn; i++)
    { /* Compute the delete-1 jackknife value */
      double newrho, xtemp[MAXN], ytemp[MAXN];
      /* First construct the delete-1 dataset */
      for (r=j=0; j<nn; j++) if (j!=i)
       { xtemp[r]=xx[j];  ytemp[r]=yy[j];  r++; }
      /* The delete-1 pearson rho */
      /* Note that we must use nn-1 instead of nn */
      newrho=del1value[i]=get_pearson_corr(nn-1,xtemp,ytemp);
      /* The i-th psuedo value */
      psvalue[i]=nn*rho-(nn-1)*newrho; }
   printf ("Delete-1 rho values:\n");
   for (i=0; i<nn; i++)
    { if (i>0 && (i%5)==0) printf("\n");
      printf ("  %7.4f(%d)%s",del1value[i],i+1,(i<10)?" ":""); } 
   printf("\n");
   printf ("Jackknife pseudo-values for rho:\n");
   for (i=0; i<nn; i++)
    { if (i>0 && (i%5)==0) printf("\n");
      printf ("  %7.4f(%d)%s",psvalue[i],i+1,(i<10)?" ":""); } printf("\n");
   /* Find mean and standard deviation of pseudo-values */
   for (jsum=0.0, i=0; i<nn; i++) jsum+=psvalue[i];
   jmean=jsum/nn;
   for (jsum=0.0, i=0; i<nn; i++) jsum+=(psvalue[i]-jmean)*(psvalue[i]-jmean);
   jvar=jsum/(nn-1);
   jstd=sqrt(jvar/nn);  tt=jmean/jstd;
   /* Critical value for Student-t with df=nn-1 */
   /* Needed for 95% Student-t CI for true rho */
   /* This function is in statlib.c */
   tcrit=crit_tdist(nn-1,0.05);
   /* Half-widths of 95% Normal and Student-t confidence intervals */
   jzwid=1.960*jstd;  jtwid=tcrit*jstd;
   printf ("For original rho=%6.4f:\n",rho);
   printf ("Jackknife mean, 95%% CIs, and jackknife (two-sided) P-values:\n");
   printf (
   "  Normal:     (%6.4f, %6.4f, %6.4f)  Z=%6.4f  P=%7.5f\n",
     jmean-jzwid, jmean, jmean+jzwid, tt, normpval(tt));
   printf (
   "  Student-t:  (%6.4f, %6.4f, %6.4f)  T=%6.4f  P=%7.5f (%d df)\n",
     jmean-jtwid, jmean, jmean+jtwid, tt, pvttest(nn-1,tt), nn-1);
                           }


double prnd (double pp)
 { return (pp>1.0) ? 1.0 : pp; }


void show_bootstrap_output (const char *varname,
    int nboot, double *zboot, double val)
 { int nhits,ns, blo,bhi;
   double bootmedian, zlo, pp,pwid;
   /* First, sort zboot[]:  qsort() is a standard C function that implements
       `quick-sort' that sorts a list of size n  in  O(n log n) steps
       One of qsort's arguments is a pointer to a `comparison function'
       that qsort uses to do the sorting
       `dubsort' in statlib.c is for sorting doubles in increasing order */
   qsort(zboot,nboot, sizeof(double),dubsort);
   /* Find the median */
   if (nboot%2)  bootmedian=zboot[nboot/2];
   else bootmedian = (zboot[nboot/2] + zboot[(nboot/2)-1])/2;
   /* The 95% percentile bootstrap CI is the middle 95% of the bootstrap */
   /*   sample, based on order statistics */
   /* We follow the suggestion in the text (p388-389), following */
   /*   Efron and Tibshirani (1993), about the best way to calculate these */
   zlo=0.250*nboot;
   if (zlo-floor(zlo)<0.001)  /* If zlo is an integer */
     blo=(int)floor(zlo)-1;   /* The zlo-th array entry */
   else blo=(int)floor(0.025*(nboot+1))-1;
   bhi=nboot-1-blo;     /* Number k from the top if blo is #k from the bottom */
   printf ("  Bootstrap median and 95%% CI for %s:  (%6.4f, %6.4f, %6.4f)\n",
     varname, zboot[blo], bootmedian, zboot[bhi]);
   /* Bootstrap test of var!=0: */
   /* Count bootstrap replications where newval has the opposite sign as val */
   for (nhits=ns=0; ns<nboot; ns++)
     if (zboot[ns]*val <= 0.0)  nhits++;
   pp=(double)nhits/nboot;
   /* The 95% CI for pval is the CI for pp scaled up */
   /* prnd() defined above rounds p-values down to 1.0 if necessary */
   pwid=1.960*sqrt(pp*(1-pp)/nboot);
   printf ("  Bootstrap test for %s!=0:\n", varname);
   printf ("    One-sided P:  P=%d/%s = (%6.4f, %6.4f, %6.4f)\n",
     nhits,commnum(nboot), pp-pwid,pp,pp+pwid);
   printf ("    Two-sided P:            (%6.4f, %6.4f, %6.4f)\n",
     prnd(2*(pp-pwid)), prnd(2*pp), prnd(2*(pp+pwid)));
      }



double zboot[NBOOT+3];

void do_bootstrap_pearson(int nn, const double *xx, const double *yy)
 { int i,ns;
   double xxtemp[MAXN], yytemp[MAXN];
   double rho=get_pearson_corr(nn,xx,yy);
   printf ("Bootstrap of Pearson rho (%s replications, original rho=%6.4f)\n",
     commnum(nboot), rho);
   /* Do `nboot' bootstrap replications for the value of K */
   /* Multiply by 2/(nn*(nn-1)) to convert to Kendall tau coefficients */
   for (ns=0; ns<nboot; ns++)
    { /* Construct a bootstrap replication in (xxtemp,yytemp) */
      for (i=0; i<nn; i++)
       { int ran=nrand(nn);
         xxtemp[i]=xx[ran];  yytemp[i]=yy[ran]; }
      /* Put the Pearson correlation in zboot[ns] */
      zboot[ns]=get_pearson_corr(nn,xxtemp,yytemp); }
   /* Analyze bootstrap output */
   show_bootstrap_output("rho", nboot,zboot, rho);
      }


void do_bootstrap_kendall(int nn, const double *xx, const double *yy, 
  double zkendstat, double p1kend, double p2kend, double p3kend)
 { int i,ns, obskk;
   double tau, xxtemp[MAXN], yytemp[MAXN];
   double nfac=2/((double)nn*(nn-1));
   obskk=get_kendall_kk(nn,xx,yy,NULL);
   tau=nfac*obskk;
   printf ("\nBootstrap of Kendall tau (%s replications, tau=%g)\n",
     commnum(nboot), tau);
   printf ("  Recall: Large-sample Z=%6.4f  P=%6.4f "
     "(2-sided, tie-corrected)\n", zkendstat,normpval(zkendstat));
   printf ("  95%% CI for true permutation-test P-value:  (%6.4f, %6.4f, %6.4f)\n",
     p1kend, p2kend, p3kend);
   /* Do `nboot' bootstrap replications for the value of K */
   /* Multiply by 2/(nn*(nn-1)) to convert to Kendall tau coefficients */
   for (ns=0; ns<nboot; ns++)
    { int kk;
      /* Construct a bootstrap replication in (xxtemp,yytemp) */
      for (i=0; i<nn; i++)
       { int ran=nrand(nn);
         xxtemp[i]=xx[ran];  yytemp[i]=yy[ran]; }
      /* Put the Kendall tau correlation in zboot[ns] */
      kk=get_kendall_kk(nn,xxtemp,yytemp,NULL);
      zboot[ns]=nfac*kk; }
   qsort(zboot,nboot, sizeof(double),dubsort);
   show_bootstrap_output("tau", nboot,zboot, tau);
      }