VC学习资料收集（20）：随机数（转贴）

在论坛上看到很多人在问随机数的问题，我这里把我知道的有关随机数的一些内容介绍给大家，不正确的地方还请大家指正。
一．计算机中随机数的产生
现在，在计算机，用来产生随机数的算法是“线性同余”法。所谓线性同余，其实就是下面两个式子。假设I就是一个随机数的序列，Ij+1与Ij的关系如下：
Ij+1 =Ij * a+c  (mod m)
或是Ij+1 =Ij *a    (mod m)，
其中，不妨取a=16807，m=2147483647，以为一常数。写个简单的程序就是：
long r;

void scand( long v)//初始化随机种子数
{
r = v;
}

long rand()//产生随机数
{
r = (r*a + c)%m;//a,c,m为常数
return r;
}
再看一下稍复杂一点的：(Random () 的 Borland 的实现)
long long RandSeed = #### ;
unsigned long Random(long max)
{
long long x ;
double i ;
unsigned long final ;
x = 0xffffffff;
x += 1 ;
RandSeed *= ((long long)134775813);
RandSeed += 1 ;
RandSeed = RandSeed % x ;
i = ((double)RandSeed) / (double)0xffffffff ;
final = (long) (max * i) ;
return (unsigned long)final;
}
二．计算机产生的随机数不是真的随机数
[引:]我们建立了真正调用伪随机数生成器的 random()。但什么是伪随机数生成器？假定需要生成介于 1 和 10 之间的随机数，每一个数出现的几率都是一样的。理想情况下，应生成 0 到 1 之间的一个值，不考虑以前值，这个范围中的每一个值出现的几率都是一样的，然后再将该值乘以 10。请注意，在 0 和 1 之间有无穷多个值，而计算机不能提供这样的精度。
为了编写代码来实现类似于前面提到的算法，常见情况下，伪随机数生成器生成 0 到 N 之间的一个整数，返回的整数再除以 N。得出的数字总是处于 0 和 1 之间。对生成器随后的调用采用第一次运行产生的整数，并将它传给一个函数，以生成 0 到 N 之间的一个新整数，然后再将新整数除以 N 返回。这意味着，由任何伪随机数生成器返回的数目会受到 0 到 N 之间整数数目的限制。
在大多数的常见随机数发生器中，N 是 232? （大约等于 40 亿），对于 32 位数字来说，这是最大的值。换句话说，我们经常碰到的这类生成器能够至多生成 40 亿个可能值。而这 40 亿个数根本不算大，只是指尖这么大。
伪随机数生成器将作为“种子”的数当作初始整数传给函数。这粒种子会使这个球（生成伪随机数）一直滚下去。伪随机数生成器的结果仅仅是不可预测。由伪随机数生成器返回的每一个值完全由它返回的前一个值所决定（最终，该种子决定了一切）。如果知道用于计算任何一个值的那个整数，那么就可以算出从这个生成器返回的下一个值。
结果，伪随机数生成器是一个生成完全可预料的数列（称为流）的确定性程序。一个编写得很好的的 PRNG 可以创建一个序列，而这个序列的属性与许多真正随机数的序列的属性是一样的。例如：
PRNG 可以以相同几率在一个范围内生成任何数字。
PRNG 可以生成带任何统计分布的流。
由 PRNG 生成的数字流不具备可辨别的模式。
PRNG 所不能做的是不可预测。如果知道种子和算法，就可以很容易地推算出这个序列。

