/* Jackknife/Bootstrap for the slope of a regression */

#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 MAXN 50
#define NBOOT 10000

/* Data for single covariate */

int nn=10;
double xx[MAXN];
double xbar,xsse;
double tcrit;  /* 95% critical value for t(nn-2) */

/* Constants for regressions with normal errors */
double nmu=3.0, nbeta=3.0, nsigma=20.0;

/* Constants for regressions with exponential errors */
double exmu=3.0, exbeta=3.0, exsigma=20.0;

/* Number of bootstrap replications and default starting RN seed */
/* Nboot is also the number of simulated CIs for coverage probabilities */
int nboot=NBOOT;
unsigned long startseed=123456;


/* Prototypes */
void show_ols_analysis(const char *name,
   const double *yy, const double *xx, int nn,
   double truemu, double truebeta, double truesigma);

void find_normal_ols_ci_coverage(int nsims,
     double beta, double mu, double sigma, int nn);
void find_exponential_ols_ci_coverage(int nsims,
     double beta, double mu, double sigma, int nn);

void do_jackknife(const char *name,
     const double *yy, const double *xx, int nn);
void find_normal_jackknife_ci_coverage(int nsims,
     double beta, double mu, double sigma, int nn);
void find_exponential_jackknife_ci_coverage(int nsims,
     double beta, double mu, double sigma, int nn);

void do_bootstrap_record(const char *name,
     const double *yy, const double *xx, int nn);
void find_normal_bootstrap_record_ci_coverage(int nsims,
     double beta, double mu, double sigma, int nn);
void find_exponential_bootstrap_record_ci_coverage(int nsims,
     double beta, double mu, double sigma, int nn);

void do_bootstrap_residual(const char *name,
     const double *yy, const double *xx, int nn);
void find_normal_bootstrap_residual_ci_coverage(int nsims,
     double beta, double mu, double sigma, int nn);
void find_exponential_bootstrap_residual_ci_coverage(int nsims,
     double beta, double mu, double sigma, int nn);


/* For a final table of confidence intervals */
/*   and a final table of coverage probabilities */
void load_ci_values(const char *name, 
  double med, double cilow, double cihigh);
void show_ci_table(void);

void load_ci_coverage_values(const char *name, int ncover);
void show_ci_coverage_table(void);


