/*
  CSC354 Problem Set 1, Question 5.

  Name:    Andria Hunter
  Login:   at354hun
  Stud.ID: 555555555

  Program Description:
    This CSIM program determines the goodness of the uniform
    random number generator in CSIM by applying the chi-square
    frequency test and the runs up and runs down test, for the
    (same) first 100 numbers.
*/

#include "csim.h"
#include "math.h"

#define   NUMS    100          /* number of random numbers */
#define   INTVL    10          /* number of intervals */

#define   CHI_CRIT 16.0        /* Critical value for chi-squared */
#define   Z0_CRIT  1.96        /* Critical value for runs test */

#define   NO_EVENT  0          /* runs-up-runs-down events */
#define   RUNS_UP   1
#define   RUNS_DOWN 2

HIST      his;                 /* pointer for uniform table */

/* --- Domestic Routines --- */

static void chi_squared();
static void runs_up_runs_down();

/* sim process - Execution starts here.  Create an array "arr" */
/*    to store 100 uniformly distributed random numbers.  Then */
/*    call the chi_squared and runs_up_runs_down functions to  */
/*    determine the goodness of the distribution. */

sim()
{
   int    i=0;                 /* for loop counter */
   double arr[NUMS];           /* stores uniform random numbers */

   set_model_name("Test Random Numbers");
   create("sim");

   his = histogram("Uniform Distribution",10,0.0,1.0);   /* init. histogram */

   /* Generate 100 uniform random numbers and store in "arr". */

   for(i=0 ; i<NUMS ; i++) {
      arr[i] = uniform (0.0,1.0);    /* random num between 0.0 & 1.0 */
      record (arr[i], his);          /* record random value in a table */
   }

#ifdef DEBUG
   /* Print contents of the array. */

   printf ("These are the uniform random numbers:\n");
   for(i=0 ; i<NUMS ; i++) {
      printf ("%lf\n", arr[i]);
   }
#endif DEBUG

   /* Perform Chi-Squared Test */
   chi_squared (arr);

   /* Perform "runs up runs down" test. */
   runs_up_runs_down (arr);

   /* Show the tables and histograms. */
   printf ("\n\nShow CSIM Statistics\n");
   printf ("====================\n\n");

   report();
}

/* chi_squared - This procedure performs the chi-square test.   */
/*    It splits the interval 0.0-1.0 into 10 equal parts and    */
/*    determines the number of random numbers in each interval. */
/*    Then it calculates chi-squared based on these intervals.  */

static
void chi_squared (arr)

double arr[NUMS];

{
   int    i;

   int    O[INTVL];             /* Store count for each interval */

   int    Ei = NUMS/INTVL;      /* length of intervals */

   int    O_sum=0;              /* Used for Chi-Test Computations */
   int    O_min_E_sum=0;
   double O_min_E_sq_sum=0.0;

   for(i=0 ; i<INTVL ; i++) {
      O[i] = 0;
   }

   /* Store interval counts in array "O". */

   for(i=0 ; i<NUMS ; i++) {

      /* The first significant digit of the random number is */
      /* used as the interval array index.  Increment the */
      /* interval array contents at this index. */

      O[(int)floor(arr[i]*10.0)]++;
   }

   /* Generate Table to show Chi-squared Computations */

   printf ("Chi-Square Frequency Test\n");
   printf ("=========================\n");

   printf ("                           2     2\n");
   printf ("   Int   Oi   Ei   O-E (O-E) (O-E)/E\n");
   printf ("  ---- ---- ---- ----- ----- ------\n");

   for(i=0 ; i<INTVL ; i++) {
      printf ("  %4d %4d %4d %5d %5d %6.2f\n", i+1, O[i], Ei, O[i]-Ei,
               (O[i]-Ei)*(O[i]-Ei), (O[i]-Ei)*(O[i]-Ei)/(double)Ei);

      O_sum          += O[i];
      O_min_E_sum    += (O[i]-Ei);
      O_min_E_sq_sum += (O[i]-Ei)*(O[i]-Ei)/(double)Ei;
   }

   printf ("       ---- ---- -----       ------\n");
   printf ("       %4d %4d %5d       %6.2f\n", O_sum, NUMS,
                         O_min_E_sum, O_min_E_sq_sum);

   printf ("\nChi-Square-0 is %lf\n", O_min_E_sq_sum);

   printf ("\nNote that for alpha=0.05, if the Chi-Square-0 value is\n");
   printf ("smaller than the critical value of 16.0, then the hypothesis of\n");
   printf ("no difference between the sample distribution and the uniform\n");
   printf ("distribution is not rejected.\n\n");

   if (O_min_E_sq_sum > CHI_CRIT)
      printf ("null hypothesis of uniformity is rejected.\n\n");
   else
      printf ("null hypothesis of uniformity is not rejected.\n\n");
}

/* runs_up_runs_down - This procedure performs the "runs up runs down" */
/*    test.  It determines the number of runs, and uses this to deter- */
/*    mine independence of the random numbers in the sequence.         */

static
void runs_up_runs_down (arr)

double arr[NUMS];

{
   int    i;
   double prev;                /* Previous random value */
   int    state = NO_EVENT;    /* current runs-up-runs-down state */
   int    num_runs = 0;        /* number of runs */

   double mean_a, var_a;       /* mean and variance of number of runs */

   double Z0;                  /* Standardized normal test statistic used */
                               /* to determine independence of numbers */

   prev = arr[0];  /* store first random value */

   /* Compare each random value to the previous value to determine */
   /* if we are starting a new run in the opposite direction.  If  */
   /* so, increase the number of runs.                             */

   for(i=1 ; i<NUMS ; i++) {
      if (arr[i]<prev && (state==RUNS_UP || state==NO_EVENT)) {
         num_runs++;
         state = RUNS_DOWN;
      }
      else if (arr[i]>prev && (state==RUNS_DOWN || state==NO_EVENT)) {
         num_runs++;
         state = RUNS_UP;
      }
      prev = arr[i];
   }

   /* Print Results */

   printf ("\n\nRuns Up Runs Down Test\n");
   printf ("======================\n");

   printf ("\nThere are %d runs in this sequence of random numbers.\n", num_runs);
   printf ("  (a=%d and N=%d)\n\n", num_runs, NUMS);

   mean_a = (2*NUMS-1)/3.0;
   var_a  = (16*NUMS-29)/90.0;
   Z0     = (num_runs - mean_a) / sqrt(var_a);

   printf ("       Mean for a = %lf\n", mean_a);
   printf ("   Variance for a = %lf\n", var_a);
   printf ("Test statistic Z0 = %lf\n", Z0);

   printf ("\nNote that for alpha=0.05, the hypothesis of independence\n");
   printf ("is rejected only if Z0<-1.96 or Z0>1.96.\n\n");

   if (Z0 > Z0_CRIT || Z0 < -1.0*Z0_CRIT)
      printf ("null hypothesis of independence is rejected.\n");
   else
      printf ("null hypothesis of independence is not rejected.\n");
}