C語言產生標准正態分布或高斯分布隨機數

C語言 產生標准正態分布或高斯分布 隨機數

 

產生正態分布或高斯分布的三種方法: 


1. 運用中心極限定理(大數定理)

#include 
#include 

#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.141592654double 
　　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); 18     }

    phase = 1 - phase;
    retrn 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)
     
          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