C语言产生标准正态分布或高斯分布随机数

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