/* Extended treatment of Miller jackknife procedure in the text.
   Test whether the variances of two samples are significantly 
     different, not assuming that the samples are normally distributed
   Uses a stratified jackknife (that is, within each sample)
     for the sample log-standard-deviation.
   Compare with a stratified bootstrap */

#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"

/* Tritiated water diffusion data across human chorioamnion */
/* Table 4.1 in Hollander & Wolfe */
/* xx[] values are at term.  yy[] values are at 12-26 weeks. */

const char *title_data=
  "Diffusion of tritiated water (that is, water containing H^3)\n"
  "  across human chorioamnion as a test of placental permeability\n"
  "The first sample is at term (approximately 39 weeks)\n"
  "  The second sample is at 12-26 weeks.\n"
  "Table 4.1 (p110) in Hollander & Wolfe, 2nd 3dn";

#define MAXN 30        /* Maximum possible sample size */
#define NBOOT 10000    /* Do NBOOT resamples in stratified bootstrap */

/* Sample at term (control sample): */

const char *xname="At term (approx 39 weeks)";
const int mm=10;
const double xx[MAXN] = { 0.80, 0.83, 1.89, 1.04, 1.45, 
   1.38, 1.91, 1.64, 0.73, 1.46 };

/* Sample at 12-26 weeks (test sample): */

const char *yname="At 12-26 weeks";
const int nn=5;
const double yy[MAXN] = { 1.15, 0.88, 0.90, 0.74, 1.21 };

/* Number of bootstrap replications */

int nboot=NBOOT;
unsigned long startseed=123456;


/* Subroutines to keep us from repeating standard blocks of code. */

/* The keyword `const' means that this function promises not to change */
/*   any of the values pointed to. (If it tries to, the compiler will */
/*   crash.)  This means that there is one less thing for the calling */
/*   function to worry about. */

void do_miller_jackknife(const double *xx, int mm,
   const double *yy, int nn, double gamma_obs);

void do_bootstrap(const double *xx, int mm,
   const double *yy, int nn, double gamma_obs);

int main(int argc, char **argv)
 { int i,j;
   double xmean,xvar,xstdev, ymean,yvar,ystdev;
   double gamma_obs;
   unsigned long wseed;
   printf (
     "Use the Miller method to test whether or not the VARIANCES of\n"
     "  two samples are the same (see Section 5.2 in text).\n"
     "Miller's idea is to apply a `stratified' jackknife to the\n"
     "  log sample variance.\n"
     "Compare Miller's P-value and 95%% confidence interval with a\n"
     "  stratified bootstrap P-value and 95%% confidence interval.\n");
   /* Title data */
   printf ("\n%s\n", title_data);
   /* Note: If ``var'' is any variable, ``&var'' is a POINTER to that */
   /*   variable, or equivalently the ADDRESS of that variable. */
   /* A function in C cannot change the value of a variable in the */
   /*   calling function unless it is given its address in memory. */
   /* The function  getmeanvar()  fills  xmean,xvar  with the */
   /*   sample mean and sample variance of the array xx[]. */
   /* getmeanvar() is in statlib.c */
   getmeanvar(xx,mm, &xmean,&xvar);  xstdev=sqrt(xvar);
   printf ("\n%s (Xs, m=%d)\n", xname,mm);
   printf ("(Xbar=%g, SampStdev=%.3f, Sampvar=%.3f)\n",
     xmean,xstdev,xvar);
   for (i=0; i<mm; i++) printf ("  %.2f",xx[i]);  printf("\n");
   /* Now the second sample */
   getmeanvar(yy,nn, &ymean,&yvar);  ystdev=sqrt(yvar);
   printf ("\n%s (Ys, n=%d)\n", yname,nn);
   printf ("(Ybar=%g, SampStdev=%.3f, Sampvar=%.3f)\n",
     ymean,ystdev,yvar);
   for (j=0; j<nn; j++) printf ("  %.2f",yy[j]);  printf("\n");

   printf ("\nSet gamma = SampVar(X)/SampVar(Y)\n");
   gamma_obs=xvar/yvar;
   printf ("Observed gamma = %g/%g = %g\n", xvar,yvar, gamma_obs);

   printf ("\nJACKKNIFE PROCEDURE:\n");
   /* Carry out the algorithm described in Section 5.2 of the text */
   do_miller_jackknife(xx,mm, yy,nn, gamma_obs);

   /* Initialize the randum-number generator for bootstrap replications */
   /* 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  will be used. */
   if (argc>=2) startseed=atol(argv[1]);
   wseed=setrand(startseed);
   printf ("\nStarting random-number seed: %lu\n",wseed);

   printf ("\nBOOTSTRAP PROCEDURE (n=%s bootstrap replications):\n",
     commnum(nboot));
   do_bootstrap(xx,mm, yy,nn, gamma_obs);

   show_randnums_used();
   /* Remind the user */
   printf ("Starting random-number seed: %lu\n",wseed);

   return 0; }


