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

产生正态分布或高斯分布的三种方法:

1. 运用中心极限定理(大数定理)

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#define NSUM 25

 double gaussrand() {
      double x = 0;
      int i; 
      for(i = 0; i < NSUM; i++) 
      {         
            x += (double)rand() / RAND_MAX;
      }  
      x -= NSUM / 2.0; 
      x /= sqrt(NSUM / 12.0);  
     return x;
 }

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为底的函数)

#include <stdlib.h>
#include <stdio.h>
#define PI 3.141592654
double gaussrand( )
{
     static double U, V; 
     static int phase = 0;
     double z;  
     if(phase == 0) 
     {           
          U = rand() / (RAND_MAX + 1.0); 
          V = rand() / (RAND_MAX + 1.0);           
          Z = sqrt(-2.0 * log(U))* sin(2.0 * PI * V);  
      } else {         
          Z = sqrt(-2.0 * log(U)) * cos(2.0 * PI * V); 
      } 
     phase = 1 - phase; 
     return Z; 
}

C++代码:

#include <cstdlib>
#include <cmath>
#include <limits>
double generateGaussianNoise(double mu, double sigma)  {
      const double epsilon = std::numeric_limits<double>::min();
      const double two_pi = 2.0*3.14159265358979323846;

      static double z0, z1; 
      static bool generate;
      generate = !generate; 
      if (!generate)
           return z1 * sigma + mu; 
      double u1, u2; 
      do { 
           u1 = rand() * (1.0 / RAND_MAX); 
           u2 = rand() * (1.0 / RAND_MAX); 
      } while ( u1 <= epsilon ); 
      z0 = sqrt(-2.0 * log(u1)) * cos(two_pi * u2);
      z1 = sqrt(-2.0 * log(u1)) * sin(two_pi * u2);
      return z0 * sigma + mu; 
}

3 使用最初有marsaglia 提供的方法

#include <stdlib.h>
#include <stdio.h>
double gaussrand()  {
       static double V1, V2, S; 
       static int phase = 0;
       double X; 
       if(phase == 0) { 
             do{ 
                    double U1 = (double)rand() / RAND_MAX; 
                    double U2 = (double)rand() / RAND_MAX;                 
                    V1 = 2 * U1 - 1;
                    V2 = 2 * U2 - 1; 
                    S  = V1 * V1 + V2 * V2; 
             }while( S >= 1 || S ==0) 19      
             X = V1 * sqrt (-2 * log(S) / S); 
       }else{            
             X = V2 * sqrt(-2 * log(S) / S); 
       } 
       phase = 1 - phase; 
       return X; 
}

参考: http://blog.chinaunix.net/uid-22666248-id-357093.html
https://en.wikipedia.org/wiki/Box%E2%80%93Muller_transform