/* RAND-TEST.C - Program to test random number generators. *//* Copyright (c) 1995 by Radford M. Neal  * * Permission is granted for anyone to copy, use, or modify this program  * for purposes of research or education, provided this copyright notice  * is retained, and note is made of any changes that have been made.  * * This program is distributed without any warranty, express or implied. * As this program was written for research purposes only, it has not been * tested to the degree that would be advisable in any important application. * All use of this program is entirely at the user's own risk. *//* Usage:     rand-test seed generator { parameters } / sample-size [ low high bins ]Using the seed given, tests the random number generator identified bythe second argument, for the parameter values specified.  The possiblegenerators and required parameters are as follows:    uniform                  Uniform from [0,1)    uniopen                  Uniform from (0,1)    int n                    Uniform from the set { 0, 1, ..., (n-1) }    gaussian                 From Gaussian with mean zero and unit variance    exp                      From exponential with mean one    cauchy                   From Cauchy centred at zero with unit width    gamma alpha              From Gamma with shape parameter (and mean) alpha    beta a b                 From Beta with parameters a and bThe size of the sample to use is also specified.  The program reportsthe mean and variance of the sample.  A histogram is also printed if alow and high range and number of bins are for it are specified.These tests are not really adequate to detect subtle forms of bias dueto use of pseudo-random numbers, but are hopefully good enough to findmost programming errors. */#include <stdlib.h>#include <stdio.h>#include <string.h>#include <math.h>#include "rand.h"#define Max_bins 1000		/* Maximum number of histogram bins */static void usage (void);/* MAIN PROGRAM. */main( int argc,  char **argv){  int seed, sample_size, bins, np;  double low, high;  char *generator;  double p1, p2;  double mean, variance;  double tmean, tvariance;  int undef_mean, undef_variance;  int count[Max_bins];  int under, over;  char **ap;  double x;  int i, n;  bins = 0;  if (argc<5) usage();  if ((seed = atoi(argv[1]))==0 && strcmp(argv[1],"0")!=0) usage();  generator = argv[2];  if      (strcmp(generator,"uniform")==0)  np = 0;  else if (strcmp(generator,"uniopen")==0)  np = 0;  else if (strcmp(generator,"int")==0)      np = 1;  else if (strcmp(generator,"gaussian")==0) np = 0;  else if (strcmp(generator,"exp")==0)      np = 0;  else if (strcmp(generator,"cauchy")==0)   np = 0;  else if (strcmp(generator,"gamma")==0)    np = 1;  else if (strcmp(generator,"beta")==0)     np = 2;  else  { fprintf(stderr,"Unknown generator: %s\n",generator);    exit(1);  }  ap = argv+3;  if (np>0)   { if (*ap==0 || (p1 = atof(*ap++))<=0) usage();  }  if (np>1)   { if (*ap==0 || (p2 = atof(*ap++))<=0) usage();  }  if (*ap==0 || strcmp(*ap++,"/")!=0) usage();  if (*ap==0 || (sample_size = atoi(*ap++))<=0) usage();  if (*ap!=0)  { low = atof(*ap++);    if (*ap==0) usage();    high = atof(*ap++);    if (high<=low) usage();    if (*ap==0 || (bins = atoi(*ap++))<=0) usage();    if (bins>Max_bins)     { fprintf(stderr,"Too many histogram bins\n");      exit(1);    }  }  if (*ap!=0) usage();  printf("\nTest of %s(",generator);  if (np>0) printf("%.4f",p1);  if (np>1) printf(",%.4f",p2);  printf(") generator using sample of size %d with seed %d\n\n",   sample_size, seed);  undef_mean = undef_variance = 0;  if (strcmp(generator,"uniform")==0)    { tmean = 0.5;    tvariance = 1.0/12.0;  }  else if (strcmp(generator,"uniopen")==0)    { tmean = 0.5;    tvariance = 1.0/12.0;  }  else if (strcmp(generator,"int")==0)        { tmean = (p1-1)/2;    tvariance = p1*p1/3.0 - p1/2.0 + 1/6.0 - tmean*tmean;  }  else if (strcmp(generator,"gaussian")==0)   { tmean = 0;    tvariance = 1;  }  else if (strcmp(generator,"exp")==0)   { tmean = 1;    tvariance = 1;  }  else if (strcmp(generator,"cauchy")==0)  { undef_mean = 1;    undef_variance = 1;  }  else if (strcmp(generator,"gamma")==0)  { tmean = p1;    tvariance = p1;  }  else if (strcmp(generator,"beta")==0)       { tmean = p1 / (p1+p2);    tvariance = (p1*p2) / ((p1+p2)*(p1+p2)*(p1+p2+1));  }  else   { abort();  }  mean = 0;  variance = 0;  if (bins>0)  { for (i = 0; i<bins; i++) count[i] = 0;    under = over = 0;  }  rand_seed(seed);  for (n = 0; n<sample_size; n++)  {    if      (strcmp(generator,"uniform")==0)  x = rand_uniform();    else if (strcmp(generator,"uniopen")==0)  x = rand_uniopen();    else if (strcmp(generator,"int")==0)      x = rand_int((int)p1);    else if (strcmp(generator,"gaussian")==0) x = rand_gaussian();    else if (strcmp(generator,"exp")==0)      x = rand_exp();    else if (strcmp(generator,"cauchy")==0)   x = rand_cauchy();    else if (strcmp(generator,"gamma")==0)    x = rand_gamma(p1);    else if (strcmp(generator,"beta")==0)     x = rand_beta(p1,p2);    else abort();    mean += x;    variance += x*x;    if (bins>0)    { if (x<low)       { under += 1;      }      else      { i = (int) ((x-low)/((high-low)/bins));        if (i>=bins)         { over += 1;        }        else         { count[i] += 1;        }      }    }  }  mean /= sample_size;  variance /= sample_size;  variance -= mean*mean;  printf("Sample mean:     %.4f",mean);  if (undef_mean)  { printf(" (true value: undefined)\n");  }  else  { printf(" (true value: %.4f)\n",tmean);  }  printf("Sample variance: %.4f",variance);  if (undef_variance)  { printf(" (true value: undefined)\n");  }  else  { printf(" (true value: %.4f)\n",tvariance);  }  printf("\n");  if (bins!=0)  { printf("Histogram:\n");    printf("                    under : %8d  %.5f\n\n",       under, (double)under / sample_size);    for (i = 0; i<bins; i++)    { printf("  %10.4f - %10.4f : %8d  %.5f\n",         i*(high-low)/bins + low, (i+1)*(high-low)/bins + low,         count[i], (double)count[i] / sample_size);    }    printf("\n                     over : %8d  %.5f\n",       over, (double)over / sample_size);    printf("\n");  }    exit(0);}/* PRINT USAGE MESSAGE AND EXIT. */static void usage (void){  fprintf(stderr,   "Usage: rand-test seed generator { parameters } / sample-size [ low high bins ]\n");  exit(1);}