not random enough for monte carlo

Question

I am trying to generate values from a normal distribution using a monte carlo method, as per the website http://math60082.blogspot.ca/2013/03/c-coding-random-numbers-and-monte-carlo.html

I modified the code a bit from the original so it calculates the variance and mean for the numbers generated directly to check if the method is working rather than do the tests separately (same difference really but just a heads up).

Question

Regardless of what I do, the variance is way above 1 and the mean is not zero. Is it possible the pseudo-random numbers generated aren't random enough?

Code

PLEASE NOTE THAT THE AUTHOR OF THE ABOVE GIVEN WEBSITE IS THE PERSON WHO WROTE THE CODE

#include <cstdlib>
#include <cmath>
#include <ctime>
#include <iostream>
using namespace std;
// return a uniformly distributed random number
double uniformRandom()
{
    return ( (double)(rand()) + 1. )/( (double)(RAND_MAX) + 1. );
}

// return a normally distributed random number
double normalRandom()
{
    double u1=uniformRandom();
    double u2=uniformRandom();
    return cos(8.*atan(1.)*u2)*sqrt(-2.*log(u1)); 
}

int main()
{
    double z;
    int N=1000;
    double array[N];
    double mean=0 ,variance=0;
    srand(time(NULL));

    for(int i=0;i<N;i++)
    {
        z=normalRandom();
        cout << i << "->"<< z<< endl;
        mean+=z;
        array[i]=z;
    }

    mean=mean/N ;
    cout << " mean = " << mean << endl;

    for(int i=0;i<N;i++)
    {
        variance = variance + (mean - array[i])*(mean - array[i]);
    }
    variance = variance/N;
    cout << " variance = " << variance << endl;

    return 0;
}

UPDATE

Apparently as pointed by users, I screwed up and the program was not working because of a very silly mistake.

I've read in the past that generating a normal distribution random from a uniform random is generally not particularly accurate -- better to use a normal distribution RNG. (Dunno exactly where to find one, though. Used to be one in the old standard Fortran libs, but don't know about C++.) — Hot Licks, May 08 '13 at 00:41
@HotLicks: I don't know where you read that but that's a very shaky statement. Generation of normal distributions form a uniform distribution is quite a [common](http://en.wikipedia.org/wiki/Marsaglia_polar_method) [thing](http://en.wikipedia.org/wiki/Box%E2%80%93Muller_transform). *However*, using a poor uniform generator (like `rand()`) will lead to a poor normal distribution. — Mike Bailey, May 08 '13 at 00:42
@HotLicks Noted, apparently the friend I am doing this for, wants it to be a monte carlo method. — metric-space, May 08 '13 at 00:44
@MikeBantegui -- Like I said, this was "in the past" -- several decades ago, I'm thinking. — Hot Licks, May 08 '13 at 01:19

score 3 · Accepted Answer · answered May 08 '13 at 00:36

3

You seems computed the mean in a wrong way. mean should be averaged over N, while you only sum over all array elements. current mean is actually sum.

mean = mean /N

answered May 08 '13 at 00:36

taocp

23,276
10
49
62

thanks for pointing that out.BUt the variance is still nowhere near 1. – metric-space May 08 '13 at 00:38
2

@nerorevenge you may try to generate more than 2000 rand numbers and see whether the variance is getting closer to 1? the rand() may not be a good random number generator (just my wild guess). – taocp May 08 '13 at 00:40
1

@nerorevenge: you're wrong. the incorrect value of the mean affects the value of the variance. after you fix the value of the mean, that fixes all the problems in your code. – carlosdc May 08 '13 at 00:46
@nerorevenge Standard deviation should be 1, not Variance. Std. Dev = square root of variance. – Adrian Torrie May 08 '13 at 01:01
@iValueValue does it matter? (not for just this case, but for any case) – metric-space May 08 '13 at 01:06
1

@iValueValue: By definition, if `standard deviation = 1` then `variance = 1` because `sqrt(1) = 1` and `1 * 1 = 1`. – Mike Bailey May 08 '13 at 01:07

Stefano Sanfilippo · Answer 2 · 2013-05-08T00:41:14.057

2

rand() is a very low quality random numbers generator in most implementations. Some Linux versions would take value from kernel entropy pool, but it is not guaranteed across platforms (e.g. on Windows?) Use a Mersenne Twister instead. Boost libraries implement one.

EDIT: taocp answer highlights a coding problem, but the RNG issue still applies.

edited May 08 '13 at 00:41

answered May 08 '13 at 00:36

Stefano Sanfilippo

32,265
7
79
80

Noted, and the code is being run on a linux machine ,so your answer could be right, let me get gback to you. – metric-space May 08 '13 at 00:40
Gback? Anyway, not all Linux libc do actually use the entropy pool. I cannot tell which do and which do not, but I'm pretty sure `rand` is good for playing "guess my number" against the PC, not for doing statistical analysis or cryptography. – Stefano Sanfilippo May 08 '13 at 00:45
1

You don't need boost for this if you're using an up to date version of gcc. `` has all of this functionality now. – Yuushi May 08 '13 at 02:20

not random enough for monte carlo

2 Answers2