logo标签列表

在C/C++中生成遵循正态分布的随机数

c++crandomdistributionnormal-distribution
128

如何在C或C ++中轻松生成正态分布的随机数？

我不想使用Boost。

我知道Knuth对此进行了详细讨论，但我现在手头没有他的书。

2
重复了 https://dev59.com/YHVD5IYBdhLWcg3wI3-L 和 https://dev59.com/pUjSa4cB1Zd3GeqPEl75 中的一个。 - dmckee --- ex-moderator kitten
18个回答
100

有许多方法可以从常规随机数生成器中生成高斯分布的数字。

Box-Muller变换是常用的方法。它可以正确地产生具有正态分布的值。这个数学公式很简单。您可以生成两个（均匀）随机数，通过对它们应用一个公式，您可以获得两个正态分布的随机数。返回一个，将另一个保存下来以供下一次请求随机数时使用。

10
如果你需要速度，那么极坐标法会更快。而Ziggurat算法甚至更快（尽管要写的复杂得多）。 - Joey
3
我找到了Ziggurat的一个实现，链接在这里：http://people.sc.fsu.edu/~jburkardt/c_src/ziggurat/ziggurat.html 它非常完整。 - dwbrito
29
注意，C++11增加了std::normal_distribution，它可以完全满足您的要求，而不需要深入研究数学细节。 - user283145
3
std::normal_distribution不能保证在所有平台上都是一致的。我正在进行测试，MSVC提供的一组值与例如Clang提供的不同。C++11引擎似乎会生成相同的序列（给定相同的种子），但是C++11分布似乎在不同平台上采用了不同的算法实现。 - Arno Duvenhage
https://people.sc.fsu.edu/~jburkardt/c_src/normal/normal.c谢谢 @dwbrito - angstyloop
54

C++11

C++11提供了std::normal_distribution，这是我今天会选择的方法。

C或早期的C++

以下是按升序排列的一些解决方案：

  1. 将0到1之间的12个均匀随机数相加并减去6。这将匹配正态变量的平均值和标准差。一个明显的缺点是范围仅限于±6，而不像真正的正态分布。

  2. Box-Muller变换。这在上面已经列出，并且相对容易实现。然而，如果需要非常精确的样本，请注意Box-Muller变换与一些均匀生成器结合使用会受到称为Neave Effect1的异常影响。

  3. 为了获得最佳精度，建议绘制均匀分布并应用逆累积正态分布来获得正态分布变量。这里有一个非常好的逆累积正态分布算法。

1. H. R. Neave, “On using the Box-Muller transformation with multiplicative congruential pseudorandom number generators,” Applied Statistics, 22, 92-97, 1973

3
原始参考文献已添加。有趣的一点是，在谷歌搜索“box muller neave”以找到该参考文献时，这个StackOverflow问题出现在第一页的搜索结果中！ - Peter G.
1
是的，它并不为某些小社区和兴趣群体以外的人所熟知。 - pyCthon
@Peter G. 为什么会有人踩你的回答呢？可能是同一个人也踩了我下面的评论，但我认为你的回答非常好。如果 Stack Overflow 要求踩的人必须留下真实的评论就好了。我猜测大多数对旧话题的踩都只是无聊和恶意的行为。 - Pete855217
将0-1之间的12个均匀数相加并减去6，这个变量的分布将是正态分布吗？您能否提供一个推导的链接，因为在推导中极限定理的中心极限定理，n->+inf非常需要假设。 - Konstantin Burlachenko
1
@bruzias，将0到1之间的12个均匀数相加并进行减法运算并不服从正态分布，它只会得到相同的平均值（即0）和标准差（即sqrt(12*1/12)=1，因为0到1的均匀分布的标准差为1/12，而独立变量的方差相加）。 - Peter G.
显示剩余3条评论
31

一种快速而简单的方法是将一些均匀分布的随机数相加，并取它们的平均值。请参阅中心极限定理，了解这种方法有效的完整解释。

4
样本数越大，平均值越接近于高斯分布。如果您的应用程序对分布的精确性有严格要求，那么最好使用一些更严格的方法，比如Box-Muller，但对于许多应用程序（例如生成音频应用程序的白噪声），您可以通过相对较少的平均样本数量（例如16个）来实现。 - Paul R
3
另外，如何参数化以获得特定的方差，比如说您想要均值为10，标准差为1？ - Morlock
1
这是一种非常低效的从正态分布中生成样本的方法。我绝对不会称其为“快速”。 - Petter
1
@Ben：你能给我指一下一个高效的算法吗？我只用过平均技术来生成大致符合高斯分布的噪声，用于音频和图像处理，并且有实时约束 - 如果有一种方法可以在更少的时钟周期内实现这一点，那将非常有用。 - Paul R
1
@Petter：在浮点数值的一般情况下，你可能是对的。但是在像音频这样的应用领域中，您需要快速的整数（或定点）高斯噪声，并且精度并不太重要，简单的平均方法更有效和有用（特别是对于嵌入式应用程序，在那里甚至可能没有硬件浮点支持）。 - Paul R
显示剩余9条评论
28

