C++生成随机数：正态分布（normal distribution）

首次写博客，见谅！

win32控制台程序

 1 #include "stdafx.h"
 2 #include <stdio.h>
 3 #include <iostream>
 4 #include <time.h>
 5 #include <stdlib.h>
 6 #include <math.h>
 7 
 8 using namespace std;
 9 
10 
11 double gaussrand()
12 {
13     static double V1, V2, S;
14  　　static int phase = 0;
15  　　double X;
16     double a[100];
17     srand((unsigned)time(NULL));
18 
19 for (int i = 0; i<100; i++)
20     {
21     　　if ( phase == 0 ) 
22    　　 {
23     　　　　do 
24     　　　　{
25     　　　　　　double U1 = (double)rand() / RAND_MAX;
26　　　　　　　　　double U2 = (double)rand() / RAND_MAX;
27     　　　　　　V1 = 2 * U1 - 1;
28     　　　　　　V2 = 2 * U2 - 1;
29     　　　　　　S = V1 * V1 + V2 * V2;
30    　　　　 } 
31    　　　　 while(S >= 1 || S == 0);
32     　　　　X = V1 * sqrt(-2 * log(S) / S);
33    　　 } 
34    　　 else
35     　　　　X = V2 * sqrt(-2 * log(S) / S);
36    　　　　 phase = 1 - phase;
37 　　　　　　//X = X * 0.04 + 40;
38     　　　　a[i] = X;
39    　　　　 cout<<a[i]<<endl; 
40     }
41     return 0;
42 }
43 
44 int main(int argc, char* argv[])
45 {
47 　　gaussrand(); 
48 　　return 0;
49 }

 

这样生成的高斯分布随机数序列的期望为0.0，方差为1.0。若指定期望为E，方差为V，则只需增加：

X = X * σ + μ;

向屏幕输出 期望为0.0，方差为1.0的正态分布数据 100个

 

在此感谢博友：yeahgis

博客链接：http://www.cnblogs.com/yeahgis/archive/2012/07/13/2590485.html#2608670

 

   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);
        X = X * std + mean;

 ////////////////////////////////////////////////////////////////////////////分割线 2014/02/19

static double V1, V2, S;
    static int phase = 0;
    double X;
    srand((unsigned)time(NULL));
    for (int i = 0; i<10000; i++)
    {
        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;
        X = X * std + mean;

        if(index != 5)
        {
            b[i] = b[i] - X;
        }
        else
        {    
            b[i] = b[i] + X;
        }
        //a[i] = X;
    }

以上代码生产均值：mean，标准差：std的10000个正态分布随机数。

（注意：在进行随机数生成过程中尽量不要使用其他代码影响随机数的生成过程，一下代码不可取）

 1 static double V1, V2, S;
 2     static int phase = 0;
 3     double X1;
 4     X1 = 0.0;
 5     int i;
 6     srand((unsigned)time(NULL));
 7     //srand((unsigned int)time(&now)/(j+1));
 8     for (i = 0; i<10000; i++)
 9     {
10         if ( phase == 0 ) 
11         {
12             do 
13             {
14                 double U1 = (double)rand() / RAND_MAX;
15                 double U2 = (double)rand() / RAND_MAX;    
16                 V1 = 2 * U1 - 1;    
17                 V2 = 2 * U2 - 1;
18                 S = V1 * V1 + V2 * V2;
19             }
20             while(S >= 1 || S == 0);
21             X1 = V1 * sqrt(-2 * log(S) / S);
22         } 
23         else
24             X1 = V2 * sqrt(-2 * log(S) / S);
25         phase = 1 - phase;
26         //X1 = X1 * 0.04 + 40;
27 
28         X1 = X1 * σ + μ;
29         if(CString("减环") == m_Clist1.GetItemText(j,5))
30         {
31             b[i] = b[i] - X1;
32         }
33         else if(CString("增环") == m_Clist1.GetItemText(j,5))
34         {    
35             b[i] = b[i] + X1;
36         }
37     }

代码中的if else判断对随机数生产过程影响较大（可能是此判断过于复杂，在此存疑！）

posted @ 2013-01-27 20:35  ONWAYO  阅读(5547)  评论(0)  编辑  收藏  举报