0

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.

4

2 回答 2

3

您似乎mean以错误的方式计算。mean应该平均超过N，而你只对所有数组元素求和。当前mean实际上是sum。

mean = mean /N
于 2013-05-08T00:36:27.480 回答
2

rand()在大多数实现中是一个质量非常低的随机数生成器。一些 Linux 版本会从内核熵池中获取价值，但不能保证跨平台（例如在 Windows 上？）使用 Mersenne Twister 代替。Boost库实现了一个。

编辑： taocp 答案突出了编码问题，但 RNG 问题仍然适用。

于 2013-05-08T00:36:12.790 回答