test random number generators

test_random.cpp

// The practice of programming
// chapter 3 : problem 3-1
//
// for testing mersenne twister and linear congruential random number generators
// prints normalized frequencies that i'th item is selected.
// statistical errors are printed as well.
#include <iostream>
#include <boost/random.hpp>

const size_t MAX = 1024;         // number of elements to be selected
const size_t AVG_COUNT = 128;
const size_t NUM_TRY = MAX * AVG_COUNT;  // number of trials
const size_t NUM_SAMPLES = 64;       // number of samples to calculate average and its error

#define RANDTYPE boost::mt19937
//#define RANDTYPE boost::minstd_rand

void print_counts();
void test_random( std::vector<int> &vec_nCount, size_t seed);

int main() {

  print_counts();
  
  return 0;
}

//
//
void print_counts() {
  std::vector<int> vec_nCountTotal( MAX, 0);       // total number of selected times
  std::vector<int> vec_nCountSquareTotal( MAX, 0); // square of selected times (to calculate variance)

  //
  // count
  //
  for( size_t i=0; i<NUM_SAMPLES; i++) {
    std::cerr << "sample : " << i << " / " << NUM_SAMPLES << std::endl;

    std::vector<int> vec_nCount( MAX, 0);
    test_random( vec_nCount, i);

    for( size_t j=0; j<MAX; j++) {
      vec_nCountTotal[j] += vec_nCount[j];
      vec_nCountSquareTotal[j] += (vec_nCount[j] * vec_nCount[j]);
    }
  }

  //
  // calculate average and statistical erros
  //
  for( size_t j=0; j<MAX; j++) {
    double dAverage = (double)vec_nCountTotal[j] / (double)NUM_SAMPLES;
    double dVariance = ( vec_nCountSquareTotal[j] - NUM_SAMPLES*dAverage*dAverage ) / (NUM_SAMPLES-1);
    double dError = std::sqrt( dVariance / NUM_SAMPLES);
    std::cout << j << ' ' << dAverage/AVG_COUNT << ' ' << dError/AVG_COUNT << std::endl;
  }
}

// -------------------------------------------------
// select index [0:MAX-1] randomly
// count the times that i'th index is selected
// -------------------------------------------------
void test_random( std::vector<int> & vec_nCount, size_t seed) {
  RANDTYPE rand_gen(seed);
  boost::uniform_01<> dist01;
  boost::variate_generator< RANDTYPE, boost::uniform_01<> > zGen(rand_gen,dist01);
  
  for( size_t i=0; i<NUM_TRY; i++) {
    long iSelected = 0;
    for( size_t j=0; j<vec_nCount.size(); j++) {
      // double dProb = 1.0 / (j+1);
      // if( zGen() < dProb ) { iSelected = j; }
      if( (rand_gen() % (j+1)) == 0 ) { iSelected = j; }
    }
    vec_nCount[iSelected]++;
  }

  /*
  for( size_t i=0; i<MAX; i++) {
    std::cout << vec_nCount[i] << std::endl;
  }
  */

}