int main(int argc, char **argv)
 { int i;
   double xsum, ysum,ybar, yy[MAXN];
   unsigned long wseed;
   /* Helpful for segmentation faults and other fatal errors */
   /* This function is in statlib.c */
   set_newsignals();
   printf (
   "\nCompare Jackknife/Bootstrap methods for the slope of a regression\n");

   /* Initialize the randum-number generator: */
   /* Allow a new seed (or 0) to be entered on the command line. */
   /* Setrand(0) initializes the seed from the system clock. */
   /* If nothing is entered on the command line, the default value */
   /*   of  startseed  is used. */
   if (argc>=2) startseed=atol(argv[1]);
   wseed=setrand(startseed);
   printf ("Starting random-number seed: %lu\n",wseed);

   for (i=0; i<nn; i++)
    { double dd=drand();  xx[i]=floor(300*dd*dd+1); }
   /* System `quick sort' of array xx[] of doubles */
   /* dubsort() is in statlib.c */
   /* See statlib.c for a discussion of qsort() and comparison */
   /*   functions such as dubsort() */
   qsort(xx,nn, sizeof(double),dubsort);

   /* Compute global variables xbar and xsse */
   for (xsum=0.0, i=0; i<nn; i++) xsum+=xx[i];
   xbar=xsum/nn;
   for (xsse=0.0, i=0; i<nn; i++) xsse+=(xx[i]-xbar)*(xx[i]-xbar);
   printf ("\nX-coordinate of regressions (n=%d, xbar=%.3f, xsse=%g):\n",
     nn,xbar,xsse);
   for (i=0; i<nn; i++)
    { if (i>0 && (i%15)==0) printf("\n");
      printf ("  %g",xx[i]); }  printf("\n");
   /* 95% critical value for t(nn-2): crit_tdist() is in statlib.c */
   /* This value is needed for Student-t confidence intervals */
   tcrit=crit_tdist(nn-2,0.05);

   printf (
   "\nREGRESSION WITH NORMAL ERRORS (n=%d, mu=%g, beta=%g, sigma=%g)\n",
     nn, nmu, nbeta, nsigma);
   for (i=0; i<nn; i++)
     yy[i] = nmu + nbeta*xx[i] + nsigma*normrand();
   for (ysum=0.0, i=0; i<nn; i++)  ysum+=yy[i];
   ybar=ysum/nn;
   printf ("Normal regression data (n=%d, ybar=%.3f, xbar=%.3f):\n", 
     nn,ybar,xbar);
   for (i=0; i<nn; i++)
    { if (i>0 && (i%4)==0) printf("\n");
      printf ("   (%g, %.2f)",xx[i],yy[i]); }  printf("\n");
   show_ols_analysis("Normal", yy,xx,nn, nmu,nbeta,nsigma);
   find_normal_ols_ci_coverage(nboot, nbeta,nmu,nsigma, nn);

   printf ("\nJACKKNIFE ESTIMATES OF BETA FOR NORMAL DATA:\n");
   do_jackknife("Normal", yy,xx,nn);
   find_normal_jackknife_ci_coverage(nboot, nbeta,nmu,nsigma, nn);

   printf ("\nBOOKSTRAP RECORD ESTIMATES OF BETA FOR NORMAL DATA:\n");
   do_bootstrap_record("Normal", yy,xx,nn);
   find_normal_bootstrap_record_ci_coverage(nboot, nbeta,nmu,nsigma, nn);

   printf ("\nBOOKSTRAP RESIDUAL ESTIMATES OF BETA FOR NORMAL DATA:\n");
   do_bootstrap_residual("Normal", yy,xx,nn);
   find_normal_bootstrap_residual_ci_coverage(nboot, nbeta,nmu,nsigma, nn);


   printf (
   "\nREGRESSION WITH EXPONENTIAL ERRORS (n=%d, mu=%g, beta=%g, sigma=%g)\n",
     nn, exmu, exbeta, exsigma);
   for (i=0; i<nn; i++)
     yy[i] = exmu + exbeta*xx[i] + exsigma*exprand();
   for (ysum=0.0, i=0; i<nn; i++)  ysum+=yy[i];
   ybar=ysum/nn;
   printf ("Exponential regression data (n=%d, ybar=%.3f, xbar=%.3f):\n", 
     nn,ybar,xbar);
   for (i=0; i<nn; i++)
    { if (i>0 && (i%4)==0) printf("\n");
      printf ("   (%g, %.2f)",xx[i],yy[i]); }
   show_ols_analysis("Exponential", yy,xx,nn, nmu,nbeta,nsigma);
   find_exponential_ols_ci_coverage(nboot, nbeta,nmu,nsigma, nn);

   printf ("\nJACKKNIFE ESTIMATES OF BETA FOR EXPONENTIAL DATA:\n");
   do_jackknife("Exponential", yy,xx,nn);
   find_exponential_jackknife_ci_coverage(nboot, exbeta,exmu,exsigma, nn);

   printf ("\nBOOKSTRAP RECORD ESTIMATES OF BETA FOR EXPONENTIAL DATA:\n");
   do_bootstrap_record("Exponential", yy,xx,nn);
   find_exponential_bootstrap_record_ci_coverage(nboot,
      exbeta,exmu,exsigma, nn);

   printf ("\nBOOKSTRAP RESIDUAL ESTIMATES OF BETA FOR EXPONENTIAL DATA:\n");
   do_bootstrap_residual("Exponential", yy,xx,nn);
   find_exponential_bootstrap_residual_ci_coverage(nboot,
      exbeta,exmu,exsigma, nn);

   show_ci_table();
   show_ci_coverage_table();

   show_randnums_used();
   printf ("Initial random-number seed: %lu\n",wseed);

   return 0; }



