C语言 产生标准正态分布或高斯分布 随机数
产生正态分布或高斯分布的三种方法:
1. 运用中心极限定理(大数定理)
1 #include
2 #include
3
4 #define NSUM 25
5
6 double gaussrand()
7 {
8 double x = 0;
9 int i;
10 for(i = 0; i < NSUM; i++)
11 {
12 x += (double)rand() / RAND_MAX;
13 }
14
15 x -= NSUM / 2.0;
16 x /= sqrt(NSUM / 12.0);
17
18 return x;
19 }
2.利用有box 和 muller 提供的,在 knuth的网上讨论过的方法 (比较常用的方法)Box-Muller,一般是要得到服从正态分布的随机数,基本思想: 先得到服从均匀分布的随机数; 然后再将服从均匀分布的随机数转变为服从正态分布.Box-Muller 是产生随机数的一种方法。Box-Muller 算法隐含的原理非常深奥,但结果却是相当简单。如果在 (0,1] 值域内有两个一致的随机数字 U1 和 U2,
可以使用以下两个等式中的任一个算出一个正态分布的随机数字 Z:
Z = R * cos( θ ) 或 Z = R * sin( θ )
其中, R = sqrt(-2 * ln(U2)), θ = 2 * π * U1
正态值 Z 有一个等于 0 的平均值和一个等于 1 的标准偏差,可使用以下等式将 Z 映射到一个平均值为 m、标准偏差为 sd 的统计量 X:
X = m + (Z * sd)
C代码: (计算机编程中, log函数==ln()函数,以e为底的自然对数, log10 才是以10为底的函数)
1 #include <stdlib.h>
2 #include <stdio.h>
3 #define PI 3.141592654double double gaussrand( )
4 {
5 static double U, V;
6 static int phase = 0;
7 double z;
8
9 if(phase == 0)
10 {
11 U = rand() / (RAND_MAX + 1.0);
12 V = rand() / (RAND_MAX + 1.0);
13 Z = sqrt(-2.0 * log(U))* sin(2.0 * PI * V);
14 }
15 else
16 {
17 Z = sqrt(-2.0 * log(U)) * cos(2.0 * PI * V);
18 }
19
20 phase = 1 - phase;
21 retrn Z;
22 }
C++代码:
1 #include <cstdlib>
2 #include <cmath>
3 #include <limits>
4 double generateGaussianNoise(double mu, double sigma)
5 {
6 const double epsilon = std::numeric_limits<double>::min();
7 const double two_pi = 2.0*3.14159265358979323846;
8
9 static double z0, z1;
10 static bool generate;
11 generate = !generate;
12
13 if (!generate)
14 return z1 * sigma + mu;
15
16 double u1, u2;
17 do
18 {
19 u1 = rand() * (1.0 / RAND_MAX);
20 u2 = rand() * (1.0 / RAND_MAX);
21 }
22 while ( u1 <= epsilon );
23
24 z0 = sqrt(-2.0 * log(u1)) * cos(two_pi * u2);
25 z1 = sqrt(-2.0 * log(u1)) * sin(two_pi * u2);
26 return z0 * sigma + mu;
27 }
3 使用最初有marsaglia 提供的方法
1 #include <stdlib.h>
2 #include <stdio.h>
3 double gaussrand()
4 {
5 static double V1, V2, S;
6 static int phase = 0;
7 double X;
8
9 if(phase == 0)
10 {
11 do{
12 double U1 = (double)rand() / RAND_MAX;
13 double U2 = (double)rand() / RAND_MAX;
14
15 V1 = 2 * U1 - 1;
16 v2 = 2 * U2 - 1;
17 S = V1 * V1 + V2 * V2;
18 }while( S >= 1 || S ==0)
19
20 X = V1 * sqrt (-2 * log(S) / S);
21 }
22 else
23 {
24 X = V2 * sqrt(-2 * log(S) / S);
25 }
26
27 phase = 1 - phase;
28 return X;
29 }
参考: http://blog.chinaunix.net/uid-22666248-id-357093.html
https://en.wikipedia.org/wiki/Box%E2%80%93Muller_transform
原文链接: https://www.cnblogs.com/tsingke/p/6194737.html
欢迎关注
微信关注下方公众号,第一时间获取干货硬货;公众号内回复【pdf】免费获取数百本计算机经典书籍
原创文章受到原创版权保护。转载请注明出处:https://www.ccppcoding.com/archives/246035
非原创文章文中已经注明原地址,如有侵权,联系删除
关注公众号【高性能架构探索】,第一时间获取最新文章
转载文章受原作者版权保护。转载请注明原作者出处!