如何在 C 或 C++ 中按照正态分布轻松生成随机数?
我不想使用 Boost。
我知道 Knuth 详细地谈到了这一点,但我现在手头没有他的书。
如何在 C 或 C++ 中按照正态分布轻松生成随机数?
我不想使用 Boost。
我知道 Knuth 详细地谈到了这一点,但我现在手头没有他的书。
有许多方法可以从常规 RNG 生成高斯分布的数字。
Box-Muller 变换是常用的。它正确地产生具有正态分布的值。数学很容易。您生成两个(均匀)随机数,并通过对它们应用公式,您得到两个正态分布的随机数。返回一个,并保存另一个用于下一个随机数请求。
C++11 提供std::normal_distribution
,这是我今天要走的路。
以下是一些按复杂度升序排列的解决方案:
从 0 到 1 加上 12 个均匀随机数并减去 6。这将匹配正态变量的均值和标准差。一个明显的缺点是范围被限制在 ±6 - 与真正的正态分布不同。
Box-Muller 变换。这在上面列出,并且实现起来相对简单。但是,如果您需要非常精确的样本,请注意 Box-Muller 变换与一些均匀生成器相结合会遭受称为 Neave 效应1的异常。
为了获得最佳精度,我建议绘制制服并应用逆累积正态分布来得出正态分布变量。这是一个非常好的逆累积正态分布算法。
1. HR Neave,“关于将 Box-Muller 变换与乘法同余伪随机数生成器结合使用”,应用统计学,1973 年 22 月 92-97 日
一种快速简便的方法是将多个均匀分布的随机数相加并取其平均值。请参阅中心极限定理,以获得有关其工作原理的完整解释。
它比较了几种算法,包括
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。
PS 基准测试使用最简单的 LCG PRNG 来生成均匀分布的随机数。因此,对于某些应用程序可能还不够。但性能比较应该是公平的,因为所有实现都使用相同的 PRNG,所以基准测试主要测试转换的性能。
这是一个基于一些参考资料的 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 值并绘制它们。向上的直线是所需的结果)。
This is how you generate the samples on a modern C++ compiler.
#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;
使用std::tr1::normal_distribution
.
std::tr1 命名空间不是 boost 的一部分。它是包含 C++ 技术报告 1 中添加的库的命名空间,可在最新的 Microsoft 编译器和 gcc 中使用,独立于 boost。
看看: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);
您可以使用许多其他分布来转换随机数引擎的输出。
我遵循了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://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 方法应用于来自不良线性同余生成器的数字,则转换后的数字提供的空间覆盖率极差。许多书中都可以找到带有螺旋尾的变换数字图,最著名的是里普利的经典著作,他可能是第一个做出这种观察的人”
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;
}
comp.lang.c 常见问题列表分享了三种不同的方法来轻松生成具有高斯分布的随机数。
蒙特卡罗方法
最直观的方法是使用蒙特卡罗方法。取一个合适的范围-X,+X。较大的 X 值将导致更准确的正态分布,但需要更长的时间才能收敛。一个。在 -X 到 X 之间选择一个随机数z 。 b.保持N(z, mean, variance)
其中 N 是高斯分布的概率。否则放弃并返回步骤(a)。
1) 生成高斯随机数的图形直观方法是使用类似于蒙特卡洛方法的方法。您将使用 C 中的伪随机数生成器在高斯曲线周围的框中生成一个随机点。您可以使用分布方程计算该点是在高斯分布内部还是之下。如果该点在高斯分布内,那么您将高斯随机数作为该点的 x 值。
这种方法并不完美,因为从技术上讲,高斯曲线趋向无穷大,并且您无法创建在 x 维度上接近无穷大的框。但是高斯曲线在 y 维度上非常快地接近 0,所以我不会担心。C 中变量大小的限制可能更多地限制了您的准确性。
2)另一种方法是使用中心极限定理,该定理指出,当添加独立随机变量时,它们会形成正态分布。牢记这个定理,您可以通过添加大量独立随机变量来近似高斯随机数。
这些方法不是最实用的,但是当您不想使用预先存在的库时,这是可以预料的。请记住,这个答案来自很少或没有微积分或统计经验的人。
Take a look at what I found.
This library uses the Ziggurat algorithm.
计算机是确定性设备。计算中没有随机性。此外,CPU 中的算术设备可以评估一些有限整数集(在有限域中执行评估)和有限实有理数集的求和。并且还进行了按位运算。数学处理更棒的集合,如 [0.0, 1.0] 具有无限点数。
您可以使用某些控制器在计算机内部收听一些电线,但它会均匀分布吗?我不知道。但是如果假设它的信号是累积值大量独立随机变量的结果,那么您将收到近似正态分布的随机变量(在概率论中得到了证明)
存在称为伪随机生成器的算法。我觉得伪随机生成器的目的是模拟随机性。好的标准是: - 经验分布收敛(在某种意义上 - 逐点,均匀,L2)到理论 - 您从随机生成器收到的值似乎是独立的。当然,从“真实的观点”来看这不是真的,但我们假设它是真的。
一种流行的方法 - 你可以对 12 irv 求和,分布均匀......但老实说,在傅里叶变换、泰勒级数的帮助下推导中心极限定理时,需要有几次 n->+inf 假设。因此,例如理论上 - 就我个人而言,我不了解人们如何以均匀分布执行 12 irv 的总和。
我在大学学过概率论。特别是对我来说,这只是一个数学问题。在大学里,我看到了以下模型:
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;
}
这种方式如何做只是一个例子,我想它存在另一种实现它的方法。
可以在 Krishchenko Alexander Petrovich ISBN 5-7038-2485-0的这本书“Moscow, BMSTU, 2004: XVI Probability Theory, Example 6.12, p.246-247”中找到它正确的证明
不幸的是,我不知道这本书是否存在英文翻译。