/* How can you estimate a variance, and also get a confidence interval
     for the variance, if the observations are not normally distributed? */

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
/* 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"

/* 
   The `jackknife' and the `bootstrap' are two general methods for obtaining
     a confidence interval for nearly any parameter based on nearly any
     estimator, with an attempt to obtain a bias-corrected estimate
     of the parameter as a bonus.

   This program illustrates both methods for finding a confidence
     interval for the standard deviation of a sample that happens to be
     from a normal distribution. In that case, the sample variance is an
     unbiased estimator of the true variance and there is a classical
     symmetric confidence interval for the standard deviation based on
     quantiles of the chi-square distribution. Hence we will have something
     to compare the `jackknife' and `bootstrap' estimators with.

   We apply the `jackknife' method to both the standard deviation itself
     and to the log standard deviation, not forgetting to exponentiate
     the jackknifed interval for the log standard deviation afterwards.
     The `log standard deviation' method is to provide some protection
     against outliers. The (percentile) bootstrap is `scale free' or
     `fully nonparametric' and gives the same results before and after
     any monotonic one-one transformation of the data. This will give us
     four confidence intervals to compare: The classical, two jackknifes,
     and the bootstrap.

   In this case, the four confidence intervals turn out to be similar,
     although the `log standard deviation' jackknife seems to do the best
     of the non-exact methods. However, the `jackknife' and `bootstrap'
     can be applied to nearly any functional, while even for the standard
     deviation the classical method is exact only for normally-distributed
     random variables.

   Enter the program name followed by `0' or an explicit seed for a
     different `original' sample to analyze.  Note that while the sample
     standard deviation is highly variable, the four estimates and the
     four confidence intervals track one another very closely. */

/* Default values for the sample size, the mean and standard deviation */
/*   of the sample, and the starting random-number seed: */

#define MAXN 100           /* Largest possible sample size value */
#define NBOOT 1000         /* Number of bootstrap replications */

int nn=30, nboot=NBOOT;    /* nn must be MAXN or less */
unsigned long wseed=23456;  /* Starting random-number seed */

double Mu=3.0;             /* The mean and standard deviation */
double Sdev=5.0;           /*   of the samples */

/* If we wanted to put in additional work, we could use these values as
     defaults and allow the option of changing them on the command line,
     with a syntax perhaps something like
                                       
        JknifeBstrap -n50 -b10000 -mu14.7 -w44444
                                       
     to set nn=50, nboot=10000, startseed=44444, etc.

  `mu' `Mu' and `MU' are different symbols to C, so that all can have
     different values in a program. Here MU is changed by the C
     preprocessor to the string `3.0', so that MU is not
     (and could not be) a variable at all. */


/* Prototypes of functions used in this program: */

/* Display an array in %7.3f format with 8 values per line: */
/* `const' means that the function promises not to change the values */
/*   in the array */
void show_array(int nn, const double *array);

/* Save the values for a particular confidence interval in a table
     to be displayed later */
void table_ci (const char *descript,
   double estvalue, double clow, double chigh);
void table_ci_tstat (const char *descript,
   double estvalue, double clow, double chigh);

/* Show the tabled CI values */
void show_summary(void);

/* The main (starting) function itself: */