/* ``Getlogvar'' returns the log of the sample variance */
/* This is the estimator that is being jackknifed and bootstrapped: */
/* Again, `const' means that the subroutine promises not to change */
/*   any of the values in the array. */

double getlogvar(const double *xx, int num)
 { int i;  double dsum,dmean, dvar;
   for (dsum=0.0, i=0; i<num; i++) dsum+=xx[i];
   dmean=dsum/num;
   for (dsum=0.0, i=0; i<num; i++) dsum+=(xx[i]-dmean)*(xx[i]-dmean);
   dvar=dsum/(num-1);
   return log(dvar);  }


/* The following two routines carry out the `Miller jackknife' */
/*   procedure in Section 5.2 of the text (Hollander+Wolfe) */

/* The Miller jackknife procedure needs pseudovalues for each sample, */
/*   calculated independently for each sample. (That is, these are */
/*   `stratified' pseudovalues.) */
/* To avoid repeating the same code, we use a subroutine `do_jackknife' */
/*   that fills an array  ps[]  with pseudovalues calculated from the */
/*   array  dat[]. */
/* The key word `const' makes it clear that `do_jackknife' reads values */
/*   from dat[] and writes values to ps[] and not vice versa. */


/* General bookstrap recipe:  Compute pseudovalues for the sample xx[] and put */
/*   the pseudovalues in psval[] */

void do_jackknife(const char *jname, double *ps, const char *sampname, 
    const double *xx, int num)
 { int i,k,r;  double gfull;
   /* Get the estimate for the full sample: */
   gfull=getlogvar(xx,num);
   printf ("%s of %s (n=%d) for full sample: %g\n", 
     jname, sampname, num, gfull);
   /* Compute and store `num' pseudovalues: */
   for (i=0; i<num; i++)
    { /* Build the i-th delete-1 resample in xsamp[]: */
      double xsamp[MAXN];
      for (r=k=0; k<num; k++) if (k!=i)  /* Drop the i-th observation */
        xsamp[r++]=xx[k];       /* End with r=num-1 */
      /* Store the i-th pseudovalue: */
      ps[i] = num*gfull - (num-1)*getlogvar(xsamp,num-1); }
   printf ("Pseudovalues of %s (%s, n=%d):\n", sampname, jname, num);
   for (i=0; i<num; i++)
    { if (i>0 && (i%5)==0) printf("\n");
      printf ("  %.4f", ps[i]); }
   printf("\n"); }




void do_miller_jackknife(const double *xx, int mm,
   const double *yy, int nn, double gamma_obs)
 { /* Two sets of pseudo-values */
   double aa[MAXN], bb[MAXN];
   double abar,bbar, savar,sbvar, qsvar,qsdev, qq,gamma;
   printf(
   "First, get the jackknife pseudovalues for the log sample variance\n"
   "  for each sample:\n");
   /* Collect the pseudovalues for xx[] into aa[] and yy[] into bb[], */
   /*   following the notation in the text */
   do_jackknife("LogSampVar",aa,"Sample 1",xx,mm);
   do_jackknife("LogSampVar",bb,"Sample 2",yy,nn);

   /* Find the means and sample variances of the two sets of pseudovalues */
   /* Remember that aa[],bb[] are pseudovalues of LOG sample variances, */
   /*   not of the sample variances themselves. */
   getmeanvar(aa,mm, &abar,&savar);
   getmeanvar(bb,nn, &bbar,&sbvar);
   /* The strategy is to use the two sets of pseudo-values as */
   /*   two independent samples with the sample variance */
   /* The unpooled variance estimator for Stdev(abar-bbar) is: */
   qsvar=(savar/mm)+(sbvar/nn);
   qsdev=sqrt(qsvar);
   printf (
   "\nJACKKNIFE TEST of H_0:gamma=1  for  gamma=Var(X)/Var(Y):\n"
   "Sample means and sample variances of the two samples of\n"
   "  pseudovalues of log sample variances:\n"
   "  Abar=%.3f  Bbar=%.3f    S2Avar=%.3f  S2Bvar=%.3f\n"
   "  Unpooled variance of Abar-Bbar=%.3f  Unpooled StDev=%.3f\n",
     abar,bbar, savar,sbvar, qsvar,qsdev);
   /* The two-sample Z statistic for pooled variances: */
   /* Recall that aa[],bb[] are log variances, not the variances */
   /*   themselves, so that gamma=1 corresponds to qq=0. */
   qq=(abar-bbar)/qsdev;
   printf ("V1=S2Avar/m (etc) and the Miller two-sample Z-statistic:\n");
   printf ("  V1=%.3f  V2=%.3f    Z=%.3f   P=%.4f (2-sided)\n",
     savar/mm, sbvar/nn, qq, normpval(qq));
   /* Convert `abar-bbar' to gamma=Var(X)/Var(Y) */
   gamma=exp(abar-bbar);
   printf (
   "\nFor gamma(X,Y)=s2(X)/s2(Y):\n"
   "Bias-corrected gamma estimator and symmetric 95%% CI for gamma:\n"
   "  %.4f   (%.4f, %.4f)   Observed gamma=%.4f\n",
     gamma, exp(abar-bbar-1.960*qsdev),
     exp(abar-bbar+1.960*qsdev), gamma_obs);
       }


