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

如何在 C 或 C + + 中按照正态分布轻松地生成随机数?

I don't want any use of Boost.

I know that Knuth talks about this at length but I don't have his books at hand right now.

219005 次浏览

一个快速简单的方法就是把一些均匀分布的随机数相加,然后取它们的平均值。请参阅 中心极限定理以获得关于此工作原理的完整解释。

使用 std::tr1::normal_distribution

Tr1名称空间不是 Boost 的一部分。这个名称空间包含了 C++ Technical Report 1中添加的库,可以在最新的微软编译器和 gcc 中使用,独立于升级。

从正则 RNG 生成高斯分布数有许多方法。

Box-Muller 变换是常用的。它正确地产生正态分布的值。数学很简单。你生成两个(统一的)随机数,通过对它们应用一个公式,你得到两个正态分布的随机数。返回一个,另一个保存到下一个随机数的请求中。

C + + 11

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

C or older C++

以下是一些按复杂性递增顺序排列的解决方案:

  1. Add 12 uniform random numbers from 0 to 1 and subtract 6. This will match mean and standard deviation of a normal variable. An obvious drawback is that the range is limited to ±6 – unlike a true normal distribution.

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

  3. 为了达到最好的精度,我建议绘制制服并应用反累积正态分布来得到正态分布的变量。给你是一种非常好的反累积正态分布算法。

1973年《应用统计学》 ,第22期,第92-97页,“关于乘同余伪随机数生成器的 Box-Muller 变换的应用”

You can use the GSL. Some 给出了完整的例子 to demonstrate how to use it.

下面是一个基于一些引用的 C + + 示例。这是快速和肮脏的,你最好不要重新发明和使用提高库。

#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 ie 总计数的比例。多少个样本,得到 z 值并绘制出来。一条向上的直线是理想的结果)。

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

我遵循了 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);
}

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

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

如果你使用 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);

您可以使用许多其他发行版来转换随机数引擎的输出。

看看我找到了什么。

这个 图书馆使用了金字塔算法。

我创建了一个 用于正态分布随机数生成基准测试的 C + + 开源项目

它比较了几种算法,包括

  • 中心极限定理法
  • Box-Muller 变换
  • 马尔萨格利亚极性法
  • 金字塔算法
  • 逆变换采样法。
  • cpp11random使用 C + + 11 std::normal_distributionstd::minstd_rand(实际上是 Clang 语中的 Box-Muller 变换)。

IMac Corei5-3330S@2.70 GHz,clang 6.1,64位:

normaldistf

为了准确起见,程序验证了样本的平均值、标准差、偏度和峰度。结果表明,4、8、16个统一数相加的 CLT 法不如其他方法具有较好的峰度。

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

因此,我将建议使用 Box-Muller 作为 SIMD 指令集的体系结构,否则可以使用 zigurat。


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

计算机是确定性设备,计算不存在随机性。 此外,CPU 中的算术设备可以计算某些有限整数集(在有限域中进行计算)和有限实有理数集上的和。并执行位操作。数学需要更多像[0.0,1.0]这样的优秀集合,这些集合的分数是无限的。

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

有一种算法叫做伪随机生成器。因为我觉得伪随机生成器的目的是为了模拟随机性。好的标准是: 经验分布是收敛的(在某种意义上,点,均匀,L2)到理论 从随机生成器获得的值似乎是独立的。当然,从“真实的观点”来看,这是不正确的,但我们假设它是正确的。

One of the popular method - you can summ 12 i.r.v with uniform distributions....But to be honest during derivation Central Limit Theorem with helping of Fourier Transform, Taylor Series, it is neededed to have n->+inf assumptions couple times. 举个理论上的例子,我个人不明白人们如何在均匀分布的情况下,进行12个静脉注射的总和。

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


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: 十六概率论,范例6.12,第246-247页”,克里先科 · 亚历山大 · 彼得罗维奇 ISBN 5-7038-2485-0

不幸的是,我不知道这本书是否存在翻译成英文的情况。

ABc0分享了三种不同的方法,可以用正态分布轻松地生成随机数。

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

反累积正态分布的计算方法有很多种,其中最常用的一种是在 http://chasethedevil.github.io/post/monte-carlo-inverse-cumulative-normal-distribution/上进行测试

In my opinion, there is not much incentive to use something else than algorithm AS241 from Wichura: it is machine precision, reliable and fast. Bottlenecks are rarely in the Gaussian random number generation.

这里的最佳答案是提倡 Box-Müller,你应该意识到它有已知的缺陷,我引用 https://www.sciencedirect.com/science/article/pii/S0895717710005935:

在文献中,博克斯-穆勒有时被认为稍逊一筹,主要有两个原因。首先,如果将 Box-Muller 方法应用于来自坏线性同余方法的数字,转换后的数字对空间的覆盖极其有限。在许多书中都可以找到带有螺旋尾巴的变换数的情节,最著名的是在瑞普利的经典著作中,他可能是第一个做出这种观察的人

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;
}

1)从图形上直观地生成高斯随机数的方法是使用类似于蒙特卡罗方法的东西。你可以用 C 语言中的伪随机数生成器在高斯曲线周围的盒子里生成一个随机点。你可以利用分布方程计算出这个点是在正态分布的内部还是下面。如果这个点在正态分布里面,那么你就得到了高斯随机数作为这个点的 x 值。

这个方法并不完美,因为从技术上讲,高斯曲线是朝向无穷大的,你不能创建一个在 x 维度上接近无穷大的盒子。但是高斯曲线在 y 维上接近0的速度很快,所以我不会担心这个。C 语言中变量大小的限制可能是影响准确性的一个限制因素。

2)另一种方法是使用中心极限定理,当独立随机变量被加入时,它们形成一个正态分布。记住这个定理,你可以通过添加大量的独立随机变量来近似高斯随机数。

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

蒙特卡罗方法 最直观的方法是使用蒙特卡罗方法。取一个合适的范围-X,+ X。 X 的值越大,正态分布越精确,但收敛需要更长的时间。 在 X 到 X 之间选择一个随机数 Z。 B 保持正态分布的概率为 N(z, mean, variance),其中 n 为。否则放弃,回到步骤(a)。