Random Number Failure


Posted
Comments None

I was recently starting to write a monte carlo simulation involving part survival, and got some very strange results. For one single time-step, 30 time-steps after the start of the simulation, parts broke less often. I made the probability of failure the same for every time-step and got the following plot. Note the large downward spike at 30 in the otherwise smooth decaying exponential curve.

example of random number generator failure

Since I was still writing code during the quick-and-dirty phase, just to get something started, I was using the POSIX system call random() to generate my random numbers. Now I have always known that the system-provided random numbers are awful, but I had never actually encountered such a striking failure. Here is the code after it was hacked down to eliminate errors on my part:


#include <cstdio>
#include <cstdlib>

int
main()
    {
    const int seed = 1 ;
    const double p_m = 0.05 ;
    const int NDIST = 200 ;
    const int NRUN = 1000000 ;

    srandom(seed) ;

    static int dist_m[NDIST] ;
    for( int d = 0 ; d < NDIST ; ++d )
	dist_m[d] = 0 ;

    for( int run = 0 ; run < NRUN ; ++run )
	{
	int M = 0 ;
	for(;;)
	    {
	    double r_m = double( random() ) / double( RAND_MAX ) ;
	    if( r_m < p_m )
		break ;
	    else
		M++ ;
	    }
	dist_m[ M < NDIST ? M : NDIST-1 ]++ ;
	}
    for( int d = 0 ; d < NDIST ; ++d )
	{
	if( dist_m[d] )
	    printf( "%3d %f\n", d, dist_m[d] / double(NRUN) ) ;
	}
    }

Replacing the POSIX random() call with something better fixed the problem nicely.


Categories Numerical Analysis, Probability

Comments

There are currently no comments on this article.

Please leave a comment

Enter your comment below. Fields marked * are required. You must preview your comment first before finally posting.





← Older Newer →