Different Results of Monte Carlo Integration Due To Output

Question

I've just recently stumbled upon a C++ bug/feature, that I can't fully understand and was hoping someone with a better knowledge of C++ here could point me in the right direction. Below you will find my attempt at finding out the area under the Gaussian curve using Monte Carlo integration. The recipe is: - Generate a large number of normally distributed random variables (with mean of zero and standard deviation of one). - Square these numbers. - Take the average of all the squares. The average will be a very close estimation of the area under the curve (in the case of the Gaussian, it is 1.0). The code below consists of two simple functions: `rand_uni`, which returns a random variable, uniformly distributed between zero and one, and `rand_norm`, which is a (rather poor, but 'good enough for government work') approximation of a normally distributed random variable. `main` runs through a loop one billion times, calling `rand_norm` each time, squaring it with `pow` and adding to the accu…

Accepted Answer

The `result` in `double rand_norm()` is not initialised to 0 before you start += random values to it.

All the cout's use and reset some part of the stack memory that is later used by the call of rand_norm, "fixing" your problem when you do the cout's.

Change the first line to `double result = 0.0;` to fix it.

```
double rand_norm() {
  double result = 0.0;
  for(int i=12; i > 0; --i) {
    result += rand_uni();
  }
  return result - 6;
};
```

Also, you should seed your random number generator, you will always get the exact same sequence of pseudo random numbers without it, add a

```
srand(time(NULL));
```

before your first call to rand(), or look at some better random number generators, e.g. Boost.Random

Different Results of Monte Carlo Integration Due To Output

Solution

Problem