double get_mean(const double *xx, int nn)
 { int i;  double xsum;
   for (xsum=0.0, i=0; i<nn; i++) xsum+=xx[i];
   return xsum/nn; }


double get_betahat(const double *yy, 
     const double *xx, int nn, double *muhatptr)
 { int i;
   double ybar,xbar, betasum,xsse, betahat;
   ybar=get_mean(yy,nn);  xbar=get_mean(xx,nn);
   betasum=xsse=0.0;
   for (i=0; i<nn; i++)
    { betasum+=(yy[i]-ybar)*(xx[i]-xbar);
      xsse+=(xx[i]-ybar)*(xx[i]-xbar); }
   betahat=betasum/xsse;
   if (muhatptr!=NULL) *muhatptr=ybar-betahat*xbar;
   return betahat; }


double get_betahat_short(const double *yy, int nn, double *muhatptr)
 { int i;
   double ysum,ybar, betasum,xsse, betahat;
   for (ysum=0.0, i=0; i<nn; i++) ysum+=yy[i];
   ybar=ysum/nn;  betasum=xsse=0.0;
   for (i=0; i<nn; i++)
    { betasum+=(yy[i]-ybar)*(xx[i]-xbar);
      xsse+=(xx[i]-ybar)*(xx[i]-xbar); }
   betahat=betasum/xsse;
   if (muhatptr!=NULL) *muhatptr=ybar-betahat*xbar;
   return betahat; }


void show_ols_analysis(const char *name,
   const double *yy, const double *xx, int nn,
   double truemu, double truebeta, double truesigma)
 { char nbuff[99];  const char *rname;
   int i;
   double muhat,betahat, rsse,rmse,betavarhat;
   double bwid, cilow,cihigh;
   betahat=get_betahat(yy,xx,nn, &muhat);
   printf ("\nLeast-squares estimates:  BetaHat=%.3f   Muhat=%.3f\n",
      betahat, muhat);
   printf ("  True values:            Beta=%.3f      Mu=%.3f \n",
      truebeta, truemu);
   for (rsse=0.0, i=0; i<nn; i++)
    { double resid=yy[i]-muhat-betahat*xx[i];  rsse+=resid*resid; }
   rmse=rsse/(nn-2);  betavarhat=rmse/xsse;
   printf ("  Residual MSE=%g  (true sigma^2=%g)\n",
     rmse,truesigma*truesigma);
   printf ("  Standard Est(Var(betahat)):  %g  True Var(betahat): %g\n",
      betavarhat, truesigma*truesigma/xsse);
   bwid=tcrit*sqrt(betavarhat);
   cilow=betahat-bwid;  cihigh=betahat+bwid;
   printf ("  t(%d,0.975)=%g   (z(0.975)=1.960)   CIwid=%.3f\n", nn-2,tcrit,bwid);
   printf ("\nFull-sample estimates for true beta:\n");
   printf ("  BetaHat=%.3f   95%% CI (%.3f, %.3f)   True beta=%.3f\n",
     betahat, cilow, cihigh, truebeta);
   sprintf (nbuff, "%s, Classical t(%d)", name,nn-2);
   /* Make a non-local copy of the string */
   rname=strdup(nbuff);
   /* Save values for final summary table */
   load_ci_values(rname, betahat, cilow, cihigh);
     }



int get_beta_ci_ols_cover(const double *yy, int nn, double truebeta)
 { int i;
   double muhat,betahat, rsse,rmse,betavarhat;
   double bwid, cilow,cihigh;
   betahat=get_betahat_short(yy,nn,&muhat);
   for (rsse=0.0, i=0; i<nn; i++)
    { double res=yy[i]-muhat-betahat*xx[i];  rsse+=res*res; }
   rmse=rsse/(nn-2);  betavarhat=rmse/xsse;
   bwid=tcrit*sqrt(betavarhat);
   cilow=betahat-bwid;
   cihigh=betahat+bwid;
   if (cilow<truebeta && truebeta<cihigh) return TRUE;
   else return FALSE;
     }