我为正态分布随机数生成创建了一个C++开源项目，它比较了几种算法，包括：

  • 中心极限定理方法
  • Box-Muller变换
  • Marsaglia polar方法
  • Ziggurat算法
  • 反转换抽样方法。
  • cpp11random使用C++11 std::normal_distribution与std::minstd_rand（实际上在clang中是Box-Muller变换）。

在iMac Corei5-3330S@2.70GHz , clang 6.1，64位单精度（float）版本的结果如下:

normaldistf

为了检验正确性，程序验证了样本的平均值、标准差、偏度和峰度。发现通过对4个、8个或16个均匀数求和的CLT方法在峰度方面不如其他方法表现好。

Ziggurat算法的性能比其他方法更好。然而，它不适用于SIMD并行处理，因为需要进行表查找和分支操作。使用带有SIMD指令集的Box-Muller比非SIMD版本的ziggurat算法快得多（x1.79，x2.99）。

因此，我建议在具有SIMD指令集的架构中使用Box-Muller，否则可以考虑使用ziggurat。

P.S.基准测试使用最简单的LCG PRNG生成均匀分布的随机数。因此，对于某些应用程序而言，这可能不足够。但是，由于所有实现都使用相同的PRNG，因此性能比较应该是公平的，因此基准测试主要测试了变换的性能。

2
“但是性能比较应该是公平的，因为所有实现都使用相同的伪随机数生成器。” 除了BM使用一个输入RN每个输出，而CLT使用更多等等... 因此，生成均匀随机数所需的时间很重要。 - greggo
14

这是一个基于一些参考文献的C++示例。这只是个快速而简单的实现，最好不要重新发明轮子，使用boost库更好。

#include "math.h" // for RAND, and rand
double sampleNormal() {
    double u = ((double) rand() / (RAND_MAX)) * 2 - 1;
    double v = ((double) rand() / (RAND_MAX)) * 2 - 1;
    double r = u * u + v * v;
    if (r == 0 || r > 1) return sampleNormal();
    double c = sqrt(-2 * log(r) / r);
    return u * c;
}
您可以使用Q-Q图来检查结果，看它如何逼近真正的正态分布（将样本排名为1..x，将排名变成总x计数的比例即样本数量，获取z值并将其绘制出来。向上的直线是期望的结果）。
1
sampleNormalManual()是什么？ - solvingPuzzles
1
这个程序注定会在某些罕见事件中崩溃（向老板展示应用程序时有所警示吗？）。这应该使用循环实现，而不是使用递归。这个方法看起来很陌生。它的来源/如何调用是什么？ - the swine
Box-Muller从Java实现转录而来。正如我所说，这只是一个快速而简单的实现，随意修改即可。 - Pete855217
看起来非常合理，而且效果很好。我在JavaScript中实现了这个（并测试了分布）。 - Octopus
2
就此而言，许多编译器将能够将特定的递归调用转换为“跳转到函数顶部”。问题是你是否想依赖它 :-)此外，迭代次数大于10的概率为480万分之1。p(>20)是其平方，以此类推。 - greggo
显示剩余2条评论
13

这是现代C++编译器生成样本的方法。

#include <random>
...
std::mt19937 generator;
double mean = 0.0;
double stddev  = 1.0;
std::normal_distribution<double> normal(mean, stddev);
cerr << "Normal: " << normal(generator) << endl;
生成器应该真正被种子化。 - Walter
2
它总是被种子化的。有一个默认的种子。 - Petter
12

使用std::tr1::normal_distribution。

std::tr1命名空间不是boost的一部分。它包含来自C++技术报告1的库添加，并可在最新的Microsoft编译器和gcc中独立于boost使用。

25
他没有要求标准，他要求的是“不增压”。 - JoeG
6
为什么这个答案没有得到很多赞？使用GSL有哪些需要注意的缺陷吗？ - noskillnoob
4

如果你使用的是C++11，你可以使用std::normal_distribution：

#include <random>

std::default_random_engine generator;
std::normal_distribution<double> distribution(/*mean=*/0.0, /*stddev=*/1.0);

double randomNumber = distribution(generator);

除了这个随机数引擎，还有许多其他发行版可以用来转换其输出。

那已经被Ben提到了（https://dev59.com/1XE95IYBdhLWcg3wadGq#11977979） - Mat
4
但是它需要很多随机数。 - gkucmierz
网页内容由stack overflow 提供, 点击上面的
可以查看英文原文，
原文链接