void do_bootstrap(const double *xx, int mm,
   const double *yy, int nn, double gamma_obs)
 { int i,j,mb, ns, nbig, kk,klow,khigh;
   double gammboot[NBOOT], bootmed, bootpval, gamma0;
   printf ("Use a `stratified bootstrap' to generate an array\n");
   printf ("  of gamma(star) bootstrap resampled values for\n");
   printf ("  gamma(star) = s(X)^2(star)/s(Y)^2(star)\n");
   printf("Filling the bootstrap ``gamma(star)'' array:\n");
   for (ns=0; ns<nboot; ns++)
    { /* Arrays for two bootstrap resamples */
      double xtemp[MAXN], ytemp[MAXN];
      double savar,sbvar;  /* Store SampVar(xxstar),SampVar(yystar) here. */
      /* Define xx1=xx(star) = bootstrap resample of xx[] */
      for (i=0; i<mm; i++)  xtemp[i]=xx[nrand(mm)];
      /* Tell getmeanvar() we only want the variance */
      getmeanvar(xtemp,mm, NULL,&savar);
      do { /* Define ytemp=yy(star) = bootstrap resample of yy[] */
           /*   BUT IGNORE resamples with ytemp[j]=constant */
           /* This loop repeats until fabs(sbvar)>1e-8, which should be safe */
           for (j=0; j<nn; j++)  ytemp[j]=yy[nrand(nn)];
           getmeanvar(ytemp,nn, NULL,&sbvar); }
      while(fabs(sbvar)<1e-8);
      /* Store gamma(star)=SampVar(xtemp)/SampVar(ytemp) */
      gammboot[ns]=savar/sbvar; }
   /* Use the system function `qsort()' to sort the gamma(star) array */
   /*   so that we can display a bootstrap percentile CI. */
   /* See discussion in statlib.c (dubsort() is in statlib.c) */
   printf("Sorting the bootstrap ``gamma(star)'' array:\n");
   qsort(gammboot,nboot, sizeof(double),dubsort);
   /* Short first and last sampled gamma(star) values */
   printf("The first 10 gamma(star) values are:\n");
   for (i=0; i<10; i++)
    { if (i>0 && (i%5)==0) printf("\n");
      printf ("   %.3f",gammboot[i]); }  printf("\n");
   printf("The last 10 gamma(star) values are:\n");
   for (i=0; i<10; i++)
    { if (i>0 && (i%5)==0) printf("\n");
      printf ("   %.3f",gammboot[nboot-10+i]); }  printf("\n");

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

   /* The median bootstrap value */
   if (nboot%2) bootmed=gammboot[nboot/2];
   else  bootmed=(gammboot[nboot/2] + gammboot[(nboot-2)/2])/2;

   printf(
   "\nFor gamma(X,Y)=s2(X)/s2(Y):\n"
   "Bootstrap median and bootstrap 95%% CI for gamma:\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 gamma=%.4f\n",
     bootmed, gammboot[klow], gammboot[khigh], gamma_obs);

   printf ("\nBootstrap P-value for H_0:gamma=gamma0=1  "
     "(gamma(observed)=%g):\n",gamma_obs);
   gamma0=1;  nbig=0;
   if (gamma0<=bootmed)  /* gamma0=1 is on lower tail of distribution */
    { for (i=0; i<nboot; i++) if (gammboot[i]<=gamma0) nbig++; }
   else                  /* gamma0=1 is on upper tail */
    { for (i=0; i<nboot; i++) if (gammboot[i]>=gamma_obs) nbig++; }
   printf ("(Twice the proportion of resampled gammas %s %g.)\n",
     (gamma0<=bootmed) ? "<=" : ">=", gamma0);
   bootpval=(double)(2*nbig)/nboot;
   if (bootpval>1.0) bootpval=1.0;
   printf ("  P = 2*%d/%s = %.4f (2-sided)\n",
     nbig,commnum(nboot), bootpval);
        }