任意分布的随机数的产生方法—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的系统仿真技术与应用》清华大学出版社，薛定宇，陈阳泉
srand函数是随机数发生器的初始化函数，原型：
srand((unsigned) time(&t));
还有一个经常用法，不需要定义time_t型t变量,即：
srand((unsigned) time(NULL));
1、srand 选用时间作为初始化的种子，是方式每次程序重启，产生的随机数相同。
2、srand 应该是在初始化的时候设置，不必在每次调用rand之前再调用srand，如果每次都调用srand，会产生在1s内，产生的随机数相同。
3、rand() 函数在每次调用完会改变全局种子变量。这样保证每次rand产生随机数不相同。
如何产生一定范围内的随机数？
直接的方法是：
rand() % N;
返回从 0 到 N - 1 的数字。但这个方法不好，因为许多随机数发生器的低位比特并不随机。一个较好的方法是：
(int)((double)rand() / ((double)RAND_MAX + 1) * N);
如果你不希望使用 double，另一个方法是：
rand() / (RAND_MAX / N + 1);
产生0-10的一个随机数
(int) (11.0*rand()/(rand_max+1.0) )
两种方法都需要知道 RAND_MAX，而且假设 N 要远远小于 RAND_MAX。ANSI 规定标准头文件 stdlib.h 中包含 RAND_MAX 的 #define。顺便提一下，RAND_MAX 是个常数，它告诉你 C 库函数 rand() 的固定范围。你不可以设 RAND_MAX 为其它的值，也没有办法要求 rand() 返回其它范围的值。如果你用的随机数发生器返回的是 0 到 1 的浮点值，要取得范围在 0 到 N - 1 内的整数，只要将随机数乘以 N 就可以了。
RAND_MAX 的宏定义是0x7fff就是 16位有符号整形的最大正值32767，或者 0x7fffffff。
j = (int)((1.0*rand() / RAND_MAX) * (rand_max - rand_min) + rand_min) 更确切吧。因为(1.0*rand() / RAND_MAX产生的是0和1之间的随机数r，而r * (rand_max - rand_min) + rand_min产生的即是rand_min到rand_max之间的数。
