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

Posted 2023-02-25

【中文标题】在 C/C++ 中按照正态分布生成随机数

【英文标题】：Generate random numbers following a normal distribution in C/C++

【发布时间】：2011-01-20 11:51:55

【问题描述】：

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

我不想使用任何 Boost。

我知道 Knuth 详细地谈到了这一点，但我现在手头没有他的书。

【问题讨论】：

***.com/questions/75677/… 和 ***.com/questions/1109446/… 中的一个或另一个的重复

【参考方案1】：

逆累积正态分布存在多种算法。在http://chasethedevil.github.io/post/monte-carlo-inverse-cumulative-normal-distribution/上测试最流行的量化金融

在我看来，除了来自Wichura 的算法 AS241 之外，没有太多动机使用其他东西：它具有机器精度、可靠和快速。高斯随机数生成中很少出现瓶颈。

这里的最佳答案提倡 Box-Müller，您应该知道它存在已知的缺陷。我引用https://www.sciencedirect.com/science/article/pii/S0895717710005935：

在文献中，Box-Muller 有时被认为略逊一筹，主要有两个原因。首先，如果将 Box-Muller 方法应用于来自不良线性同余生成器的数字，则转换后的数字提供的空间覆盖率极差。在许多书中都可以找到带有螺旋尾的变换数字图，最著名的是里普利的经典著作，他可能是第一个做出这种观察的人”

【参考方案2】：

蒙特卡洛方法 最直观的方法是使用蒙特卡罗方法。取一个合适的范围-X，+X。较大的 X 值将导致更准确的正态分布，但需要更长的时间才能收敛。 一种。在 -X 到 X 之间选择一个随机数 z。 湾。保持N(z, mean, variance) 的概率，其中 N 是高斯分布。否则放弃并返回步骤 (a)。

【参考方案3】：

C++11

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

C 或更早的 C++

这里有一些解决方案，按复杂度升序排列：

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

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

为了获得最佳精度，我建议绘制制服并应用逆累积正态分布来获得正态分布变量。 Here 是一个非常好的逆累积正态分布算法。

^{1. H. R. Neave，“关于将 Box-Muller 变换与乘法同余伪随机数生成器结合使用”，应用统计，1973 年 22 月 92-97 日}

【讨论】：

您是否有另一个指向有关 Neave 效果的 pdf 的链接？还是原始期刊文章参考？谢谢

@stony***nick 添加了原始参考。很酷的评论：在谷歌搜索“box muller neave”以查找参考时，这个 *** 问题出现在第一个结果页面上！

是的，在某些小社区和利益集团之外，它并不是每个人都知名

@Peter G. 为什么有人会否决你的答案？ - 可能同一个人也在下面发表了我的评论，我很好，但我认为你的回答非常好。如果 SO 让反对票强制发表真正的评论，那就太好了。我怀疑大多数反对老话题的投票都是轻浮和无稽之谈。

"从 0 到 1 的 12 个统一数字相加并减 6。" ——这个变量的分布会有正态分布吗？你能提供一个推导的链接吗，因为在推导中心极限定理时，n->+inf 是非常需要假设的。

【参考方案4】：

generate Gaussian-distributed numbers from a regular RNG有很多方法。

Box-Muller transform 是常用的。它正确地产生具有正态分布的值。数学很容易。您生成两个（均匀）随机数，并通过对它们应用公式，您得到两个正态分布的随机数。返回一个，并保存另一个用于下一个随机数请求。

【讨论】：

如果你需要速度，那么极地方法会更快。 Ziggurat 算法甚至更多（尽管编写起来要复杂得多）。

在这里people.sc.fsu.edu/~jburkardt/c_src/ziggurat/ziggurat.html找到了一个Ziggurat的实现，已经很完整了。

注意，C++11 添加了std::normal_distribution，它完全符合您的要求，无需深入研究数学细节。

std::normal_distribution 不能保证在所有平台上都是一致的。我现在正在做测试，MSVC 提供了一组不同的值，例如 Clang。 C++11 引擎似乎生成相同的序列（给定相同的种子），但 C++11 发行版似乎是在不同平台上使用不同算法实现的。

【参考方案5】：

1) 生成高斯随机数的图形直观方法是使用类似于蒙特卡洛方法的方法。您将使用 C 中的伪随机数生成器在高斯曲线周围的框中生成一个随机点。您可以使用分布方程计算该点是在高斯分布内部还是之下。如果该点在高斯分布内，那么您将高斯随机数作为该点的 x 值。