int main(int argc, char *argv[])
 { /* The syntax `*xvar' below means that `xvar' is a POINTER to a double, */
   /*   and is not a double itself. We will later store the ADDRESS of an */
   /*   array that is large enough to store NN doubles in `xvar', so that */
   /*   xvar[0], xvar[1], ..., xvar[NN-1] will be valid doubles. */
   /* If `xvar' is a pointer to a double value, then we can fetch that */
   /*   value by saying either  *xvar  or  xvar[0] . */
   /* The language C usually does not distinguish between the name of */
   /*   an array and a pointer to the first element in that array, */
   /*   but there are rare instances in which it does. */ 
   int i,j;
   double xvar[MAXN];
   /* Variables for calculating standard deviations and their CIs */
   double s1,s2,ss, Obs_xbar,Obs_sdev,Obs_logsdev;
   double sdev,logsdev, tcrit, sbar,swid;
   double p1,pmed,p2, clow,cmed,chigh, expclow,expsbar,expchigh;
   /* Variables used in jacknife calculations: */
   double jdel_one[MAXN], jdel_logone[MAXN];
   double pseudoval[MAXN], pseudologval[MAXN];
   /* Variables used in bootstrap calculations: */
   int ns,m, kk,klow,khigh;
   double bootarray[NBOOT];
   double bmed, blow, bhigh;

   printf (
   "An example of using the `jackknife' and `bootstrap' estimation:\n"
   "  Finding bias-corrected estimates and confidence intervals for\n"
   "  the true standard deviation from a sample.\n");

   /* Enter `ProgName NNNN' (for example, Progname 1234) to initialize
      the random-number seet at NNNN instead of the default value */
   if (argc>=2) wseed=atol(argv[1]);
   wseed=setrand(wseed);   /* Returns `wseed' if wseed!=0 */
   printf ("The starting random-number seed is %lu\n", wseed);

   /* Given an example of `nn' random N(0,1)s (standard normals): */
   /* normrand() is a function in statlib.c */
   printf ("\nA sample of %d standard normal R.V.s is:\n",nn);
   for (i=0; i<nn; i++)
    { double xx=normrand();
      if (i>0 && (i%8)==0) printf("\n");
      printf ("  %7.4f", xx); }  printf("\n");

   /* Fill xvar[] with `nn' N(Mu,Sdev^2)s and then display them.*/
   /* If X is N(0,1), recall that Y = mu+C*X is N(mu,C^2). That is, */
   /*   Y has mean E(Y)=mu and standard deviation  C. */
   printf (
   "\nA sample of %d normal variables with mean=%g and st.dev.=%g is\n",
       nn,Mu,Sdev);
   /* Generate the random sample that will be analyzed */
   for (i=0; i<nn; i++)  xvar[i]=Mu+Sdev*normrand();

   /* Display the sample */
   show_array(nn,xvar);

   /* Find and display the sample mean and sample variance: */
   s1=s2=0.0;
   for (i=0; i<nn; i++) { s1+=xvar[i];  s2+=xvar[i]*xvar[i]; }
   Obs_xbar=s1/nn;
   ss=(1.0/(nn-1))*(s2 - s1*s1/nn);  Obs_sdev=sqrt(ss);
   printf ("OBSERVED sample mean: %.4f  Sample standard deviation: %.4f\n",
     Obs_xbar, Obs_sdev);

   /* Assuming that the sample is normal, find the `normal-theory'
        95% symmetric confidence interval for the true standard deviation: */
   printf ("\nCLASSICAL NORMAL-THEORY CONFIDENCE INTERVAL:\n");

   /* Find the 0.025 and 0.975 quantiles of Chi^2(nn-1)  */
   /* invchisq() is in statlib.c */
   p1=  invchisq(nn-1, 0.025);
   pmed=invchisq(nn-1, 0.500);
   p2=  invchisq(nn-1, 0.975);
   printf ("The 0.025, 0.500, and 0.975 quantiles of ChiSquare(%d) are\n"
     "  %g  %g  %g\n", nn-1, p1, pmed, p2);
   /* The endpoints of the classical 95% CI for the standard deviation:
      Note the reversal of endpoints due to the `upside down' derivation
        of exact confidence intervals.
      Note  1/p2 < 1/pmed < 1/p1  since  p1 < pmed < p2  */
   clow= Obs_sdev*sqrt((nn-1)/p2);
   cmed= Obs_sdev*sqrt((nn-1)/pmed);
   chigh=Obs_sdev*sqrt((nn-1)/p1);
   printf ("The classical symmetric 95%% CI for sdev is (%7.4f, %7.4f)\n",
     clow,chigh);
   printf ("\n(This is an exact CI, but assumes normal variables. It is based\n"
     "  on quantiles of the chi-square distribution with\n"
     "  df=%d degrees of freedom.)\n", nn-1);
   /* Save these values in the normal-theory SAMPEST structure
        for a later summary */
   table_ci("NormalTheory", cmed, clow,chigh);

   /* Now do a delete-1 jackknife for both sdev and log(sdev) */
   printf ("\nJACKKNIFE CONFIDENCE INTERVALS:\n");
   /* Within the loop below, we will first apply the functionals */
   /*   (sdev and log(sdev)) to delete-1 samples and store the results */
   /*   in `jdel_one' and `jdel_logone': */
   /* At the same time, we calculate the pseudovalues */
   /*    ps[i] = nn*(full value) - (nn-1)*(delete-1 value) */
   /* Full value for log(sdev): */
   Obs_logsdev=log(Obs_sdev);
   /* Compute the delete-1 values and pseudovalues */
   /*   for both sdev and log(sdev): */
   for (i=0; i<nn; i++)
    { double s1=0.0, s2=0.0, ss, sdev;
      /* Compute the sample standard deviation (sdev) of the sample with */
      /*   the i-th value dropped: */
      for (j=0; j<nn; j++) if (j!=i)
       { s1+=xvar[j];  s2+=xvar[j]*xvar[j]; }
      /* The delete-one sample size is nn-1, not nn */
      ss=(1.0/(nn-2))*(s2-s1*s1/(nn-1));
      /* Store both `sdev' and `log(sdev)' in jdel_one[] and jdel_logone[]: */
      jdel_one[i] = sdev = sqrt(ss);
      jdel_logone[i] = log(sdev);
      /* Store the pseudovalues as well */
      pseudoval[i] = nn*Obs_sdev - (nn-1)*jdel_one[i]; 
      pseudologval[i] = nn*Obs_logsdev - (nn-1)*jdel_logone[i]; }
   /* Display the delete-1 sample values and psuedovalues in both cases: */
   printf ("For the jackknife of sdev itself:\n");
   printf ("Est Sdev=%.4f: The sdev delete-1 jackknife values are:\n",
       Obs_sdev);
   show_array(nn,jdel_one);
   printf ("The sdev jackknife pseudo-values are:\n");
   show_array(nn,pseudoval);
   /* Find the mean and a classical 95% CI for the pseudovalues: */
   s1=s2=0.0;
   for (i=0; i<nn; i++)
    { s1+=pseudoval[i];  s2+=pseudoval[i]*pseudoval[i]; }
   sbar=s1/nn;  ss=(1.0/(nn-1))*(s2-s1*s1/nn);  sdev=sqrt(ss);
   printf ("The mean and sample sdev of the %d pseudovalues are:  %g  %g\n",
     nn, sbar,sdev);
   swid=1.960*sqrt(ss/nn);
   clow=sbar-swid;  chigh=sbar+swid;
   printf ("The pseudovalue 95%% symmetric CI for sdev is "
      "(%7.4f, %7.4f)\n", clow,chigh);
   /* Remember these values */
   table_ci("Sdev Jackknife", sbar, clow,chigh);

   /* What if we use a Student-t confidence interval instead? */
   tcrit=crit_tdist(nn-1,0.05);
   swid=tcrit*sqrt(ss/nn);
   clow=sbar-swid;  chigh=sbar+swid;
   printf ("\nUsing Student-t (df=%d) 95%% CI instead: "
      "(%7.4f, %7.4f)\n", nn-1, clow,chigh);
   /* Remember these values also */
   table_ci_tstat("Sdev Jackknife", sbar, clow,chigh);


   printf ("\nFor the jackknife of log(sdev):\n");
   printf ("Log(sdev)=%.4f: The %d logsdev delete-1 jackknife values are:\n", 
     logsdev, nn);
   show_array(nn,jdel_logone);
   printf ("The log(sdev) jackknife pseudo-values are:\n");
   show_array(nn,pseudologval);
   /* Find the mean and a classical 95% CI for the logged pseudovalues: */
   s1=s2=0.0;
   for (i=0; i<nn; i++)
    { s1+=pseudologval[i];  s2+=pseudologval[i]*pseudologval[i]; }
   sbar=s1/nn;  ss=(1.0/(nn-1))*(s2-s1*s1/nn);  sdev=sqrt(ss);
   printf ("The mean and sample sdev of the %d log pseudovalues are:  %g  %g\n",
     nn,sbar,sdev);
   swid=1.960*sqrt(ss/nn);
   clow=sbar-swid;  chigh=sbar+swid;
   printf ("The pseudovalue 95%% symmetric CI for logsdev is "
      "(%7.4f, %7.4f)\n", clow,chigh);
   expclow=exp(clow);  expsbar=exp(sbar);  expchigh=exp(chigh);
   printf ("The exponentiated mean and exponentiated 95%% symmetric CI are\n"
      "  %7.4f  (%7.4f, %7.4f)\n", expsbar, expclow,expchigh);
   /* Remember these values */
   table_ci("LogSdev Jackknife", expsbar, expclow,expchigh);

   /* What if we use a Student-t confidence interval instead? */
   swid=tcrit*sqrt(ss/nn);
   clow=sbar-swid;  chigh=sbar+swid;
   expclow=exp(clow);  expchigh=exp(chigh);
   printf ("\nUsing Student-t (df=%d) exponentiated 95%% CI instead: "
      "(%7.4f, %7.4f)\n", nn-1, expclow, expchigh);
   /* Remember these values also */
   table_ci_tstat("LogSdev Jackknife", expsbar, expclow,expchigh);

   /* Now do a bootstrap: */
   printf ("\nBOOTSTRAP CONFIDENCE INTERVAL:\n");
   /* Recall that commnum() (short for `comma-num') formats (for example) */
   /*   1234567 as 1,234,567. It requires %s (not %d) in a printf statement */
   printf (
   "Doing %s bootstrap replications for the sample standard deviation:\n",
     commnum(nboot));
   /* Recall that a bootstrap (re)sample for a sample of size n consists */
   /*   of n values sampled WITH REPLACEMENT from the original sample */
   /* Loop over the number of bootstrap replications: */
   for (ns=0; ns<nboot; ns++)
    { /* Build a bootstrap replicated sample in `tempval[]' */
      double tempval[MAXN];
      for (i=0; i<nn; i++)
        tempval[i] = xvar[nrand(nn)];
      /* Compute the sample variance */
      s1=s2=0.0;
      for (i=0; i<nn; i++)
       { s1+=tempval[i];  s2+=tempval[i]*tempval[i]; }
      ss=(1.0/(nn-1))*(s2-s1*s1/nn);
      /* Store the bootstrap `sdev' in bootarray[] at offset ns */
      bootarray[ns] = sqrt(ss); }
   /* Sort the array to find the median and quantiles using the */
   /*   standard C library `quick sort' function, `qsort()'. */
   /*   The `qsort' comparison function `dubsort()' is in statlib.c */
   /* If we used a `bubble sort' instead of a `quick sort', there */
   /*   could be a noticeable delay at this point. */
   printf ("Sorting the %s bootstrap replications:\n", commnum(nboot));
   qsort(bootarray,nboot, sizeof(double),dubsort);
   /* Show the first nn and last nn values of the bootstrap array. */
   /* These sets of values should be noticeably different: */
   printf ("The first %d sorted bootstrap sdev values are (#%d-#%d):\n",nn,1,nn);
   show_array(nn,bootarray);
   printf ("The last  %d sorted bootstrap sdev values are (#%d-#%d):\n",
     nn,nboot-nn+1,nboot);
   /* If xvar[] is an array, then `xvar' is a pointer to the beginning */
   /*   of the array (that is, to xvar[0]) and `xvar+m' is a pointer to */
   /*   the m-th element of the array, that is, to xvar[m]. */
   show_array(nn,bootarray+(nboot-nn));

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

   /* Assuming that `nboot' is even, calculate the median: */
   bmed = (bootarray[nboot/2] + bootarray[(nboot-2)/2])/2;
   /* Now the bootstrap `percentile' CI endpoints and display: */
   /* 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 */
   blow  = bootarray[klow];
   bhigh = bootarray[khigh];
   printf ("Bootstrap median: %.4f   95%% Percentile CI: (%7.4f, %7.4f)\n",
     bmed, blow,bhigh);
   /* Remember these values */
   table_ci("Bootstrap", bmed, blow,bhigh);

   /* Display the information that we stored in the 4 SAMPESTS structures */
   /*   in a single table: */
   printf ("\nSUMMARY  (Estimator of sdev, 95%% symmetric CI for sdev):\n");
   show_summary();

   /* Repeat basic parameters for clarity */
   printf ("\nSample size: %d   Number of bootstrap replications: %s\n",
     nn, commnum(nboot));
   printf ("Starting random-number seed: %lu\n", wseed);

   printf ("\nThe `Sdev Jackknife' CI seems to track the exact CI best,\n"
    "  even for different starting seeds and different numbers of\n"
    "  bootstrap replication.\n");

   show_randnums_used();

   return 0; }  /* Return with no errors */




