/******* na-na na-na na-na na-na  STAT-MAN!!! *************************************/
/* stats.c, by Randy Wilson, September 26, 1995.
   Provides routines to compute t-test values, statistical significants, standard
   deviations, etc. */

#include <stdlib.h>
#include <stdio.h>
#include <math.h>

/** t[d][p] = probability that T < 't' (computed from a function like TPair(x[],y[],n)
       with confidence prob[p], given that you have 'd' degrees of freedom
       (or n=d+1 samples.). */
double prob[10] = {0.60, 0.70, 0.80, 0.90, 0.95, 0.975, 0.99, 0.995, 0.999, 0.9995};
int    freedom[38] = { 1,  2,  3,  4,  5,  6,  7,  8,  9, 10,
                      11, 12, 13, 14, 15, 16, 17, 18, 19, 20,
                      21, 22, 23, 24, 25, 26, 27, 28, 29, 30,
                      40, 50, 60, 80, 100, 200, 500, 0};
double t[38][10] = {
       {0.325, 0.727, 1.376, 3.078, 6.314, 12.710, 31.820, 63.660, 318.300, 636.600}, /* 1 */
       {0.289, 0.617, 1.061, 1.886, 2.920, 4.303, 6.965, 9.925, 22.330, 31.600},     /* 2 */
       {0.277, 0.584, 0.978, 1.638, 2.353, 3.182, 4.541, 5.841, 10.220, 12.940},    /* 3 */
       {0.271, 0.569, 0.941, 1.533, 2.132, 2.776, 3.747, 4.604, 7.173, 8.610},     /* 4 */
       {0.267, 0.559, 0.920, 1.476, 2.015, 2.571, 3.365, 4.032, 5.893, 6.859},    /* 5 */
       {0.265, 0.553, 0.906, 1.440, 1.943, 2.447, 3.143, 3.707, 5.208, 5.959},     /* 6 */
       {0.263, 0.549, 0.896, 1.415, 1.895, 2.365, 2.998, 3.499, 4.785, 5.405},    /* 7 */
       {0.262, 0.546, 0.889, 1.397, 1.860, 2.306, 2.896, 3.355, 4.501, 5.041},   /* 8 */
       {0.261, 0.543, 0.883, 1.383, 1.833, 2.262, 2.821, 3.250, 4.297, 4.781},  /* 9 */
       {0.260, 0.542, 0.879, 1.372, 1.812, 2.228, 2.764, 3.169, 4.144, 4.587}, /* 10 */
       {0.260, 0.540, 0.876, 1.363, 1.796, 2.201, 2.718, 3.106, 4.025, 4.437},     /* 11 */
       {0.259, 0.539, 0.873, 1.356, 1.782, 2.179, 2.681, 3.055, 3.930, 4.318},    /* 12 */
       {0.259, 0.538, 0.870, 1.350, 1.771, 2.160, 2.650, 3.012, 3.852, 4.221},   /* 13 */
       {0.258, 0.537, 0.868, 1.345, 1.761, 2.145, 2.624, 2.977, 3.787, 4.140},  /* 14 */
       {0.258, 0.536, 0.866, 1.341, 1.753, 2.131, 2.602, 2.947, 3.733, 4.073}, /* 15 */
       {0.258, 0.535, 0.865, 1.337, 1.746, 2.120, 2.583, 2.921, 3.686, 4.015},     /* 16 */
       {0.257, 0.534, 0.863, 1.333, 1.740, 2.110, 2.567, 2.898, 3.646, 3.965},    /* 17 */
       {0.257, 0.534, 0.862, 1.330, 1.734, 2.101, 2.552, 2.878, 3.611, 3.922},   /* 18 */
       {0.257, 0.533, 0.861, 1.328, 1.729, 2.093, 2.539, 2.861, 3.579, 3.883},  /* 19 */
       {0.257, 0.533, 0.860, 1.325, 1.725, 2.086, 2.528, 2.845, 3.552, 3.850}, /* 20 */
       {0.257, 0.532, 0.859, 1.323, 1.721, 2.080, 2.518, 2.831, 3.527, 3.819},     /* 21 */
       {0.256, 0.532, 0.858, 1.321, 1.717, 2.074, 2.508, 2.819, 3.505, 3.792},    /* 22 */
       {0.256, 0.532, 0.858, 1.319, 1.714, 2.069, 2.500, 2.807, 3.485, 3.767},   /* 23 */
       {0.256, 0.531, 0.857, 1.318, 1.711, 2.064, 2.492, 2.797, 3.467, 3.745},  /* 24 */
       {0.256, 0.531, 0.856, 1.316, 1.708, 2.060, 2.485, 2.787, 3.450, 3.725}, /* 25 */
       {0.256, 0.531, 0.856, 1.315, 1.706, 2.056, 2.479, 2.779, 3.435, 3.707},     /* 26 */
       {0.256, 0.531, 0.855, 1.314, 1.703, 2.052, 2.473, 2.771, 3.421, 3.690},    /* 27 */
       {0.256, 0.530, 0.855, 1.313, 1.701, 2.048, 2.467, 2.763, 3.408, 3.674},   /* 28 */
       {0.256, 0.530, 0.854, 1.311, 1.699, 2.045, 2.462, 2.756, 3.396, 3.659},  /* 29 */
       {0.256, 0.530, 0.854, 1.310, 1.697, 2.042, 2.457, 2.750, 3.385, 3.646}, /* 30 */
       {0.255, 0.529, 0.851, 1.303, 1.684, 2.021, 2.423, 2.704, 3.307, 3.551}, /* 40 */
       {0.255, 0.528, 0.849, 1.298, 1.676, 2.009, 2.403, 2.678, 3.262, 3.495}, /* 50 */
       {0.254, 0.527, 0.848, 1.296, 1.671, 2.000, 2.390, 2.660, 3.232, 3.460}, /* 60 */
       {0.254, 0.527, 0.846, 1.292, 1.664, 1.990, 2.374, 2.639, 3.195, 3.415}, /* 80 */
       {0.254, 0.526, 0.845, 1.290, 1.660, 1.984, 2.365, 2.626, 3.174, 3.389}, /* 100 */
       {0.254, 0.525, 0.843, 1.286, 1.653, 1.972, 2.345, 2.601, 3.131, 3.339}, /* 200 */
       {0.253, 0.525, 0.842, 1.283, 1.648, 1.965, 2.334, 2.586, 3.106, 3.310}, /* 500 */
       {0.253, 0.524, 0.842, 1.282, 1.645, 1.960, 2.326, 2.576, 3.090, 3.291}}; /* Infinity */

