C++ TR1：如何使用normal_distribution？

Question

我正在尝试使用C++ STD TechnicalReport1扩展生成遵循正态分布的数字，但是这段代码（改编自此文章）：

mt19937 eng;
eng.seed(SEED);

normal_distribution<double> dist;
// XXX if I use the one below it exits the for loop
// uniform_int<int> dist(1, 52);

for (unsigned int i = 0; i < 1000; ++i) {
  cout << "Generating " << i << "-th value" << endl;
  cout << dist(eng) << endl;
}

仅打印一条“正在生成…”日志消息，然后永远不会退出for循环！ 如果我使用我注释掉的distribution，则它会终止，所以我想知道我做错了什么。 有任何想法吗？

非常感谢！

- puccio

Answer 1

我曾遇到过与最初发布的代码相同的问题，并调查了GNU实现。

首先，有些观察结果如下：
使用g++-4.4和上述代码会停顿；
使用g++-4.5并使用-std=c++0x（即非TR1而是真正的）则以上代码可行。

我的看法是，TR1和c++0x之间在随机数生成和消耗方面的适配器发生了变化——mt19937生成整数，normal_distribution消耗doubles。
c++0x自动适应，而g++ TR1代码则不是这样的。

为了使您的代码在g++-4.4和TR1中可用，请执行以下操作：

std::tr1::mt19937 prng(seed);
std::tr1::normal_distribution<double> normal;
std::tr1::variate_generator<std::tr1::mt19937, std::tr1::normal_distribution<double> > randn(prng,normal);
double r = randn();

- pmeerw

Answer 2

这肯定不会使程序挂起。但是，不确定它是否真正满足您的需求。

 #include <random>
 #include <iostream>

 using namespace std;

 typedef std::tr1::ranlux64_base_01 Myeng; 

 typedef std::tr1::normal_distribution<double> Mydist; 

 int main() 
 { 
      Myeng eng; 
      eng.seed(1000);
      Mydist dist(1,10); 

      dist.reset(); // discard any cached values 
      for (int i = 0; i < 10; i++)
      {
           std::cout << "a random value == " << (int)dist(eng) << std::endl; 
      }

 return (0); 
 }

- Jagannath

---

谢谢，老兄。它运行得非常好，但我想知道为什么它可以用这个引擎运行，而另一个引擎却不行。 - puccio

显然，唯一的区别在于您使用了mt19937数字生成器，而Jagannath使用std :: tr1 :: ranlux64_base_01。从逻辑上讲，我猜测错误可能在您对mt19937对象的实现中（这是我在您之前从未听说过的算法，感谢您 :-)），因为它不是标准库的一部分。 - Stephane Rolland

当绘制随机数时，是否可能将这样的for循环向量化？ 我记得你不能向量化包含函数调用的循环。 - Lindon

Answer 3

如果你的TR1随机数生成实现有缺陷，可以通过编写以下自己的正态分布生成器来避免TR1。

使用任何你信任的随机生成器生成两个均匀（0, 1）随机样本u和v。然后让r = sqrt(-2 log(u))，并返回x = r sin(2 pi v)。（这被称为Box-Mueller方法。）

如果需要均值为mu，标准差为sigma的正态分布样本，则返回sigma * x + mu而不仅仅是x。

- John D. Cook

---

刚试了一下，运行速度超级快——测试了100万个样本，对于1-sigma、2-sigma等范围内的样本，几乎呈现出完美的统计数据。 - Arno Duvenhage

Answer 4

虽然这似乎是一个错误，但快速确认的方法是传递默认的0.0、1.0参数。normal_distribution<double>::normal_distribution()应该等于normal_distribution<double>::normal_distribution(0.0, 1.0)

- MSalters

---

它仍然无法工作，仍然停留在执行第一个计算的状态。 - puccio