//======================================================== file = lec5a.c =====
//=  Program to demonstrate inversion method of generating an RV              =
//=============================================================================
//=  Build: bcc32 lec5a.c, cl lec5a.c, gcc lec5a.c                            =
//=---------------------------------------------------------------------------=
//=  Execute: lec5a                                                           =
//=---------------------------------------------------------------------------=
//=  History: KJC (05/29/11) - Refresh from old version                       =
//=============================================================================
//----- Include files ---------------------------------------------------------
#include <stdio.h>              // Needed for printf()
#include <stdlib.h>             // Needed for rand() and RAND_MAX
#include <math.h>               // Needed for log()

//----- Include files ---------------------------------------------------------
#define NUM_SAMPLES    100000    // The number of samples to generate

//----- Function prototypes ---------------------------------------------------
double exponential(double x);    // Returns an exponential random variable
double compute_mean(double *x);
double compute_variance(double *x, double mean);

//===== Main program ==========================================================
void main()
{
  double lambda;                // Exponential rate parameter
  double x[NUM_SAMPLES];        // Number of samples
  double mean, variance;        // Mean and variance of x
  int    i;                     // Loop counter

  // Set lambda value
  lambda = 0.5;

  // Generate samples
  for (i=0; i<NUM_SAMPLES; i++)
    x[i] = exponential(lambda);

  // Compute mean and variance for series x
  mean = compute_mean(x);
  variance = compute_variance(x, mean);

  // Output mean and variance of series x
  printf("lambda      = %f \n", lambda);
  printf("NUM_SAMPLES = %d \n", NUM_SAMPLES);
  printf("mean        = %f \n", mean);
  printf("variance    = %f \n", variance);
}

//=============================================================================
//=  Function to generate exponentially distributed RVs                       =
//=    - Input:  lamnda (exponential rate parameter)                          =
//=    - Output: Returns with exponential RV                                  =
//=============================================================================
double exponential(double lambda)
{
  double z;     // Uniform random number from 0 to 1

  // Pull a uniform RV (0 < z < 1)
  do
  {
    z = ((double) rand() / RAND_MAX);
  }
  while ((z == 0) || (z == 1));

  return(-(1.0 / lambda) * log(z));
}

//===========================================================================
//=  Function to compute mean for a series x                                =
//===========================================================================
double compute_mean(double *x)
{
  double   mean;        // Computed mean value to be returned
  int      i;           // Loop counter

  // Loop to compute mean
  mean = 0.0;
  for (i=0; i<NUM_SAMPLES; i++)
    mean = mean + (x[i] / NUM_SAMPLES);

  return(mean);
}

//===========================================================================
//=  Function to compute variance for a series x                            =
//===========================================================================
double compute_variance(double *x, double mean)
{
  double   variance;    // Computed variance value to be returned
  int      i;           // Loop counter

  // Loop to compute variance
  variance = 0.0;
  for (i=0; i<NUM_SAMPLES; i++)
    variance = variance + (pow((x[i] - mean), 2.0) / NUM_SAMPLES);

  return(variance);
}