testbino.cpp

很好用的随机数产生器

/*************************** TESTBINO.CPP ******************* AgF 2001-11-15 *
*                                                                            *
*  Test binomial distribution                                                *
*  Necessary files are in stocc.zip                                          *
*                                                                            *
*****************************************************************************/

#include <time.h>                      // define time()
#include "randomc.h"                   // define classes for random number generators
#include "mersenne.cpp"                // code for random number generator
#define RANDOM_GENERATOR TRandomMersenne // define which random number generator to use
#include "stocc.h"                     // define random library classes
#include "stoc1.cpp"                   // random library source code
#include "stoc2.cpp"                   // random library source code
#include "userintf.cpp"                // define system specific user interface


// main  
void main() {

  // parameters. You may change these
  long int n = 10;                     // number of balls
  double p = 0.4;                      // probability
  long int nn = 1000000;               // number of samples
  
  // other variables
  double sum;                          // sum
  double ssum;                         // squaresum
  long int min, max;                   // minimum, maximum
  double mean;                         // mean
  double var;                          // variance
  long int i;                          // loop counter
  long int x;                          // random variate
  const int DSIZE = 18;                // size of array
  long int dis[DSIZE] = {0};           // number of samples in each range
  int c;                               // index into dis list
  double f;                            // calculated function value
  double xme;                          // expected mean
  double xva;                          // expected variance
  double xsd;                          // expected standard deviation
  long int delta;                      // step in list
  long int xa;                         // list minimum
  long int x1, x2;                     // category range
  
  // make random library
  long int seed = time(0);             // random seed
  StochasticLib sto(seed);             // make instance of random library

  // calculate mean and variance
  xme = n*p;                           // calculated mean
  xva = n*p*(1-p);                     // calculated variance

  // calculate appropriate list divisions
  xsd = sqrt(xva);                     // calculated std.dev.
  xa = int(xme - 6.*xsd + 0.5);        // calculated minimum
  if (xa < 0) xa=0;
  delta = int(12.*xsd/(DSIZE-1));      // list step
  if (delta < 1) delta = 1; 

  // initialize
  sum = ssum = 0; min = 2000000000; max = -1;

  // start message
  printf("taking %li samples from binomial distribution...", nn);
  
  // sample nn times
  for (i=0; i < nn; i++) {

    // generate random number with desired distribution
    x = sto.Binomial(n,p);

    // update sums
    sum += x;
    ssum += (double)x*x;
    if (x < min) min = x;
    if (x > max) max = x;
    c = (x-xa)/delta;
    if (c < 0) c = 0;
    if (c >= DSIZE) c = DSIZE-1;
    dis[c]++;}
    
  // calculate sampled mean and variance
  mean = sum / nn; 
  var = (ssum - sum*sum/nn) / nn;

  // print sampled and theoretical mean and variance
  printf("\n\nparameters: n=%li, p=%.3G", n, p);
  printf("\n                mean        variance");
  printf("\nsampled:   %12.6f %12.6f", mean, var);
  printf("\nexpected:  %12.6f %12.6f", xme,  xva);

  // print found and expected frequencies
  printf("\n\n      x              found     expected");
  for (c = 0; c < DSIZE; c++) {
    x1 = c*delta + xa;
    if (x1 > n) break;
    x2 = x1+delta-1;
    if (c==0 && min<x1) x1 = min;
    if (c==DSIZE-1 && max>x2) x2 = max;
    if (x2 > n) x2 = n;

    // calculate expected frequency
    for (x=x1, f=0.; x <= x2; x++) {
      f += exp(LnFac(n) - LnFac(x) - LnFac(n-x) + x*log(p) + (n-x)*log(1-p));}

    // output found and expected frequency
    if (dis[c] || f*nn > 1E-4) {
      if (x1==x2) {
        printf("\n%7li       %12.6f %12.6f", x1, (double)dis[c]/nn, f);}
      else {
        printf("\n%6li-%-6li %12.6f %12.6f", x1, x2, (double)dis[c]/nn, f);}}}
    
  EndOfProgram();                      // system-specific exit code
  }