void find_normal_ols_ci_coverage(int nsims,
     double beta, double mu, double sigma, int nn)
 { int i,ns,ncover;
   printf ("\nCoverage probability of beta for "
     "95%% normal CIs (true beta=%g):\n", beta);
   printf ("Normal samples of size %d:\n",nn);
   for (ncover=ns=0;  ns<nsims;  ns++)
    { double yy[MAXN];
      for (i=0; i<nn; i++)  yy[i] = mu + beta*xx[i] + sigma*normrand();
      if (get_beta_ci_ols_cover(yy,nn,beta)) ncover++; }
   printf ("Estimated coverage = %s/%s = %.2f%%\n",
     commnum1(ncover),commnum2(nsims), (100.0*ncover)/nsims);
   load_ci_coverage_values("Normal, Classical OLS", ncover); }


void find_exponential_ols_ci_coverage(int nsims,
     double beta, double mu, double sigma, int nn)
 { int i,ns,ncover;
   printf ("\nCoverage probability of beta for "
     "95%% normal CIs (true beta=%g):\n", beta);
   printf ("Exponential samples of size %d:\n",nn);
   for (ncover=ns=0;  ns<nsims;  ns++)
    { double yy[MAXN];
      for (i=0; i<nn; i++)  yy[i] = mu + beta*xx[i] + sigma*exprand();
      if (get_beta_ci_ols_cover(yy,nn,beta)) ncover++; }
   printf ("Estimated coverage = %s/%s = %.2f%%\n",
     commnum1(ncover),commnum2(nsims), (100.0*ncover)/nsims);
   load_ci_coverage_values("Exponential, Classical OLS", ncover); }



void do_jackknife(const char *name, const double *yy, const double *xx, int nn)
 { char nbuff[88];  const char *rname;
   int i,j,r;
   double obs_beta, psval[MAXN], jmean,jvar,jwid;
   obs_beta=get_betahat(yy,xx,nn,NULL);
   for (i=0; i<nn; i++)
    { double ydel[MAXN], xdel[MAXN];
      for (r=j=0; j<nn; j++) if (j!=i)
       { ydel[r]=yy[j];  xdel[r]=xx[j];  r++; }
      psval[i] = nn*obs_beta - (nn-1)*get_betahat(ydel,xdel,nn-1,NULL); }
   getmeanvar(psval,nn, &jmean,&jvar);
   jwid=1.960*sqrt(jvar);
   printf ("Jackknife record mean: %.3f   95%% CI: (%.3f, %.3f)\n",
     jmean, jmean-jwid, jmean+jwid);
   sprintf(nbuff, "%s, Jackknife", name);
   /* Make a non-local copy of the string */
   rname=strdup(nbuff);
   load_ci_values(rname, jmean,jmean-jwid,jmean+jwid); }



int get_beta_ci_jackknife_cover(const double *yy, int nn, double truebeta)
 { int i,j,r;
   double obs_beta, psval[MAXN], jmean,jvar,jwid, cilow,cihigh;
   obs_beta=get_betahat_short(yy,nn,NULL);
   for (i=0; i<nn; i++)
    { double ydel[MAXN], xdel[MAXN];
      for (r=j=0; j<nn; j++) if (j!=i)
       { ydel[r]=yy[j];  xdel[r]=xx[j];  r++; }
      psval[i] = nn*obs_beta - (nn-1)*get_betahat(ydel,xdel,nn-1,NULL); }
   getmeanvar(psval,nn, &jmean,&jvar);
   jwid=1.960*sqrt(jvar);
   cilow=jmean-jwid;  cihigh=jmean+jwid;
   if (cilow<truebeta && truebeta<cihigh) return TRUE;
   else return FALSE;
     }


void find_normal_jackknife_ci_coverage(int nsims,
     double beta, double mu, double sigma, int nn)
 { int i,ns,ncover;
   printf ("\nCoverage probability of beta for "
     "95%% jackknife CIs (true beta=%g):\n", beta);
   printf ("Normal samples of size %d:\n",nn);
   for (ncover=ns=0;  ns<nsims;  ns++)
    { double yy[MAXN];
      for (i=0; i<nn; i++)  yy[i] = mu + beta*xx[i] + sigma*normrand();
      if (get_beta_ci_jackknife_cover(yy,nn,beta)) ncover++; }
   printf ("Estimated coverage = %s/%s = %.2f%%\n",
     commnum1(ncover),commnum2(nsims), (100.0*ncover)/nsims);
   load_ci_coverage_values("Normal, Jackknife", ncover); }