/* Display an array in %7.3f format with 8 values per line: */

void show_array(int nn, const double *array)
 { int i;
   for (i=0; i<nn; i++)
    { if (i>0 && (i%8)==0) printf("\n");
      printf ("  %7.3f", array[i]); }
   printf("\n"); }


#define MAXT 20

int ntab;

const char *table_descript[MAXT];
double table_estvalue[MAXT];
double table_ci_low[MAXT];
double table_ci_high[MAXT];
int table_tstat[MAXT];

void table_ci (const char *descript,
   double estvalue, double clow, double chigh)
 { if (ntab<MAXT)
    { /* Copy values to a new row in the Summary table */
      table_descript[ntab]=descript;
      table_estvalue[ntab]=estvalue;
      table_ci_low[ntab]=clow;  table_ci_high[ntab]=chigh;
      table_tstat[ntab]=FALSE;  ntab++; }}


void table_ci_tstat (const char *descript,
   double estvalue, double clow, double chigh)
 { if (ntab<MAXT)
    { /* Copy values to a new row in the Summary table */
      table_descript[ntab]=descript;
      table_estvalue[ntab]=estvalue;
      table_ci_low[ntab]=clow;  table_ci_high[ntab]=chigh;
      table_tstat[ntab]=TRUE;  ntab++; }}


/* %-17s means use 17 spaces for the string and left justify. */
/* In this way, the descriptions and numbers will line up in columns. */

void show_summary(void)
 { int i;
   /* First, normal-theory (1.960) CIs: */
   for (i=0; i<ntab; i++) if (table_tstat[i]==FALSE)
     printf ("  %-17s  %7.4f  (%7.4f, %7.4f)\n",
       table_descript[i], table_estvalue[i],
       table_ci_low[i], table_ci_high[i]);
   printf ("Student-t CIs:\n");
   for (i=0; i<ntab; i++) if (table_tstat[i]==TRUE)
     printf ("  %-17s  %7.4f  (%7.4f, %7.4f)\n",
       table_descript[i], table_estvalue[i],
       table_ci_low[i], table_ci_high[i]); }