/*-------------------------------------------------------------

Copyright (C) 2000 Peter Clote, Sebastian Will 
All Rights Reserved.

written by Peter Clote
modified by Sebastian Will

Permission to use, copy, modify, and distribute this
software and its documentation for NON-COMMERCIAL purposes
and without fee is hereby granted provided that this
copyright notice appears in all copies.


THE AUTHOR AND PUBLISHER MAKE NO REPRESENTATIONS OR
WARRANTIES ABOUT THE SUITABILITY OF THE SOFTWARE, EITHER
EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE
IMPLIED WARRANTIES OF MERCHANTABILITY, FITNESS FOR A
PARTICULAR PURPOSE, OR NON-INFRINGEMENT. THE AUTHORS
AND PUBLISHER SHALL NOT BE LIABLE FOR ANY DAMAGES SUFFERED
BY LICENSEE AS A RESULT OF USING, MODIFYING OR DISTRIBUTING
THIS SOFTWARE OR ITS DERIVATIVES.

-------------------------------------------------------------*/


/*------------------------------

normalDist.c

This illustrates application of the central
limit theorem and the polar method (Box-Muller)
for trannsforming the uniform distribution on
(0,1) to the normal distribtion N(mean,variance)
having given mean and variance.

--------------------------------*/


#include <stdio.h>
#include <stdlib.h>  // for RAND_MAX
#include <math.h>    // for sqrt

#define M 100000.0
    // M sample points are generated
#define N 30.0
    // N is number of independent
    // random variables (used in normal())
    // SN = X1 + ... + XN

#define PI 3.14159

double normal (double, double);
double boxMuller (double, double);
double boxMuller2 (double, double);
void test_distrib(char *descr, double g(double,double),
		  double mean, double variance);

main(int argc, char *argv[]){
    double mean;
    double variance;

    if (argc != 3) {
       fprintf(stderr,"%s mean variance\n",argv[0]);
       exit(1);
    }

    mean = atof(argv[1]);
    variance = atof(argv[2]);

    test_distrib("Central limit - method",normal,mean,variance);
    
    test_distrib("BoxMuller - method 1",boxMuller,mean,variance);

    test_distrib("BoxMuller - method 2",boxMuller2,mean,variance);
}

/* test_distrib(char *descr, double g(double,double),
		  double mean, double variance)
   tests the distribution generated by then function g by determining
   the mean and variance of the distribution
   
   prints description string descr and results
*/

void test_distrib(char *descr, double g(double,double),
		  double mean, double variance)
{
    int i;

    double sum,sumSq, sample;
    double sampleMean, sampleVariance;
   
    sum = sumSq = 0;
    printf("\n\n%s\n",descr);
    for (i=0;i<M;i++) {
       sample = g(mean,variance);
       sum += sample;
       sumSq += sample*sample;
    }
    sampleMean = sum/M;
    sampleVariance = sumSq/M - sampleMean*sampleMean;

    printf("Mean\tVariance");
    printf("\n\n%lf\t%lf\n\n\n",sampleMean,sampleVariance);
}

/* normal ( double mean, double variance )
   
   returns random number from approximated gaussian distribution with
   given mean and variance approximates distribution by adding N
   random numbers from uniform distribution. (central-limit method)
*/
double normal ( double mean, double variance )
{
    int i;
    double x=0;

    for (i=0;i<N;i++) x += (double) rand()/ (RAND_MAX+1.0);
    return ( mean + sqrt(variance)*(x - N/2.0)/sqrt(N/12.0) );
}

/* boxMuller ( double mean, double variance ) 

   returns random number from gaussian distribution with given mean and
   variance using the polar method due to Box-Muller 
*/
double boxMuller ( double mean, double variance )
{
    double a,b,u,v,x,y,s;

    do {
       a = (double) rand()/(RAND_MAX+1.0);
       b = (double) rand()/(RAND_MAX+1.0);
                   // a,b uniform in (0,1)
       u = 2*a-1;
       v = 2*b-1;  // u,v uniform in (-1,1)
       s = u*u + v*v;
    } while (s >= 1);
  
                   // (u,v) uniform in unit circle 
    x = u*sqrt(-2*log(s)/s); 
    y = v*sqrt(-2*log(s)/s); 
    return ( mean + sqrt(variance)*x );
}

/* boxMuller ( double mean, double variance ) 

   returns random number from gaussian distribution with given mean and
   variance using the polar method due to Box-Muller 

   version of Box-Muller corresponding to the one presented in the book
*/
double boxMuller2 ( double mean, double variance )
{
  double u,v,x,y;
  x = (double) rand() / (RAND_MAX+1.0);
  y = (double) rand() / (RAND_MAX+1.0);
  u = mean + sqrt(variance) * (sqrt(-2*log(x))*cos(2*PI*y));
  return u;
}