随机数也不是那么好得到的(转)

任意分布的随机数的产生方法—VC程序实现方法

下载源代码

摘要：
    随机数在实际运用中非常之多，如游戏设计，信号处理，通常我们很容易得到平均分布的随机数。但如何根据平均分布的随机数进而产生其它分布的随机数呢？本文提出了一种基于几何直观面积的方法，以正态分布随机数的产生为例讨论了任意分布的随机数的产生方法。

正文：
一、平均分布随机数的产生
    大家都知道，随机数在各个方面都有很大的作用，在vc的环境下，为我们提供了库函数rand()来产生一个随机的整数。该随机数是平均在0~RAND_MAX之间平均分布的，RAND_MAX是一个常量，在VC6.0环境下是这样定义的：

#define RAND_MAX 0x7fff

   它是一个short 型数据的最大值，如果要产生一个浮点型的随机数，可以将rand()/1000.0这样就得到一个0~32.767之间平均分布的随机浮点数。如果要使得 范围大一点，那么可以通过产生几个随机数的线性组合来实现任意范围内的平均分布的随机数。例如要产生-1000~1000之间的精度为四位小数的平均分布 的随机数可以这样来实现。先产生一个0到10000之间的随机整数。方法如下 ：

int a = rand()%10000;

然后保留四位小数产生0~1之间的随机小数：

double b = (double)a/10000.0;

然后通过线性组合就可以实现任意范围内的随机数的产生，要实现-1000~1000内的平均分布的随机数可以这样做：

double dValue = (rand()%10000)/10000.0*1000-(rand()%10000)/10000.0*1000;

则dValue就是所要的值。
   到现在为止，你或许以为一切工作都已经完成了，其实不然，仔细一看，你会发现有问题的，上面的式子化简后就变为：

double dValue = (rand()%10000)/10.0-(rand()%10000)/10.0;

   这样一来，产生的随机数范围是正确的，但是精度不正确了，变成了只有一位正确的小数的随机数了，后面三位的小数都是零，显然不是我们要求的，什么原因呢，又怎么办呢。
   先找原因，rand()产生的随机数分辨率为32767，两个就是65534，而经过求余后分辨度还要减小为10000，两个就是20000而要求的分辨率为1000*10000*2=20000000，显然远远不够。下面提供的方法可以实现正确的结果：

double a = (rand()%10000) * (rand()%1000)/10000.0;

        double b = (rand()%10000) * (rand()%1000)/10000.0;

        double dValue = a-b;

   则dValue就是所要求的结果。在下面的函数中可以实现产生一个在一个区间之内的平均分布的随机数，精度是4位小数。

double AverageRandom(double min,double max)

        {

        int minInteger = (int)(min*10000);

        int maxInteger = (int)(max*10000);

        int randInteger = rand()*rand();

        int diffInteger = maxInteger - minInteger;

        int resultInteger = randInteger % diffInteger + minInteger;

        return resultInteger/10000.0;

        }

   但是有一个值得注意的问题，随机数的产生需要有一个随机的种子，因为用计算机产生的随机数是通过递推的方法得来的，必须有一个初始值，也就是通常所说的随 机种子，如果不对随机种子进行初始化，那么计算机有一个确省的随机种子，这样每次递推的结果就完全相同了，因此需要在每次程序运行时对随机种子进行初始 化，在vc中的方法是调用srand（int）这个函数，其参数就是随机种子，但是如果给一个常量，则得到的随机序列就完全相同了，因此可以使用系统的时 间来作为随机种子，因为系统时间可以保证它的随机性。
   调用方法是srand(GetTickCount()),但是又不能在每次调用rand()的时候都用srand(GetTickCount())来初始 化，因为现在计算机运行时间比较快，当连续调用rand()时，系统的时间还没有更新，所以得到的随机种子在一段时间内是完全相同的，因此一般只在进行一 次大批随机数产生之前进行一次随机种子的初始化。下面的代码产生了400个在-1~1之间的平均分布的随机数。

double dValue[400];

        srand(GetTickCount());

        for(int i= 0;i < 400; i++)

        {

        double dValue[i] = AverageRandom(-1,1);

        }

用该方法产生的随机数运行结果如图1所示：

                  图1 400个-1~1之间平均分布的随机数
二、任意分布随机数的产生
   下面提出了一种已知概率密度函数的分布的随机数的产生方法，以典型的正态分布为例来说名任意分布的随机数的产生方法。
   如果一个随机数序列服从一维正态分布，那么它有有如下的概率密度函数：
（1-1）
其中μ，σ（ >0）为常数，它们分别为数学期望和均方差，如果读者对数学期望和均方差的概念还不大清楚，请查阅有关概率论的书。如果取μ =0，σ =0.2,则其曲线为

                  图2 正态分布的概率密度函数曲线
   从图中可以看出，在μ附近的概率密度大，远离μ的地方概率密度小，我们要产生的随机数要服从这种分布，就是要使产生的随机数在μ附近的概率要大，远离μ处 小，怎样保证这一点呢，可以采用如下的方法：在图2的大矩形中随机产生点，这些点是平均分布的，如果产生的点落在概率密度曲线的下方，则认为产生的点是符 合要求的，将它们保留，如果在概率密度曲线的上方，则认为这些点不合格，将它们去处。如果随机产生了一大批在整个矩形中均匀分布的点，那么被保留下来的点 的横坐标就服从了正态分布。可以设想，由于在μ处的f(x)的值比较大，理所当然的在μ附近的点个数要多，远离μ处的少，这从面积上就可以看出来。我们要 产生的随机数就是这里的横坐标。
   基于以上思想，我们可以用程序实现在一定范围内服从正态分布的随机数。程序如下：

