test_rng_wrapper.cpp

一个很全的matlab工具包

// $Id: test_rng_wrapper.cpp,v 1.14 2004/11/10 08:58:09 kgadeyne Exp $

// Copyright (C) 2002 Klaas Gadeyne <klaas dot gadeyne at mech dot kuleuven dot ac dot be>
//  
// This program is free software; you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation; either version 2 of the License, or
// (at your option) any later version.
//  
// This program is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
// GNU General Public License for more details.
//  
// You should have received a copy of the GNU General Public License
// along with this program; if not, write to the Free Software
// Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
//  

// KG test program of RNG wrapper (currently only libscythe is
// wrapped) 
// Very basic is the least you could say for now :-)

#include <wrappers/rng/rng.h>
#include <iostream>
#include <fstream>
#include <cmath>

using namespace std;
using namespace BFL;

#define NUM_EXPERIMENTS 20000

int main()
{
  // All samples, their sum, and the sum of quadratic differences
  float sample[NUM_EXPERIMENTS], sum=0, diffsum=0;
  // Population mean and std. dev
  double mean=0, standard_deviation=0;
  ofstream results("results.txt");
  
  cout << "Sampling from a Uniform RNG (marsaglia)" << endl;
  for (int num_exp = 11; num_exp < NUM_EXPERIMENTS ; num_exp *= 1.1)
    {
      sum = 0; diffsum = 0.0; mean = 0.0; standard_deviation = 0.0;
      for (int i = 0; i < num_exp; i++)
	{
	  sample[i] = runif();
	  sum += sample[i];
	}
      // The mean is
      mean=sum/num_exp;
  
      // Calculating the stand. dev.
      for (int j = 0; j < num_exp; j++)
	{
	  diffsum += std::pow(sample[j] - mean,2);
	}
  standard_deviation = diffsum / (num_exp-1);

//   cout << "number of experiments = " << num_exp << endl;
//   cout << "mean  = " << mean << endl;
//   cout << "std. deviation = " << standard_deviation << endl;
  results << num_exp << " " << mean << " " << standard_deviation << endl;
    }


  cout << "========================" << endl;
  cout << "Sampling from a Gaussian" << endl;


  const float MU=0.0;
  const float SIGMA=1.0;
  sum = 0; diffsum = 0.0; mean = 0.0; standard_deviation = 0.0;
  // Sampling and calculating the cum. sum
  for (int i = 0; i < NUM_EXPERIMENTS; i++)
    {
      sample[i] = rnorm(MU,SIGMA);
      sum += sample[i];
    }
  
  // The mean is
  mean=sum/NUM_EXPERIMENTS;
  
  // Calculating the stand. dev.
  for (int j = 0; j < NUM_EXPERIMENTS; j++)
    {
      diffsum += std::pow(sample[j] - mean,2);
    }

  standard_deviation = diffsum / (NUM_EXPERIMENTS-1);

  cout << "number of experiments = " << NUM_EXPERIMENTS << endl;
  cout << "mean  = " << mean << endl;
  cout << "std. deviation = " << standard_deviation << endl;

  cout << "test rng_wrapper done" << endl;
  
  return 0;
}