void find_exponential_jackknife_ci_coverage(int nsims,
     double beta, double mu, double sigma, int nn)
 { int i,ns,ncover;
   printf ("\nCoverage probability of beta for "
     "95%% jackknife CIs (true beta=%g):\n", beta);
   printf ("Exponential samples of size %d:\n",nn);
   for (ncover=ns=0;  ns<nsims;  ns++)
    { double yy[MAXN];
      for (i=0; i<nn; i++)  yy[i] = mu + beta*xx[i] + sigma*exprand();
      if (get_beta_ci_jackknife_cover(yy,nn,beta)) ncover++; }
   printf ("Estimated coverage = %s/%s = %.2f%%\n",
     commnum1(ncover),commnum2(nsims), (100.0*ncover)/nsims);
   load_ci_coverage_values("Exponential, Jackknife", ncover); }







void do_bootstrap(double *betastar, int nboot, double obs_beta,
    const char *name, const char *boottype)
 { char nbuff[88];  const char *rname;
   int i,mb, kk,klow,khigh;
   double bootmed;
   printf("Sorting the bootstrap ``beta(star)'' array:\n");
   /* System sort of double array */
   /* See comments in statlib.c around `dubsort()' */
   qsort(betastar,nboot, sizeof(double),dubsort);

   /* Calculate the bootstrap median bias */
   for (mb=i=0; i<nboot; i++) if (betastar[i]<=obs_beta) mb++;
   printf ("Bootstrap median bias (Obs_Beta=%g):  %d/%d = %.3f\n",
     obs_beta, mb,nboot, (double)mb/nboot);

   /* The median bootstrap value */
   /* getmedian_sorted() is in statlib.c */
   bootmed=getmedian_sorted(betastar,nboot);
   printf(
   "\nBootstrap median and bootstrap 95%% CI for beta:\n");
   /* Assume that 0.025*nboot is an integer */
   kk=klow=0.025*nboot;
   if (kk>0) klow=kk-1;    /* k-th from the bottom */
   khigh=nboot-1-klow;     /* k-th from the top */
   printf (
   "  %.4f   (%.4f, %.4f)   Observed beta=%.4f\n",
     bootmed, betastar[klow], betastar[khigh], obs_beta);
   sprintf(nbuff, "%s, Bootstrap %s", name, boottype);
   /* Make a non-local copy of the string */
   rname=strdup(nbuff);
   load_ci_values(rname, bootmed, betastar[klow], betastar[khigh]);
        }




void do_bootstrap_record(const char *name, const double *yy, const double *xx, int nn)
 { int i,ns;
   double obs_beta, betastar[NBOOT];
   obs_beta=get_betahat(yy,xx,nn,NULL);
   printf("Filling the bootstrap ``beta(star)'' array:\n");
   for (ns=0; ns<nboot; ns++)
    { double ystar[MAXN], xstar[MAXN];
      for (i=0; i<nn; i++)
       { int k=nrand(nn);
         ystar[i]=yy[k];  xstar[i]=xx[k]; }
      betastar[ns]=get_betahat(ystar,xstar,nn,NULL); }
   /* For safety, do the rest of the analysis in a common subroutine */
   do_bootstrap(betastar,nboot, obs_beta, name,"record"); }