/************************ TPair **************************************************/
double TPair(double *x, double *y, int n)
{ /* Given two arrays x[0..n-1] and y[0..n-1], returns the one-tailed paired t-test value,
     using a null hypothesis that the average of x is not greater than the average of y. */
  double average; /* average of differences between x[i] and y[i] */
  double sum; 
  double sd;      /* standard deviation of diff. betw. xi and yi */
  double diff;
  int i;
  
  if (n<=1)
    { printf("Come on, why do a t-test on one value??");
      return(0.0);
    }
  sum=0.0;
  for (i=0; i<n; i++)
    sum = sum + (x[i]-y[i]);
  average = sum/((double)n);
  sd = 0.0;
  for (i=0; i<n; i++)
    { diff = (x[i]-y[i]) - average;
      sd = sd + diff*diff;
    }
  sd = sd / ((double)(n-1));
  sd = sqrt(sd);
  
  if (average==0.0)
    return(0.0);
  if (sd==0.0)
    { printf("Standard deviation of diff=0.0, average diff=%.2lf => Don't trust t-test.\n",
              average);
      if (average > 0)
        return(sqrt((double)n));
      else return(sqrt(-(double)n));
    }
  return((average - 0.0) / (sd / sqrt((double)n)));
} /* TPair */

/*********************** Confidence *********************************************/
double Confidence(double *x, double *y, int n)
{ /* Using a one-tailed student t-test, returns the statistical confidence (0..1) that
       average(x) is really greater than average(y). If instead y has the greater
       average, then a NEGATIVE number is returned, indicating the confidence (0..-1)
       that average(x) is really less than average(y). */
  int v,i;
  double tp;
  double sign;
  
  if (n<2)
    { printf("You need at least two values to compute a confidence.\n");
      return(0.0);
    }
  tp = TPair(x,y,n);
  if (tp<0.0)
     sign = -1.0;
  else sign = 1.0;
  
  tp = tp * sign; /* Make it positive for now */
  /* Find the index v into 'freedom[]', the proper number of degrees of freedom. */
  for (v=0; v<37 && freedom[v]<n-1; v++)
    ; /* do nothing */

  /* Find the index 'i' into the probability array. */
  for (i=0; i<10 && t[v][i]<tp; i++)
    ; /* do nothing */
  
  if (i==0)
    return(0.0);
  else return(sign * prob[i-1]);
} /* Confidence */
    
/******************** Significant *********************************************/     
int Significant(double *x, double *y, int n, double p)
{ /* Using a confidence value 'p' and a one-tailed student t-test...
      Returns 1 if the average of x is significantly greater than the average of y;
      Returns -1 if the average of x is significantly lower than the average of y;
      Returns 0 if there is no significant difference. */
  double c;
  c = Confidence(x,y,n);
  if (c<0.0)
    { c = -c;
      if (c>=p)
        return(-1);
      else return(0);
    }
  else
    { if (c>=p)
        return(1);
      else return(0);
    }
  printf("How did this line get executed??");
  return(0);
} /* Significant */
   
/*********************** StandardDeviation *******************************************/
double StandardDeviation(double *x, int n)
{ /* Returns the Standard Deviation of x[0]..x[n-1]. Returns '100.0' if n<=1. */
  double sum, average, sd, diff;
  int i;
  sum = 0.0;
  if (n<=1)
     return(100.0);
  for (i=0; i<n; i++)
    sum = sum + x[i];
  average = sum / ((double) n);
  sd = 0.0;
  for (i=0; i<n; i++)
  { diff = x[i] - average;
    sd = sd + (diff * diff);
  }
  sd = sd/((double)(n-1));
  sd = sqrt(sd);
  return(sd);
} /* Standard Deviation */

/**************** Average ****************************/
double Average(double *x, int n)
{ /* Returns the average of n numbers x[0]..x[n-1] */
  double sum=0.0;
  int i;
  for (i=0; i<n; i++)
    sum = sum + x[i];
  if (n<1)
    return(0);
  else return(sum/((double)n));
} /* Average */