/* Comparison of 3 different methods of linear regression:
     (1) Least-squares
     (2) ``Theil'' linear regression
     (3) Rank regression       */

#define PROG_NAME "RANKREGRESS"
#define YOUR_NAME "YOURNAME"

#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)";
const char *xname="Lubricants";
const char *yname="Survival Time";
const char *ynameshort="Survival";

/* Number of treatment groups and subjects */
const int nn=20;

const double xx[MAXN] =
 {  29,  21,  30,  45,  48,     92,   74,   61,   99,   24,
    84,  37,  94,  56,  67,     79,   70,   93,   76,   72 };
 
const double yy[MAXN] =
 {  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 */
/*   number of replications in bootstrap procedures */
unsigned long wseed=123456;

#define NBOOT 10000     /* This many bootstrap replications */

int nboot=NBOOT;

/* Describe and display the data */
void display_data(const char *title_data, const double *xx,
    const double *yy, int nn, const char *xname, const char *yname);

/* Display three regression lines */
void do_regressions(const char *Ytype, const double *xx, const double *yy,
    int nn, const char *xname, const char *yname);

/* Ordinary (least-squares) regression and P-values */
/*   and bootstrap output */
void ci_pvals_lsq_regr(const char *Ytype,
   const double *xx, const double *yy, int nn);
void bootstrap_beta_lsq_regression (const char *Ytype,
   const double *xx, const double *yy, int nn);

/* Theil regression and P-values */
/*   and bootstrap output */
void ci_pvals_theilregr(const char *Ytype,
   const double *xx, const double *yy, int nn);
void bootstrap_beta_theil_regression (const char *Ytype,
   const double *xx, const double *yy, int nn);

/* Rank regression and P-values */
/*   and bootstrap output */
void bootstrap_beta_rank_regression (const char *Ytype,
   const double *xx, const double *yy, int nn);

/* Display the summary lines so far */
void Dump_Regression_Summary(void);


int main(int argc, char **argv)
 { int i;
   double logyy[MAXN];
   /* Helpful for segmentation faults and some other fatal errors */
   /* This function is in statlib.c */
   set_newsignals();

   printf (
   "Comparison of 3 different methods of linear regression:\n"
   "  (1) Least-squares\n"
   "  (2) ``Theil'' linear regression\n"
   "  (3) Rank regression\n");
   printf ("Output of program   %s\n", PROG_NAME);
   printf ("Written or modified by  %s\n", YOUR_NAME);
   /* 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 the data */
   display_data(title_data, xx,yy,nn, xname,ynameshort);

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

   printf (
   "\n(The last two columns below are the errors\n"
   "    MADf = (1/n)Sum(i=1,n) |Y_i-beta*X_i-mu|  and\n"
   "    RMS  = Root( (1/n)Sum(i=1,n)(Y_i-beta*X_i-mu)^2 )\n");

   /* Display regression lines for all three regressions */
   do_regressions("Y", xx,yy,nn, xname,yname);

   /* Try log-transformed Y? */
   /* Put log(Y) values in an array */
   for (i=0; i<nn; i++)  logyy[i]=log(yy[i]);
   do_regressions("log(Y)", xx,logyy,nn, xname,yname);

   printf ("\nCONFIDENCE INTERVALS and P-VALUES for SLOPE:\n");

   /* Ordinary (least-squares) regression */
   ci_pvals_lsq_regr("Y", xx,yy,nn);
   bootstrap_beta_lsq_regression("Y", xx,yy,nn);
   /* For log-transformed Y-values as well */
   ci_pvals_lsq_regr("log(Y)", xx,logyy,nn);

   /* Theil regression and P-values and related bootstraps */
   printf("\n");
   ci_pvals_theilregr("Y", xx,yy,nn);
   bootstrap_beta_theil_regression("Y", xx,yy,nn);

   /* Bootstraps of rank regression and related bootstraps */
   printf("\n");
   bootstrap_beta_rank_regression("Y", xx,yy,nn);

   /* Display a summary */
   printf ("\nSUMMARY of estimates, confidence intervals, and P-values:\n");
   Dump_Regression_Summary();

   show_randnums_used();

   return 0; }


/* Draw a dotted line */
void pline(void)
 { int j;
   printf ("---");
   for (j=0; j<6; j++) printf ("---------");
   printf ("\n"); }


/* Describe and display the data */

void display_data(const char *title_data, const double *xx,
    const double *yy, int nn, const char *xname, const char *yname)
 { int i, nhalf;
   printf ("\n%s\n", title_data);
   nhalf=nn/2;  if ( (nn%2)!=0) nhalf++;
   /* Show the data in two columns */
   printf ("\n    %-9s  %-9s     %-9s  %-9s\n", xname,yname, xname,yname);
   pline();
   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(); }


/* Display regression line with errors */

/* Least-squares regression */
double getlsqbeta(const double *xx, const double *yy,
    int nn, double *muptr);

/* Theil regression */
double gettheilbeta(const double *xx, const double *yy,
    int nn, double *muptr);

/* Rank regression */
double getrankbeta(const double *xx, const double *yy,
    int nn, double *muptr);


/* Display results for one simple linear regression */
void regline (const char *regname, const char *Ytype, 
    double beta, double mu, const double *xx, const double *yy, int nn)
 { char eqbuff[288];
   int i;  double sumabs,sum2, madf,rms;
   /* Write the regress line to eqbuff[] */
   sprintf (eqbuff, "%s = %7.4f X + %8.3f", Ytype, beta, mu);
   /* Finding the sum of the absolute values of the residuals */
   /*  and of the squares of the residuals */
   sumabs=sum2=0.0;
   for (i=0; i<nn; i++)
    { double res = yy[i] - beta*xx[i] - mu;
      sumabs+=fabs(res);  sum2+=res*res; }
   madf=sumabs/nn;  rms=sqrt(sum2/nn);
   printf ("  %-10s  %-30s  %7.2f   %7.2f\n",
     regname, eqbuff, madf, rms); }



void do_regressions(const char *Ytype, const double *xx, const double *yy,
    int nn, const char *xname, const char *yname)
 { double beta,mu;
   printf ("\nREGRESSIONS for %s(%s) on X(%s):\n", Ytype,yname,xname);
   printf ("%46s   MADf      RMS\n", "");
   /* Ordinary (least-squares) regression */
   beta=getlsqbeta(xx,yy,nn, &mu);
   regline ("LeastSqrs:", Ytype, beta,mu, xx,yy,nn);
   /* Theil regression */
   beta=gettheilbeta(xx,yy,nn, &mu);
   regline ("Theil:", Ytype, beta,mu, xx,yy,nn);
   /* Rank regression */
   beta=getrankbeta(xx,yy,nn, &mu);
   regline ("RankRegr:", Ytype, beta,mu, xx,yy,nn); }


/* Return the slope of the least-squares regression */

double getlsqbeta(const double *xx, const double *yy, int nn, double *muptr)
 { int i;
   double sx,sy, xmean,ymean, sxx,sxy, beta;
   sx=sy=sxx=sxy=0.0;
   /* Find the X and Y means */
   for (i=0; i<nn; i++) { sx+=xx[i];  sy+=yy[i]; }
   xmean=sx/nn;  ymean=sy/nn;
   /* Find the X-variance and covariance */
   /* Note that cx and cy below are centered at their means */
   for (i=0; i<nn; i++)
    { double cx=xx[i]-xmean,  cy=yy[i]-ymean;
      sxx+=cx*cx;  sxy+=cx*cy; }
   if (sxx<DTOL)
    { printf("The variance of X is zero!\n");  exit(1); }
   beta=sxy/sxx;
   if (muptr!=NULL) *muptr=ymean-beta*xmean;
   return  beta; }


/* Return the slope of the Theil regression line */

double getmed (const double *array, int num)
 { if ((num%2)==1)  return array[num/2];
   else return (array[num/2] + array[(num/2)-1])/2; }


double set_theil_diffratios (double *rbuff,
  const double *xx, const double *yy, int nn)
 { int i,j,nd;
   for (nd=i=0; i<nn; i++)  for (j=i+1; j<nn; j++)
    { /* Ignore pairs with X-X ties */
      if (fabs(xx[i]-xx[j])>DTOL)  /* Unless an X-X tie */
        rbuff[nd++] = (yy[j] - yy[i])/(xx[j] - xx[i]); }
   if (nd<=0) /* All X-values are the same */
    { printf ("Find Theil values: All X values are constant!\n");
      exit(1); }
   /* Sort the difference ratios */
   /* qsort() is a standard C function */
   /* sortdub() is in statlib.c */
   qsort(rbuff,nd, sizeof(double),sortdub);
   /* Return the median */
   return  getmed(rbuff,nd); }


double get_MedianWalshSumsResids(const double *xx, const double *yy,
    int nn, double beta)
 { int i,j,r;
   double resid[MAXN], walsh[MAXN*MAXN];
   /* Put residuals in resid[] */
   for (i=0; i<nn; i++)
     resid[i] = yy[i]-beta*xx[i];
   /* Put Walsh sums of residuals in walsh[] */
   for (r=i=0; i<nn; i++) for (j=i; j<nn; j++)
     walsh[r++] = (resid[i]+resid[j])/2;
   /* Return median of Walsh sums */
   /* First, sort walsh[] */
   qsort(walsh,r, sizeof(double),sortdub);
   return getmed(walsh,r); }


double gettheilbeta(const double *xx, const double *yy,
    int nn, double *muptr)
 { double beta, rdiff[MAXN*MAXN];
   beta = set_theil_diffratios(rdiff, xx,yy,nn);
   if (muptr!=NULL) *muptr = get_MedianWalshSumsResids(xx,yy,nn,beta);
   return beta; }


/* Return the slope of the rank regression line */

/* A structure used in getrankbeta() */

typedef struct { double yratio, xdiff; }  WRATIO, *PWRATIO;

/* Sort function for used with qsort() */
/* Sorts WRATIO structions in order of increasing yratios, */
/*   then in order of increasing xdiffs */

int sortwratio(const void *p, const void *q)
 { double xx;
   if ( (xx = ((PWRATIO)p)->yratio - ((PWRATIO)q)->yratio)>DTOL)  return 1;
   else if (xx<-DTOL) return -1;
   if ( (xx = ((PWRATIO)p)->xdiff - ((PWRATIO)q)->xdiff)>DTOL)  return 1;
   else if (xx<-DTOL) return -1;
   return 0; }

/* Return the rank regression estimate of beta */
/* This is the value of beta that minimizes  */
/*   D(beta) = Sum (R_i(beta)-(n+1)/2)*(Y_i-beta*X_i)   */
/* where R_i(beta) are the ranks of the residuals       */
/*   Y_i^* = Y_i-beta*X_i  */

double getrankbeta(const double *xx, const double *yy,
    int nn, double *muptr)
 { int i,j,nd;
   double qsum, beta, lslope, nbar;
   /* An array of WRATIO structures */
   WRATIO ww[MAXN*MAXN];  double midr[MAXN];
   /* Put midranks of X_i in midr[] */
   /* makeranks() is in statlib.c */
   makeranks(midr, xx,nn, NULL,NULL);
   /* Compute Q = Sum (Rank(X_i)-(n+1)/2) X_i  */
   qsum=0.0;  nbar=(double)(nn+1)/2;
   for (i=0; i<nn; i++) qsum += (midr[i]-nbar)*xx[i];
   /* Collect differences and difference ratios in ww[] */
   for (nd=i=0; i<nn; i++)  for (j=i+1; j<nn; j++)
     if (fabs(xx[j]-xx[i])>DTOL)
      { ww[nd].yratio = (yy[j]-yy[i])/(xx[j]-xx[i]);
        ww[nd].xdiff = fabs(xx[j]-xx[i]);  nd++; }
   if (nd<=0) /* All X-values are the same */
    { printf ("While finding Rank Regression slope: "
          "All X values are constant!\n");  exit(1); }
   /* Sort ww[] by increasing yratio, then by increasing xdiff */
   qsort(ww,nd, sizeof(WRATIO),sortwratio);
   /* Initialize to find the low point */
   /* lslope is iniitially the slope just before ww[0] */
   lslope=-qsum;
   for (i=0; i<nd; i++)
    { /* rslope is the derivative just after ww[i] */
      double rslope=lslope+ww[i].xdiff;
      /* Make sure that flat spots are treated as flat spots. */
      if (lslope<-DTOL && DTOL<rslope)
       { beta=ww[i].yratio;  break; }
      /* In case of a flat spot */
      else if (fabs(rslope)<DTOL)
       { beta=(ww[i].yratio+ww[i+1].yratio)/2;  break; }
      /* The derivative after is the next step's derivative before */
      lslope=rslope; }
   /* Beta is the rank-regression slope */
   /* If requested, Mu is the rank-regression intercept */
   if (muptr!=NULL) *muptr = get_MedianWalshSumsResids(xx,yy,nn,beta);
   return beta; }


/* Prototypes of routines to add comments to the Summary Table: */

void addrcomm (const char *RegName, const char *Ytype);
void addrsumm (const char *RegName, double beta,
   double cilow, double cihigh, double pval, int df);

/* Ordinary (least-squares) regression */

void ci_pvals_lsq_regr(const char *Ytype,
   const double *xx, const double *yy, int nn)
 { int i;
   double beta,mu, sx,xbar, sse,sxx, mse;
   double bstd,cwid,clow,chigh, tstat, pval;
   /* Write a paragraph heading to the summary list */
   addrcomm("Least Squares", Ytype);
   printf (
   "\nLEAST-SQUARES REGRESSION (normal-theory) for %s = beta*X + mu\n"
   "  Confidence interval and P-value for beta\n", Ytype);
   beta=getlsqbeta(xx,yy,nn, &mu);
   /* Find Xbar */
   sx=0.0;  for (i=0; i<nn; i++) sx+=xx[i];  xbar=sx/nn;
   /* Find SSE=Sum of squared residulas and SXX=Sum(X_i-Xbar)^2 */
   sse=sxx=0.0;
   for (i=0; i<nn; i++)
    { double xdif=xx[i]-xbar, res=yy[i]-beta*xx[i]-mu;
      sxx+=xdif*xdif;  sse+=res*res; }
   /* Bstd = sqrt(Var(betahat))  */
   mse=sse/(nn-2);  bstd=sqrt(mse/sxx);
   /* Student-t statistic and P-value */
   /* pvttest() is in statlib.c */
   tstat=beta/bstd;  pval=pvttest(nn-2,tstat);
   /* Half-width of a Student-t confidence interval */
   /* crit_tdist() is in statlib.c */
   cwid=crit_tdist(nn-2,0.05)*bstd;
   clow=beta-cwid;  chigh=beta+cwid;
   printf ("  Beta=%g   95%% CI (%6.4f, %6.4f)  (%d df)\n",
     beta, clow,chigh, nn-2);
   printf ("  Pval(beta!=0): %g  (two-sided, %d df)\n", pval, nn-2); 
   /* Write a note to the summary list */
   addrsumm("Classical", beta, clow,chigh, pval, nn-2);
       }


/* Analyze bootstrap output */

void show_bootstrap_output (const char *boot_type, 
   double *zboot, int nboot, double estvalue);

/* Show least-squares bootstrap results */

void bootstrap_beta_lsq_regression (const char *Ytype,
     const double *xx, const double *yy, int nn)
 { int i,ns;
   double beta, resid[MAXN], xxtemp[MAXN],yytemp[MAXN];
   double zboot[NBOOT+3];
   printf ("\n(LEAST SQUARES) BOOTSTRAP for beta in %s = beta*X + mu:\n", 
     Ytype);
   printf ("  Medians and 95%% confidence intervals, %s replications.\n",
      commnum(nboot));
   /* NULL means that we aren't interested in mu */
   beta=getlsqbeta(xx,yy,nn, NULL);
   printf ("BOOTSTRAP ON OBSERVATIONS (least squares, beta=%g):\n", beta);
   for (ns=0; ns<nboot; ns++)
    { for (i=0; i<nn; i++)
       { int ran=nrand(nn);
         xxtemp[i]=xx[ran];  yytemp[i]=yy[ran]; }
      /* Compute beta for (xxtemp,yytemp) */
      zboot[ns] = getlsqbeta(xxtemp,yytemp,nn,NULL);  }
   /* Show bootstrap output */
   show_bootstrap_output("BootObs", zboot,nboot, beta);
   
   printf ("\nBOOTSTRAP ON RESIDUALS (least squares, beta=%g):\n", beta);
   /* Put residuals in resid[] */
   for (i=0; i<nn; i++)  resid[i] = yy[i]-beta*xx[i];
   /* Bootstrap of residuals */
   for (ns=0; ns<nboot; ns++)
    { /* Replace yytemp[] is Ys with new residuals */
      for (i=0; i<nn; i++)  
         yytemp[i] = beta*xx[i] + resid[nrand(nn)];
      /* Compute beta for (xx,yytemp) */
      zboot[ns] = getlsqbeta(xx,yytemp,nn,NULL); }
   /* Show bootstrap output */
   show_bootstrap_output("BootResid", zboot,nboot, beta);
      }


/* Display output from a bootstrap */

void show_bootstrap_output (const char *boot_type, 
   double *zboot, int nboot, double estvalue)
 { int nlsz,nbias, ns;
   double zmed,zlow,zhigh;
   double pval_one,pval,pwid,pvlow,pvhigh;
   /* Sort bootstrap replications */
   qsort (zboot, nboot, sizeof(double), sortdub);
   for (nlsz=nbias=ns=0; ns<nboot; ns++)
    { double newvar=zboot[ns];
      if (newvar<=estvalue) nbias++;
      if (newvar*estvalue<=0.0) nlsz++; }
   zmed=getmed(zboot,nboot);
   zlow=zboot[(int)(0.025*nboot)-1];
   zhigh=zboot[(int)(0.975*nboot)];
   printf ("  Beta=%6.4f  95%% CI (%6.4f, %6.4f)   Bias: %g\n",
     zmed, zlow,zhigh, (double)nbias/nboot);
   printf ("Bootstrap Pvalue for H_0:beta=0 vs. H_1:beta!=0:\n");
   /* One-sided P-value */
   pval_one=(double)nlsz/nboot;
   /* Twice the width of a 95% CI for the one-sided P-value */
   pwid=2*1.960*sqrt(pval_one*(1-pval_one)/nboot);
   /* Two-sided P-value = 2*One-sided P-value, */
   /*   so its CI is twice as large as for the one-sided P-value */
   pval=2*pval_one;  if (pval>1.0) pval=1.0;
   /* Lower and upper limits of 95% CI for pval */
   pvlow=pval-pwid;  if (pvlow<0.0) pvlow=0.0;
   pvhigh=pval+pwid;  if (pvhigh>1.0) pvhigh=1.0;
   printf ("  One-sided P:  P=%d/%s=%7.5f\n",
     nlsz,commnum(nboot), pval_one);
   printf ("  Two-sided P:  P=%7.5f   95%% CI: (%8.6f, %8.6f)\n", 
     pval, pvlow, pvhigh);
   /* Write a note to the summary list */
   addrsumm(boot_type, zmed,zlow,zhigh, pval, -1); }



/* Theil regression and P-values */

/* Return Kendall K statistic */

int getkend(const double *xx, const double *yy, int nn)
 { int i,j, kk;
   for (kk=i=0; i<nn; i++)
    { double zx=xx[i], zy=yy[i];
      for (j=i+1; j<nn; j++)
       { double ww;
         if ( (ww = (xx[j]-zx)*(yy[j]-zy))!=0.0)
           { if (ww > 0.0)  kk++;  else kk--; }}}
   return kk; }


void ci_pvals_theilregr(const char *Ytype,
     const double *xx, const double *yy, int nn)
 { int mm, j1,j2, k1,k2;
   double vark, cc,dd, th1,thmed,th2, zstat,pval;
   double rdiff[MAXN*MAXN];
   printf("\nDISTRIBUTION-FREE EXACT 95%% CI for Theil's estimator:\n");
   printf ("  Large-sample approximations,  %s = beta*X + mu:\n",
     Ytype);
   /* Write a note to the summary list */
   addrcomm("Theil Regression", Ytype);
   /* Construct sorted difference ratios */
   /* Returns the median of the array of difference ratios */
   thmed=set_theil_diffratios(rdiff, xx,yy,nn);
   /* The variance of the Kendall K statistic given H_0 */
   /* If there are any Y-Y or X-X ties, this has a more complex form, */
   /*   and rdiff[] has to be defined differently. */
   vark = (nn*(nn-1)*(2*nn+5))/18;
   /* Find the associated exact nonparametric CI */
   mm=(nn*(nn-1))/2;     /* The number of difference ratios */
   cc=1.960*sqrt(vark);  /* For -mm<=K<=+mm */
   dd=(mm-cc)/2;         /* Convert to  0<=K<=mm, lower 2.5% tail  */
   k1=(int)floor(dd);    /* Round down to next integer */
   k2=mm+1-k1;           /* Upper limit, (k1,k2) for range (1,d) */
   j1=k1-1;  j2=k2-1;    /* (j1,j2) for range (0,d-1) */
   if (j1<0) j1=0;  if (j2<0) j2=0;
   /* Upper and lower limits */
   th1=rdiff[j1];  th2=rdiff[j2];
   printf ("    Beta=%6.4f  95%% CI (%6.4f, %6.4f)\n"
           "      (attained at diff.ratio offsets %d,%d in (1,%d))\n",
              thmed, th1,th2, k1,k2,mm);
   /* Find Kendall P-value for H_0:beta=0 */
   zstat = getkend(xx,yy,nn)/sqrt(vark);
   pval=normpval(zstat);
   printf ("    Large-sample P-value for H_0:beta=0:\n");
   printf ("      Z=%g  P=%6.4f (2-sided)\n", zstat,pval);
   /* Write a note to the summary list */
   addrsumm("Approx", thmed, th1,th2, pval, -1);
        }


void bootstrap_beta_theil_regression (const char *Ytype,
     const double *xx, const double *yy, int nn)
 { int i,ns;
   double beta, resid[MAXN], zboot[NBOOT+3];
   printf ("\n(THEIL REGRESSION) BOOTSTRAP for beta in %s = beta*X + mu:\n", 
     Ytype);
   printf ("  Medians and 95%% confidence intervals, %s replications.\n",
      commnum(nboot));
   /* NULL means that we aren't interested in mu */
   beta=gettheilbeta(xx,yy,nn, NULL);
   printf ("BOOTSTRAP ON OBSERVATIONS (Theil, beta=%g):\n", beta);
   for (ns=0; ns<nboot; ns++)
    { double xxtemp[MAXN],yytemp[MAXN];
      for (i=0; i<nn; i++)
       { int ran=nrand(nn);
         xxtemp[i]=xx[ran];  yytemp[i]=yy[ran]; }
      zboot[ns] = gettheilbeta(xxtemp,yytemp,nn,NULL);  }
   /* Show bootstrap output */
   show_bootstrap_output("BootObs", zboot,nboot, beta);
   
   printf ("\nBOOTSTRAP ON RESIDUALS (Theil, beta=%g):\n", beta);
   /* Put residuals in resid[] */
   for (i=0; i<nn; i++)  resid[i] = yy[i]-beta*xx[i];
   /* Bootstrap of residuals */
   for (ns=0; ns<nboot; ns++)
    { double yytemp[MAXN];
      /* Replace yytemp[] with Ys with scrambled residuals */
      for (i=0; i<nn; i++)  
         yytemp[i] = beta*xx[i] + resid[nrand(nn)];
      /* Compute beta for (xx,yytemp) */
      zboot[ns] = gettheilbeta(xx,yytemp,nn,NULL); }
   /* Show bootstrap output */
   show_bootstrap_output("BootResid", zboot,nboot, beta);
      }



void bootstrap_beta_rank_regression (const char *Ytype,
     const double *xx, const double *yy, int nn)
 { int i,ns;
   double beta, resid[MAXN], xxtemp[MAXN],yytemp[MAXN];
   double zboot[NBOOT+3];
   printf ("\n(RANK REGRESSION) BOOTSTRAP for beta in %s = beta*X + mu:\n", 
     Ytype);
   printf ("  Medians and 95%% confidence intervals, %s replications.\n",
      commnum(nboot));
   /* Write a note to the summary list */
   addrcomm("Rank Regression", Ytype);
   /* Estimate beta for (xx,yy) */
   /* NULL means that we aren't interested in mu */
   beta=getrankbeta(xx,yy,nn, NULL);
   printf ("BOOTSTRAP ON OBSERVATIONS (rank regression, beta=%g):\n", beta);
   for (ns=0; ns<nboot; ns++)
    { for (i=0; i<nn; i++)
       { int ran=nrand(nn);
         xxtemp[i]=xx[ran];  yytemp[i]=yy[ran]; }
      zboot[ns] = getrankbeta(xxtemp,yytemp,nn,NULL);  }
   /* Show bootstrap output */
   show_bootstrap_output("BootObs", zboot,nboot, beta);
   
   printf ("\nBOOTSTRAP ON RESIDUALS (rank regression, beta=%g):\n", beta);
   /* Put residuals in resid[] */
   for (i=0; i<nn; i++)  resid[i] = yy[i]-beta*xx[i];
   /* Bootstrap of residuals */
   for (ns=0; ns<nboot; ns++)
    { /* Replace yytemp[] with Ys with scrambled residuals */
      for (i=0; i<nn; i++)  
         yytemp[i] = beta*xx[i] + resid[nrand(nn)];
      /* Compute beta for (xx,yytemp) */
      zboot[ns] = getrankbeta(xx,yytemp,nn,NULL); }
   /* Show bootstrap output */
   show_bootstrap_output("BootResid", zboot,nboot, beta);
      }



/* Routines to display a summary of CIs and P-values */

/* Set aside a buffer of MAXLINES lines */
/* Each summlines[i] is a string buffer of length 200 */

#define MAXLINES 50

char summlines[MAXLINES][200];
int nsumm;    /* Number of summary lines in use */

void addrcomm (const char *RegName, const char *Ytype)
 { /* Write title to a line and increase the count */
   sprintf(summlines[nsumm++], "%s (%s):", RegName, Ytype); }

void addrsumm (const char *RegName, double beta,
   double cilow, double cihigh, double pval, int df)
 { char dbuff[55];
   /* Is there an entry for degrees of freedom? */
   /* -1 means no entry */
   if (df>=0) sprintf(dbuff, " (df=%d)", df);
   else dbuff[0]=0;
   /* Write data to a summary line and increase the count */
   sprintf(summlines[nsumm++], "%-12s  %9.5f  (%9.5f, %9.5f)  P=%7.5f%s",
     RegName, beta, cilow,cihigh, pval, dbuff); }

/* Display the summary lines */

void Dump_Regression_Summary(void)
 { int i;
   for (i=0; i<nsumm; i++)
     printf ("%s\n", summlines[i]); }