在C/C++中按照正态分布生成随机数

Question

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

我不想使用Boost.

我知道Knuth详细谈论了这个问题,但我现在还没有他的书.

Answer 1

有许多方法可以从常规RNG生成高斯分布数.

该箱穆勒变换是常用的.它正确生成具有正态分布的值.数学很容易.您生成两个(统一)随机数,并通过对它们应用公式,您将获得两个正态分布的随机数.返回一个,并保存另一个以获取随机数的下一个请求.

Answer 2

C++ 11

C++ 11提供的std::normal_distribution,这就是我今天要走的路.

C或更旧的C++

以下是按复杂程度递增的顺序解决方案:

1. 从0到1添加12个均匀随机数并减去6.这将匹配正常变量的均值和标准差.一个明显的缺点是范围限制在±6 - 与真正的正态分布不同.

2. Box-Muller变换.这在上面列出,并且实现起来相对简单.但是,如果您需要非常精确的样本,请注意Box-Muller变换与一些统一的生成器相结合会遭受称为Neave Effect ¹的异常现象.

3. 为了获得最佳精度,我建议绘制制服并应用反向累积正态分布来得到正态分布的变量.这是一个非常好的逆累积正态分布算法.

^{1. HR Neave,"使用具有乘法同余伪随机数发生器的Box-Muller变换",Applied Statistics,22,92-97,1973}

Answer 3

一种快速简便的方法是将一些均匀分布的随机数相加并取其平均值.有关其工作原理的完整说明,请参阅中心极限定理.

Answer 4

我为正态分布随机数生成基准创建了一个C++开源项目.

它比较了几种算法,包括

• 中心极限定理方法
• Box-Muller变换
• Marsaglia极地方法
• Ziggurat算法
• 逆变换采样方法.
• cpp11random使用C++ 11 std::normal_distribution与std::minstd_rand(它实际上是箱穆勒变换在铛).

floatiMac Corei5-3330S@2.70GHz上的单精度()版本的结果,clang 6.1,64位:

为了正确,程序验证样品的平均值,标准偏差,偏度和峰度.结果发现,通过求和4,8或16个均匀数的CLT方法与其他方法一样没有良好的峰度.

Ziggurat算法比其他算法具有更好的性能.但是,它不适合SIMD并行,因为它需要表查找和分支.具有SSE2/AVX指令集的Box-Muller比非SIMD版本的ziggurat算法快得多(x1.79,x2.99).

因此,我建议使用Box-Muller作为具有SIMD指令集的体系结构,否则可能是ziggurat.

PS基准测试使用最简单的LCG PRNG来生成均匀分布的随机数.因此对某些应用程序来说可能还不够.但性能比较应该是公平的,因为所有实现都使用相同的PRNG,因此基准测试主要测试转换的性能.

Answer 5

这是一个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;
}

您可以使用QQ图来检查结果并查看其与实际正态分布的近似程度(对样本进行排序1..x,将等级转换为x总数的比例,即多少样本,得到z值并绘制它们.向上的直线是期望的结果).

Answer 6

使用std::tr1::normal_distribution.

std :: tr1命名空间不是boost的一部分.它是包含C++技术报告1中的库添加的命名空间,可以在最新的Microsoft编译器和gcc中使用,与boost无关.

Answer 7

这是在现代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;

Answer 8

您可以使用GSL。给出了一些完整的例子来演示如何使用它。