此方法并不完美，因为从技术上讲，高斯曲线趋向无穷大，而您无法创建在 x 维度上接近无穷大的框。但是高斯曲线在 y 维度上非常快地接近 0，所以我不会担心。 C 中变量大小的限制可能更多地限制了您的准确性。

2) 另一种方法是使用中心极限定理，该定理指出，当添加独立随机变量时，它们会形成正态分布。牢记这个定理，您可以通过添加大量独立随机变量来近似高斯随机数。

这些方法不是最实用的，但是当您不想使用预先存在的库时，这是可以预料的。请记住，此答案来自很少或没有微积分或统计经验的人。

【参考方案6】：

Box-Muller 实现：

#include <cstdlib>
#include <cmath>
#include <ctime>
#include <iostream>
using namespace std;
 // return a uniformly distributed random number
double RandomGenerator()

  return ( (double)(rand()) + 1. )/( (double)(RAND_MAX) + 1. );

 // return a normally distributed random number
double normalRandom()

  double y1=RandomGenerator();
  double y2=RandomGenerator();
  return cos(2*3.14*y2)*sqrt(-2.*log(y1));


int main()
double sigma = 82.;
double Mi = 40.;
  for(int i=0;i<100;i++)
double x = normalRandom()*sigma+Mi;
    cout << " x = " << x << endl;
  
  return 0;

【参考方案7】：

计算机是确定性设备。计算中没有随机性。 此外，CPU 中的算术设备可以评估一些有限整数集（在有限域中执行评估）和有限实有理数集的求和。并且还进行了按位运算。数学处理更出色的集合，例如具有无限点数的 [0.0, 1.0]。

你可以用一些控制器来监听计算机内部的一些线，但它会有统一的分布吗？我不知道。但是如果假设它的信号是大量独立随机变量累积值的结果，那么你会得到近似正态分布的随机变量（概率论证明了）

存在称为伪随机生成器的算法。我觉得伪随机生成器的目的是模拟随机性。善良的标准是： - 经验分布收敛（在某种意义上 - 逐点，均匀，L2）到理论 - 您从随机生成器收到的值似乎是独立的。当然，从“真实的观点”来看这不是真的，但我们假设它是真的。

一种流行的方法 - 你可以对 12 个具有均匀分布的 irv 求和....但是说实话，在傅里叶变换、泰勒级数的帮助下推导中心极限定理时，需要有 n->+inf 假设几次。 例如理论上 - 我个人不理解人们如何执行 12 i.r.v. 的总和。分布均匀。

我在大学里学过概率论。特别是对我来说，这只是一个数学问题。在大学里我看到了以下模型：

---

double generateUniform(double a, double b)

  return uniformGen.generateReal(a, b);


double generateRelei(double sigma)

  return sigma * sqrt(-2 * log(1.0 - uniformGen.generateReal(0.0, 1.0 -kEps)));

double generateNorm(double m, double sigma)

  double y2 = generateUniform(0.0, 2 * kPi);
  double y1 = generateRelei(1.0);
  double x1 = y1 * cos(y2);
  return sigma*x1 + m;

---

这种方式只是一个例子，我猜它存在另一种实现方式。

证明它是正确的可以在这本书中找到 “莫斯科，BMSTU，2004：XVI 概率论，示例 6.12，p.246-247”，Krishchenko Alexander Petrovich ISBN 5-7038-2485-0

很遗憾，我不知道有没有将这本书翻译成英文。

【讨论】：

我投了几票。让我知道这里有什么不好的地方？

问题是如何在计算机中生成伪随机数（我知道，这里语言比较松散），不是数学存在的问题。

是的，你是对的。答案是如何基于具有均匀分布的生成器生成具有正态分布的伪随机数。已提供源代码，您可以用任何语言重写。

当然，我认为这个人正在寻找例如“C/C++ 中的数值配方”。顺便说一句，为了补充我们的讨论，最后一本书的作者提供了一些有趣的参考资料，介绍了几个符合“体面”生成器标准的伪随机生成器。

我在这里做了备份：sites.google.com/site/burlachenkok/download

【参考方案8】：

comp.lang.c 常见问题列表分享了三种不同的方法来轻松生成具有高斯分布的随机数。

你可以看看：http://c-faq.com/lib/gaussian.html

【参考方案9】：

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

【讨论】：

generator 真的应该播种。

总是播种的。有一个默认种子。

【参考方案10】：

我创建了一个C++ open source project for normally distributed random number generation benchmark。

它比较了几种算法，包括

中心极限定理法 Box-Muller 变换 Marsaglia 极地法 Ziggurat 算法 逆变换采样方法。 cpp11random 使用 C++11 std::normal_distribution 和 std::minstd_rand（它实际上是 clang 中的 Box-Muller 变换）。

