pseudo_test.cpp

随机数类,产生随机数

// Examples of random number generation
# include "MersenneTwister.h"
# include "PseudoRandom.h"
# include "matrix.hpp"
#include <fstream>

using namespace std;
using namespace ldas;

int main(void)
{
    int i,j;
      
    // Initialize with the system clock
    PseudoRandom test1;
    
    double a = test1.rand();
    double b = test1.rand();
	
	cout << showpos << left << fixed <<endl;
	
	cout << "Two real numbers in the range [0,1]:  " << a << "," << b << endl;
	
	// Those calls produced the default of floating-point
	// numbers in the range 0 to 1 inclusive.  We can also
	// get integers in the range 0 to 2^32 - 1 (4294967295)
	// with
	
	unsigned long c = test1.randInt();

	cout << "An integer in the range [0," << 0xffffffffUL << "]:  " << c << endl;
	
	// Or get an integer in the range 0 to n (for n < 2^32) with
	
	int d = test1.randInt( 42 );
	
	cout << "An integer in the range [0,42]:  " << d << endl;
	
	// We can get a real number in the range 0 to 1 exclusive of 1 with
	
	double e = test1.randExc();
	
	cout << "A real number in the range [0,1):  " << e << endl;
	
	// The functions rand() and randExc() can also have
	// ranges defined just like randInt()
	
	double f = test1.rand( 2.5 );
	double g = test1.randExc( 10.0 );
	
	cout << "A real numer in the range [0,2.5]:  " << f << endl;
	cout << "And one in the range [0,10.0):  " << g << endl;
    
    // Random number generators need a seed value to start
	// producing a sequence of random numbers.  We gave no
	// seed in our declaration of mtrand1, so one was
	// automatically generated from the system clock (or the
	// operating system's random number pool if available).
	// Alternatively we could provide our own seed. Each
	// seed uniquely determines the sequence of numbers that
	// will be produced.  We can replicate a sequence by
	// starting another generator with the same seed.
	
	PseudoRandom test2a( 5489 );   // makes new PseudoRandom with given seed
    
	double h1 = test2a.rand();   // gets the first number generated
	
	PseudoRandom test2b( 5489 );  // makes an identical PseudoRandom
	
	double h2 = test2b.rand();   // and gets the same number
	
	cout << "These two numbers are the same:  " << h1 << ", " << h2 << endl;
	
	// We can also restart an existing PseudoRandom with a new
	// seed
	
	test2a.seed( 1776 );
	test2b.seed( 1941 );
	
	double i1 = test2a.rand();
	double i2 = test2b.rand();
	
	cout << "Re-seeding gives different numbers:  " << i1 << ", " << i2 << endl;
    
    // But there are only 2^32 possible seeds when we pass a
	// single 32-bit integer.  Since the seed dictates the
	// sequence, only 2^32 different random number sequences
	// will result.  For applications like Monte Carlo
	// simulation we might want many more.  We can seed with
	// an array of values rather than a single integer to
	// access the full 2^19937-1 possible sequences.
	
	unsigned long seed[ MTRand::N ];
	for( int s = 0; s < MTRand::N; ++s )
	{
		seed[s] = 23 * s;  // fill with anything
	}
	PseudoRandom test3( seed );
	
	double j1 = test3();
	double j2 = test3();
	double j3 = test3();
	
	cout << "We seeded this sequence with 19968 bits:  " << j1 << ", " << j2 << ", " << j3 << endl;
	
	cout << endl;
	
	fstream gauss("gauss.txt",ios::out);
	if (!gauss)
	{ 
  		cout<< "MAIN: Error while opening output file "<< "gauss.txt" <<endl; 
  		return 1;
	}
	cout << "Test gaussian with mean 0" << endl;
    for (i=0; i<10000; i++) 
    {
		//cout << setw(12) << test1.gaussian(0,1);
		//if (i%5==4) cout << endl;
		double randomnumber = test1.poisson(4.0);
		gauss << randomnumber << endl;
    }
	gauss.close();
	
   
    // Generate multivariate gaussian distribution samples
    double buf1[2]={0.0,5.0};// 均值向量
    double buf2[4]={1.0,0.45,0.45,0.25};//协方差矩阵
    Matrix<double> m(buf1,2,1),sigma(buf2,2,2),mnorm(2,100);
    mnorm = test1.mnormrnd(100,m,sigma);
    for (j=1;j<=100;j++)
    	cout << mnorm(1,j) << " " << mnorm(2,j) << endl;
    
    cout << "test truncated multivariate normal distribution" << endl;
    cout << "first:between   -0.5 and  0.5" <<endl;
    cout << "second:between 4.5 and 5.5" <<endl;
    cout << endl;
    int buf3[3]={1,2};
    double buf4[3]={0.0,5.0}; // 均值向量
    double buf5[9]={1.0,0.45,0.45,0.25};//协方差矩阵
    double buf6[9]={1.0,0.0,0.0,1.0};//必须是非奇异的
    double buf7[3]={-0.5,4.5};//下边界
    double buf8[6]={0.5,5.5};//上边界
    Matrix<int> ks(buf3,2,1);
    Matrix<double> mu(buf4,2,1), a2(buf7,2,1),b2(buf8,2,1),la(2,1),lb(2,1);
    Matrix<double> sigma2(buf5,2,2),d2(buf6,2,2),xd(2,100);
	
	//xd 包含所生成的随机数 subject to a < d*x < b
	xd = test1.tnorm_rnd(100,mu,sigma2,a2,b2,la,lb,d2,ks);
	for (i=1;i<=100;i++)
	{
		cout << xd(1,i) <<" "<< xd(2,i)<< endl;
	}
}