void do_bootstrap_residual(const char *name, const double *yy, const double *xx, int nn)
 { int i,ns;
   double obs_beta, obs_mu;
   double resid[MAXN], betastar[NBOOT];
   /* Find the residuals */
   obs_beta=get_betahat(yy,xx,nn,&obs_mu);
   printf ("Residuals of  Y = %.3f*X + %.3f:\n", obs_beta,obs_mu);
   for (i=0; i<nn; i++)  resid[i] = yy[i] - obs_mu - (obs_beta)*xx[i];
   for (i=0; i<nn; i++)
    { if (i>0 && (i%5)==0) printf("\n");
      printf ("  %.3f",resid[i]); }  printf("\n");

   printf("Filling the bootstrap ``beta(star)'' array:\n");
   for (ns=0; ns<nboot; ns++)
    { double ystar[MAXN];
      for (i=0; i<nn; i++)
       { double residstar=resid[nrand(nn)];
         ystar[i] = obs_mu + (obs_beta)*xx[i] + residstar; }
      betastar[ns]=get_betahat(ystar,xx,nn,NULL); }
   /* For safety, do the rest of the analysis in a common subroutine */
   do_bootstrap(betastar,nboot, obs_beta, name,"residual"); }



int get_beta_ci_bootstrap_cover(double *betastar, int nboot,
      double truebeta)
 { int kk,klow,khigh;
   qsort(betastar,nboot, sizeof(double),dubsort);
   kk=klow=0.025*nboot;
   if (kk>0) klow=kk-1;    /* k-th from the bottom */
   khigh=nboot-1-klow;     /* k-th from the top */
   if (betastar[klow]<truebeta && truebeta<betastar[khigh]) return TRUE;
   else return FALSE; }



int get_beta_ci_bootstrap_record_cover(const double *yy,
     int nn, double truebeta)
 { int i,ns;
   double betastar[NBOOT];
   for (ns=0; ns<nboot; ns++)
    { double ystar[MAXN], xstar[MAXN];
      for (i=0; i<nn; i++)
       { int k=nrand(nn);
         ystar[i]=yy[k];  xstar[i]=xx[k]; }
      betastar[ns]=get_betahat(ystar,xstar,nn,NULL); }
   return get_beta_ci_bootstrap_cover(betastar,nboot,truebeta); }


void find_normal_bootstrap_record_ci_coverage(int nsims,
     double beta, double mu, double sigma, int nn)
 { int i,ns,ncover;
   printf ("\nCoverage probability of beta for "
     "95%% bootstrap record CIs (true beta=%g):\n", beta);
   printf ("Normal samples of size %d:\n",nn);
   for (ncover=ns=0;  ns<nsims;  ns++)
    { double yy[MAXN];
      for (i=0; i<nn; i++)  yy[i] = mu + beta*xx[i] + sigma*normrand();
      if (get_beta_ci_bootstrap_record_cover(yy,nn,beta)) ncover++; }
   printf ("Estimated coverage = %s/%s = %.2f%%\n",
     commnum1(ncover),commnum2(nsims), (100.0*ncover)/nsims);
   load_ci_coverage_values("Normal, Bootstrap Record", ncover); }


void find_exponential_bootstrap_record_ci_coverage(int nsims,
     double beta, double mu, double sigma, int nn)
 { int i,ns,ncover;
   printf ("\nCoverage probability of beta for "
     "95%% bootstrap record CIs (true beta=%g):\n", beta);
   printf ("Exponential samples of size %d:\n",nn);
   for (ncover=ns=0;  ns<nsims;  ns++)
    { double yy[MAXN];
      for (i=0; i<nn; i++)  yy[i] = mu + beta*xx[i] + sigma*exprand();
      if (get_beta_ci_bootstrap_record_cover(yy,nn,beta)) ncover++; }
   printf ("Estimated coverage = %s/%s = %.2f%%\n",
     commnum1(ncover),commnum2(nsims), (100.0*ncover)/nsims);
   load_ci_coverage_values("Exponential, Boostrap Record", ncover); }