单精度(float)版本在iMac Corei5-3330S@2.70GHz, clang 6.1, 64-bit上的结果：

为了正确性，程序会验证样本的均值、标准差、偏度和峰度。结果表明，将 4、8 或 16 个均匀数相加的 CLT 方法不像其他方法那样具有良好的峰度。

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

因此，我建议将 Box-Muller 用于带有 SIMD 指令集的体系结构，否则可能是 ziggurat。

---

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

【讨论】：

"但是性能比较应该是公平的，因为所有实现都使用相同的 PRNG" .. 除了 BM 每个输出使用一个输入 RN，而 CLT 使用更多，等等......所以时间生成一个统一的随机 # 事项。

【参考方案11】：

我遵循了http://www.mathworks.com/help/stats/normal-distribution.html 中给出的 PDF 的定义并提出了这个：

const double DBL_EPS_COMP = 1 - DBL_EPSILON; // DBL_EPSILON is defined in <limits.h>.
inline double RandU() 
    return DBL_EPSILON + ((double) rand()/RAND_MAX);

inline double RandN2(double mu, double sigma) 
    return mu + (rand()%2 ? -1.0 : 1.0)*sigma*pow(-log(DBL_EPS_COMP*RandU()), 0.5);

inline double RandN() 
    return RandN2(0, 1.0);

这可能不是最好的方法，但很简单。

【讨论】：

-1 不适用于例如RANDN2(0.0, d + 1.0)。宏因此而臭名昭著。

如果 rand() 的 RANDU 返回零，宏将失败，因为 Ln(0) 未定义。

你真的试过这段代码吗？看起来您已经创建了一个生成数字 Rayleigh distributed 的函数。与Box–Muller transform 相比，它们与cos(2*pi*rand/RAND_MAX) 相乘，而您与(rand()%2 ? -1.0 : 1.0) 相乘。

【参考方案12】：

看看我发现了什么。

此library 使用 Ziggurat 算法。

【参考方案13】：

如果你使用的是 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 (***.com/a/11977979/635608) 已经提到了这一点

【参考方案14】：

看看：http://www.cplusplus.com/reference/random/normal_distribution/。这是产生正态分布的最简单方法。

【参考方案15】：

这是一个基于一些参考资料的 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 值并绘制它们。向上的直线是所需的结果）。

【讨论】：

什么是sampleNormalManual()？

@solvingPuzzles - 抱歉，已更正代码。这是一个递归调用。

这肯定会在一些罕见的事件中崩溃（向你的老板展示应用程序敲响了警钟？）。这应该使用循环来实现，而不是使用递归。这个方法看起来很陌生。来源是什么/怎么称呼？

Box-Muller 转录自 java 实现。正如我所说，它又快又脏，请随时修复它。

FWIW，许多编译器将能够将特定的递归调用转换为“跳转到函数顶部”。问题是您是否要指望它:-) 此外，它需要 > 10 次迭代的概率是 480 万分之一。 p(>20) 是那个的平方，等等。

【参考方案16】：

您可以使用GSL。一些complete examples are given 来演示如何使用它。

【参考方案17】：

使用std::tr1::normal_distribution。

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

【讨论】：

他没有要求标准，他要求“不提升”。

【参考方案18】：

一种快速简便的方法是将多个均匀分布的随机数相加并取其平均值。请参阅Central Limit Theorem，了解其工作原理的完整说明。

【讨论】：

@Morlock 平均样本数越大，越接近高斯分布。如果您的应用程序对分布的准确性有严格的要求，那么您最好使用更严格的东西，例如 Box-Muller，但对于许多应用程序，例如为音频应用程序生成白噪声，您可以使用相当少量的平均样本（例如 16 个）。

另外，你如何参数化它以获得一定的方差，比如你想要一个标准差为 1 的平均值 10？

这是一种从正态分布生成样本的非常低效的方法。我绝对不会称之为“快速”。

@Ben：你能给我指出一个有效的算法吗？我只使用过平均技术来为具有实时约束的音频和图像处理生成近似高斯噪声 - 如果有一种方法可以在更少的时钟周期内实现这一点，那么这可能非常有用。

@Petter：在一般情况下，对于浮点值，您可能是正确的。尽管如此，仍然有像音频这样的应用领域，您需要快速整数（或定点）高斯噪声，而准确性不是太重要，简单的平均方法更有效和有用（特别是对于嵌入式应用，甚至可能没有是硬件浮点支持）。