double Normal(double x,double miu,double sigma) //概率密度函数

        {

        return 1.0/sqrt(2*PI*sigma) * exp(-1*(x-miu)*(x-miu)/(2*sigma*sigma));

        }

        double NormalRandom(double miu,
                  double sigma,double min,double max)//产生正态分布随机数

        {

        double x;

        double dScope;

        double y;

        do

        {

        x = AverageRandom(min,max); 

        y = Normal(dResult, miu, sigma);

        dScope = AverageRandom(0, Normal(miu,miu,sigma));

        }while( dScope > y);

        return x;

        }

参数说明：double miu：μ，正态函数的数学期望
                  double sigma：σ，正态函数的均方差
                  double min,double max，表明产生的随机数的范围
用如上方法，取 μ=0，σ=0.2，范围是-1~1产生400个正态随机数如图3所示：

                   图3 μ=0， σ=0.2，范围在-1~1时的400个正态分布的随机数分布图
取 μ=0， σ=0.05，范围是-1~1产生400个正态随机数如图4所示：

                   图4 μ=0，σ=0.05，范围在-1~1时的400个正态分布的随机数分布图
从图3和图4的比较可以看出， 越小，产生的随机数靠近 的数量越多，也说明了产生的随机数靠近 的概率越大。
我们，先产生4000个在0到4之间的正态分布的随机数，取μ=0，σ=0.2，再把产生的数据的数量做个统计，画成曲线，如下图5所示：

                   图5 μ=0， σ=0.2，范围在0~4时的4000个正态分布的随机数统计图
从图5中也可以看出，在靠近 处的产生的个数多，远离 处的产生的数量少，该图的轮廓线和概率密度曲线的形状刚好吻合。也就验证了该方法的正确性。
有了以上基础，也就用同样的方法，只要知道概率密度函数，也就不难产生任意分布的随机数，方法都是先产生一个点，然后进行取舍，落在概率密度曲线下方的点就满足要求，取其横坐标就是所要获取的随机数。

参考文献:
1、《概率论与数理统计》高等教育出版社，盛骤，谢式千等
2、《基于Matlab/Simulink的系统仿真技术与应用》清华大学出版社，薛定宇，陈阳泉

 

 

 

最新评论 [发表评论] [文章投稿]

查看所有评论 推荐给好友 打印

完全错误的算法，如果rand()是生成的是0到32767上平均分布的随机数，那么，rand()*rand()完全不是平均分布。何况rand()只是用线性同余法“模拟”随机数。

你这个算法只能说是生成4位小数精度的随机数，不过离平均分布差的远 ( bloodywind 发表于 2009-4-5 19:00:00)
 
楼主的本帖相信对我们都很有用。但时间紧有一处不是很明白。
    如楼主所说："实现任意范围内的平均分布的随机数。例如要产生-1000~1000之间的精度为四位小数的平均分布的随机数可以这样来实现。先产生一个0到10000之间的随机整数。方法如下 ： 
int a = rand()%10000;"

    我觉得这样得到的a不是在0到10000之间的平均分布的随机数，而用a作为平均分布的随机数的输入是有问题的。
    既然rand()得到一个0~32767之间平均分布的数，那么很容易证明：a的值属于0~2767与2768~9999的频率的比值应该是4：3。即P(0~2767):P(2768~9999)=4:3。
    因此int a = rand()%10000;我认为应改成:
//<--------------------
    int a, b;
    do{
        b = rand();
    }while(b > 29999);//截去大于29999的值。
    a = b%10000;
//-------------------->
tasbbon@gmail.com
http://gkbc.bokee.com
    回复我啊。
( tasbbon 发表于 2006-12-3 15:12:00)
 
y = Normal(dResult, miu, sigma);
中的dResult确实应该是x，这是作者的疏忽，多谢各位指出。关于公式应该是
return 1.0/(sqrt(2*PI)*sigma) 
* exp(-1*(x-miu)*(x-miu)/(2*sigma*sigma));
因此公式（1－1）也是不正确的。
( bsrong 发表于 2006-8-10 19:09:00)
 
首先非常感谢作者能给我们一个这么宝贵的算法，，谢谢。我对dresult含义也不是很清楚，但我后来用vb照着编了一个正态随机的算法，我个人认为这个 dresult应该是x，不知是否正确。。谢谢。。 ( xz403 发表于 2005-4-20 20:57:00)
 
请问文章中的dresult指的是什么?
( seeyour 发表于 2005-4-18 16:57:00)
 
原文的公式确实写错了，你的是正确的 ( pipipipi 发表于 2005-4-15 23:45:00)
 
请问你的程序中的一维正态概率密度函数公式是不是写错了，应该是return 1.0/    (sqrt(2*PI)*sigma) 
* exp(-1*(x-miu)*(x-miu)/(2*sigma*sigma));
还是你故意写成那样的啊，小弟不明特请指教 ( songzhiguan1981 发表于 2005-4-12 13:33:00)
 ===================================================