输出导致蒙特卡罗积分的不同结果

Question

我最近偶然发现了一个C++错误/功能,我无法完全理解,希望有更好的C++知识的人能指出我正确的方向.

下面你会发现我尝试使用蒙特卡罗积分找出高斯曲线下的区域.食谱是:

1. 生成大量正态分布的随机变量(平均值为零,标准差为1).
2. 将这些数字平方.
3. 取所有正方形的平均值.平均值将是对曲线下面积的非常接近的估计(在高斯的情况下,它是1.0).

下面的代码由两个简单的函数组成:rand_uni它返回一个随机变量,均匀分布在0和1之间,并且rand_norm,它是一个(相当差,但"对于政府工作来说足够好")近似正态分布的随机变量.

main通过一个循环运行十亿次,rand_norm每次调用,将其平方pow并添加到累积变量.在此循环之后,累计结果仅除以运行次数并打印到终端Result=<SOME NUMBER>.

问题在于以下代码的非常古怪的行为:当每个生成的随机变量被打印到cout(是,十亿次)时,无论使用何种编译器,最终结果都是正确的(1.0015,这与我的非常接近)想).如果我没有在每次循环迭代打印随机变量,我得到inf下gcc和448314下clang.

坦率地说,这只是博格尔斯头脑和,因为它是一个用C++我第一次(-ish)的遭遇,我真的不知道,这个问题可能是什么:它是用的东西pow？是否cout采取行动怪异？

任何暗示都将非常感激!

来源

// Monte Carlo integration of the Gaussian curve

#include <iostream>
#include <cstdlib>
#include <cmath>


using namespace std;


enum {
  no_of_runs = 1000000
};


// uniform random variable
double rand_uni() {
  return ((double) rand() / (RAND_MAX));
};


// approximation of a normaly distributed random variable
double rand_norm() {
  double result;
  for(int i=12; i > 0; --i) {
    result += rand_uni();
  }
  return result - 6;
};


int main(const int argc,
         const char** argv) {

  double result = 0;
  double x;

  for (long i=no_of_runs; i > 0; --i) {
    x = pow(rand_norm(), 2);
    #ifdef DO_THE_WEIRD_THING
    cout << x << endl;      // MAGIC?!
    #endif
    result += x;
  }

  // Prints the end result
  cout << "Result=" 
       << result / no_of_runs
       << endl << endl;

}

Makefile

CLANG=clang++
GCC=g++
OUT=normal_mc

default: *.cpp
    $(CLANG) -o $(OUT).clang.a *.cpp
    $(CLANG) -o $(OUT).clang.b -DDO_THE_WEIRD_THING *.cpp
    $(GCC) -o $(OUT).gcc.a *.cpp
    $(GCC) -o $(OUT).gcc.b -DDO_THE_WEIRD_THING *.cpp

Answer 1

result在开始之前，in未double rand_norm()初始化为 0 += 随机值。

所有 cout 的使用并重置堆栈内存的某些部分，该部分稍后由 rand_norm 的调用使用，在执行 cout 时“解决”您的问题。

将第一行更改为double result = 0.0;以修复它。

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

另外，您应该为随机数生成器播种，如果没有它，您将始终获得完全相同的伪随机数序列，添加一个

srand(time(NULL));

在第一次调用 rand() 之前，或者查看一些更好的随机数生成器，例如Boost.Random