计算机产生的随机数一般都只是一个周期很长的数列，不是真的随机数。也就是说，随机数一般是伪随机数，每个随机数都是由随机种子开始的一个已定的数列（周期很长）。一般地，为了随机数更真一点，随机种子在系统中通常是参照系统时钟生成的。
看看下面这个例子,从中，读者应能有所体会。
main()
{
  int i;
  scand(time(NULL)); /＊可向计算机读取其时钟值，并把值自动设为随机数种子＊／
  for(i=0;i<10;i++){
    printf("%10d",1+(rand()%6));/*这里１是移动值，他等于所需的连续整数＊／
　｝　　　　　　　　　　　　　　／＊数值范围的第一个数；b是比例因子；他＊／
　return(0);    　　　　　　　　／＊等于所需的连续整数值的范围宽度；＊／
}
从数学上讲，为了得到一个周期性更长的随机序列，我们可以使用如下方法：（这是我在一个书上看到的，详细的情况大家可以查一查，我自己也记不清了，呵呵）
float rand(long* idum)
{
    int j;    long k;    static long iy=0;    static long iv[NTAB];
    float temp;    
    if(*idum<=0||!iy)
    {
        if(-(*idum)<1)            *idum=1;
        else            *idum=-(*idum);
        for(j=NTAB+7;j>=0;j--){
            k=(*idum)/IQ;
            *idum=IA*(*idum-k*IQ)-k*IR;
            if(*idum<0)
                *idum+=IM;
            if(j<NTAB)
                iv[j]=*idum;
        }    iy=iv[0];    }    
    k=(*idum)/IQ;
    *idum=IA*(*idum-k*IQ)-k*IR;
    if(*idum<0)        *idum+=IM;
     j=iy/NDIV;
    iy=iv[j];
    iv[j]=*idum;    
    if((temp=float(AM*iy))>RNMX)        return float(RNMX);
    else         return temp;
}

[注:]有文讲 ，根据Randomize的工作原理，利用函数Timer得到从午夜开始到现在经过的秒数，然后再根据要得到的随机数值大小对该数值进行"“衰减”处理，这样得到的数值则可称得上是真正意义的随机数值，我认为，这也是人为的方法，仍有它的确定性和周期性，仍称不上是真正的随机数，单纯改变伪随机数的生成逻辑计算方法并不能达到目的，最有效的办法只能是改变rand_seed，就是种子。而且，改变后的rand_seed不应该是人为的。（注：目前，在 Intel 的PIII处理器中内置了一个与CPU温度相关的随机数生成器，算是一个比较有效的种子生成器。）更好的办法是根据“随机事件”生成随机数，如键盘和鼠标输入值、中断、磁盘读取等等。然而，许多服务器没有键盘和鼠标，许多“黑盒”产品也不带有硬盘，因此很难找到好的随机数源，当然，通讯密钥也就一样不安全。而如网络状态等也不能很好地保证随机数的“随机性”。电器噪声和声音频率也许是很好的随机数源，但大多数人恐怕并不愿意在计算机中增加这种设备，而且也可能出现设备失灵和外部操纵（影响）等问题。对于要处理大量连接的网关服务器，是必须要考虑的问题。如果可以通过，精确检测机器cpu的通电电流强度，来作为随机数种子，或是其他一些没有人为因素的干扰的，且瞬间变化快的方法获得种子，必将能产生符和要求的随机数。

三．几种随机数的获取办法
1.产生一个范围内的随机数
一般地，我们可用j=1+(int)(n*rand()/(RAND_MAX+1.0))来生成一个0到ｎ之间的随机数。
若用int x = rand() % 100;来生成 0 到 100 之间的随机数这种方法是不可取的，比较好的做法是： j=(int)(100.0*rand()/(RAND_MAX+1.0)) ，当然，我们也可是使用random(100)。下面的例子都是用了random(n).
2、筛选型随机数 如希望取0-99的随机数，但不能是6。
解决方法：
x = random(100);
while (x==6) {
x = random(100);
}
又如希望取0-99的随机数，但不要5的倍数 解决方法：
x = random(100);
while ((x ％ 5)==0) {
x = random(100);
}
3、从连续的一段范围内取随机数。
如从40--50的范围内取随机数。 解决方法： x=random(11)+40
4、从一组乱数中取随机数。 如：从 67, 87, 34, 78, 12, 5, 9, 108, 999, 378十个数中随机取数。 解决方法：可以用数组将些十个数存贮，然后把0－－9中取出的随机数作为序号，实现随机取数。
a = new Array(67, 87, 34, 78, 12, 5, 9, 108, 999, 378);
j = random(10);
x = a[j];
四．产生具有一定分布的随机数
关于这点，有一篇文章写得很好，请看：
非均匀分布随机数的产生及其在计算机模拟研究中的应用
作者：胡性本　刘向明　方积乾
单位：胡性本（湖北职工医学院　荆州434000）； 刘向明（湖北职工医学院　荆州434000　现为华中理工大学生物工程系博士研究生）； 方积乾（湖北职工医学院　荆州434000　中山医科大学）
关键词：非均匀分布随机数；计算机模拟；程序设计
摘　要：非均匀分布随机数在进行计算机模拟研究中有重要作用，但计算机高级语言通常都只提供产生均匀分布随机数的函数，给研究工作带来困难。该文提出的方法，较好地解决了这一问题，有很强的适用性。

中图分类号：O 242.1
文章编号：1004-4337(2000)01-0059-02▲
1　问题的提出
常用的计算机高级程序设计语言，大多提供了产生在［0，1］区间内连续均匀分布的独立随机数r的函数。若将产生的随机数作简单的变换X=a+(b-a)r，就能得到在区间［a,b］上均匀分布的随机数X。如果与取整或舍入函数结合运用，还可得到离散均匀分布的随机数，给计算机模拟提供了方便。但在很多情况下进行计算机模拟需要用到按某些特定规律分布的非均匀随机数。例如，在离子通道门控模型的模拟研究中，就经常要用到服从指数分布、对数分布、正态分布、普畦松分布等非均匀分布的随机数。这时可以在源程序中编写一个函数或过程来产生。我们在近年来所进行的计算机模拟研究工作［1,2］中，大量使用了非均们分布的随机数，获得成功。由于离散随机数可以由连续随机数通过取整或舍入等途径转换而得，本文只侧重介绍非均匀分布连续随机数的产生原理和应用实例，并给出了用类PASCAL语言［3］编写的相应的函数和过程。掌握了任何一种计算机高级语言编程技术和数据结构知识的人，都能很方便地转换成能上机运行的程序，有很大的灵活性与很强的实用性。
2　产生原理
假定连续随机变量X是连续随机变量R的函数
x= (r)　　　　　　　　　　　　　　　　　　(1)
由概率统计知识可知，X与R在对应点的分布函数的数值应当相等，即产生R的反变换
F(x)=F(r)　　　　　　　　　　　　　　　　　(2)
其中，r为(1)式的解，即r= -1(x)=Ψ(x)
　　在特殊情况下，若R是在［0，1］区间内的均匀分布的连续随机变量，那么R<r的概率P(R<r)=r。此时，(2)式可以简化为
F(x)=F(r)=r　　　　　　　　　　　　　　　　(3)
其反函数
x=F-1(r)　　　　　　　　　　　　　　　　　　(4)
具有分布函数F，现予以证明：因为分布函数是单调递增函数，其反函数存在。由概率知识可知，对任意实数X，有
　　P(X<x)=P(F-1(R)<x)=P(R<F(x))
　　因为R是在［0，1］区间上的均匀分布的连续随机变量，故P(R<F(x))=F(x)，即F-1(r)具有分布函数F。
　　(4)式表明，要产生一个服从某种分布的随机数，可以先求出其分布函数的反函数的解析式，再将一个在［0，1］区间内的均匀随机数的值代入其中计算就可以了。
3　应用实例
例1　产生密度函数满足指数分布

的随机数，其中α>0。
因为　 ，故分布函数为：

由(4)式可得
　　　　　　　　　　　　　(5)
当r为［0，1］区间内的均匀随机数时，1-r也是［0，1］区间内的均匀随机数，故(5)式还可以简化为：
　　　　　　　　　　　　　　　　(6)
假定计算机高级语言已提供了产生在［0，1］区间内均匀分布的随机数的函数RND(0)，则根据(6)式可以写出产生服从指数分布的随机数的类PASCAL语言的函数如下：
　　FUNC get_exp_random(alpha:real):real;
　　｛产生服从指数分布的随机数，值参alpha为比例参数｝
　　get_exp_random=-ln(RND(0))/alpha;
ENDF； ｛get_exp_random｝
　　例2　产生服从正态分布的随机数
　　正态分布X～N(a,σ)都可以由标准的正态分布X′～N(0,1)经过变换X=a+σX′得到，在此只讨论服从标准正态分布的随机数的产生方法。标准正态分布的密度函数为：
　　　　　　　　　　　　　　　(7)
但其分布函数不能用有限的解析形式来表示，给转换带来困难，但可以用以下的变换来产生随机数。
　　假定r1与r2是［0,1］区间的两个独立的均匀分布随机数，现将其作如下的变换，令
　　　　　　　　　　　　　(8)
再将其反变换为：
　　　　　　　　　　　　　　(9)
(9)式中的C为常数。由此可导出其密度函数为：
　　　　　　　　　　　(10)
比较(10)式与(7)式，可知X1与X2是两个相互独立的服从正态分布的随机变量，因而(8)式是与(4)相当的据以产生随机数的表达式。由此可以写出产生服从标准正态分布随机数的类PASCAL语言的过程如下：
　　PROC gauss(VAR x1,x2:real);
　　｛产生服从正态分布的随机数对，用变量参数x1与x2传递随机数｝
　　pi:=3.14159265;　　｛赋π值｝
　　r1:=RND(0);　　r2:=RND(0);
　　s:=SQRT(-2*ln(r1));　　｛中间变量赋值｝
　　x1:=s*cos(2*pi*r2);
　　x2:=s*sin(2*pi*r2);
　　｛正态分布随机数赋给变量参数｝
ENDP；　｛gauss｝
　　此过程每调用一次能产生两个相互正交的标准正态分布随机数。