int get_beta_ci_bootstrap_residual_cover(const double *yy,
     int nn, double truebeta)
 { int i,ns;
   double obs_beta, obs_mu;
   double resid[MAXN], betastar[NBOOT];
   /* Find the residuals */
   obs_beta=get_betahat_short(yy,nn,&obs_mu);
   for (i=0; i<nn; i++)  resid[i] = yy[i] - obs_mu - (obs_beta)*xx[i];
   /* Find the bootstrap replications */
   for (ns=0; ns<nboot; ns++)
    { double ystar[MAXN];
      for (i=0; i<nn; i++)
       { double residstar=resid[nrand(nn)];
         ystar[i] = obs_mu + (obs_beta)*xx[i] + residstar; }
      betastar[ns]=get_betahat_short(ystar,nn,NULL); }
   return get_beta_ci_bootstrap_cover(betastar,nboot,truebeta); }


void find_normal_bootstrap_residual_ci_coverage(int nsims,
     double beta, double mu, double sigma, int nn)
 { int i,ns,ncover;
   printf ("\nCoverage probability of beta for "
     "95%% bootstrap residual CIs (true beta=%g):\n", beta);
   printf ("Normal samples of size %d:\n",nn);
   for (ncover=ns=0;  ns<nsims;  ns++)
    { double yy[MAXN];
      for (i=0; i<nn; i++)  yy[i] = mu + beta*xx[i] + sigma*normrand();
      if (get_beta_ci_bootstrap_residual_cover(yy,nn,beta)) ncover++; }
   printf ("Estimated coverage = %s/%s = %.2f%%\n",
     commnum1(ncover),commnum2(nsims), (100.0*ncover)/nsims);
   load_ci_coverage_values("Exponential, Boostrap Residual", ncover); }


void find_exponential_bootstrap_residual_ci_coverage(int nsims,
     double beta, double mu, double sigma, int nn)
 { int i,ns,ncover;
   printf ("\nCoverage probability of beta for "
     "95%% bootstrap residual CIs (true beta=%g):\n", beta);
   printf ("Exponential samples of size %d:\n",nn);
   for (ncover=ns=0;  ns<nsims;  ns++)
    { double yy[MAXN];
      for (i=0; i<nn; i++)  yy[i] = mu + beta*xx[i] + sigma*exprand();
      if (get_beta_ci_bootstrap_residual_cover(yy,nn,beta)) ncover++; }
   printf ("Estimated coverage = %s/%s = %.2f%%\n",
     commnum1(ncover),commnum2(nsims), (100.0*ncover)/nsims);
   load_ci_coverage_values("Normal, Boostrap Residual", ncover); }




/* For a final summary table of confidence intervals */

#define RR 23

const char *rname[RR];
double rmed[RR], rcilow[RR], rcihigh[RR];
int nr;

void load_ci_values(const char *name, 
  double med, double cilow, double cihigh)
 { rname[nr]=name;  rmed[nr]=med;  rcilow[nr]=cilow;  rcihigh[nr]=cihigh;
   nr++; }

void show_ci_table(void)
 { int i;  size_t len, mxwid;
   printf (
   "\nSummary of estimates and 95%% confidence intervals:\n"
   "True values of beta:  Normal: %g   Exponential: %g\n", nbeta, exbeta);
   for (mxwid=0, i=0; i<nr; i++)
     if ( (len=strlen(rname[i]))>mxwid) mxwid=len;
   for (i=0; i<nr; i++)
     printf ("  %-*s  %7.4f   (%7.4f, %7.4f)\n",
       mxwid,rname[i], rmed[i], rcilow[i], rcihigh[i]); }


/* For a final summary table of coverage probabilities */

const char *rcname[RR];
int rcover[RR];
int nc;

void load_ci_coverage_values(const char *name, int ncover)
 { rcname[nr]=name;  rcover[nc]=ncover;  nc++; }

void show_ci_coverage_table(void)
 { int i;  size_t len, mxwid;
   printf (
   "\nSummary of coverage probabilities of 95%% confidence intervals:\n"
   "True values of beta:  Normal: %g   Exponential: %g\n", nbeta, exbeta);
   for (mxwid=0, i=0; i<nr; i++)
     if ( (len=strlen(rname[i]))>mxwid) mxwid=len;
   for (i=0; i<nr; i++)
     printf ("  %-*s  %d/%s = %.2f%%\n",
       mxwid,rname[i], rcover[i],commnum(nboot), (100.0*rcover[i])/nboot); }