随机数法

2024年1月31日发(作者：天一大联考二数学试卷)

用 “线性同余法”产生随机数，公式如下： x=(x*29+7) mod 1000

为什么是29和37？为什么用MOD ？

mod 1000的意思就是如果超过1000就把该数减去1000.还可以打个比方，x mod 12 代表的是钟表刻度，它的值永远在0..11之间。这个概念在数学里叫”同余”。一般：

在计算机上可以用物理方法来产生随机数，但价格昂贵，不能重复，使用不便。另一种方法是用数学递推公式产生，这样产生的序列与真正的随机数序列不同，所以称为伪随机数或伪随机序列，只要方法和参数选择合适，所产生的伪随机数就能满足均匀性和独立性，与真正的随机数具有相近的性质。

产生随机数的方法是先用一定的方法产生[0,1]均匀分布的随机数，然后通过一个适当的变换就可以得到符合某一概率模型的随机数。

常用的产生[0,1]均匀分布的随机数的方法有乘同余法和混合同余法。

用乘同余法产生[0,1]均匀分布的随机数递推公式为：

Xi=λxi-1(modM) (i=1,2,„„n) (2-13)

式中λ为乘子，M为模，当i=1时，xi=xi-1为初始参数，x0可取1或任意奇数。利用上式算出序列x1,x2„„后再取

γ i=xi / M (2-14)

即为第i个[0,1] 均匀分布的随机数。

混合同余法又称线性同余法，产生随机数的递推公式为：

Xi=(λxi-1+C)(modM) (i=1,2,„„n) (2-15)

γ i=xi / M (2-16)

根据概率论可知，[0,1]均匀分布的随机变量ζ与其它的变量ζ之间存在如下的重要关系：

η=F-1(ζ) (2-17)

式中F-1( . )为变量η的分布函数。对于F-1( . )不能表示为显函数的情况，可采用线性插值法变换，此时η与ζ的关系，可用图2-6来表示。在[0,1]均匀分布随机变量轴上，在0-1之间可划成m等份。给定分布的概率分布曲线也同时被划分成m段，将每一段曲线用折线来代替就可以得到对应于[0,1]均匀分布随机数γi的给定分布的随机数xi。

式中：γi：为[0,1]均匀分布随机数。

γj .γj+1为γi的所在区段的上、下限。

xj. xj+1为与γi对应的区段的x的上、下限。

图2-6 线性插入法变换

许多情况下，需要产生多维随机向量的样本，各个分量又是相关的，这时可以分别求出它们的边缘分布和条件分布，化成几个单变量的分布，然后依次进行抽样。

———————————————————————

一．计算机中随机数的产生

现在，在计算机，用来产生随机数的算法是”线性同余”法。所谓线性同余，其实就是下面两个式子。假设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));/*这里1是移动值，他等于所需的连续整数＊／

｝ ／＊数值范围的第一个数；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

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到n之间的随机数。

若用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= ® (1)

由概率统计知识可知，X与R在对应点的分布函数的数值应当相等，即产生R的反变换

F(x)=F® (2)

其中，r为(1)式的解，即r= -1(x)=Ψ(x)

在特殊情况下，若R是在〔0，1〕区间内的均匀分布的连续随机变量，那么R

F(x)=F®=r (3)

其反函数

x=F-1® (4)

具有分布函数F，现予以证明：因为分布函数是单调递增函数，其反函数存在。由概率知识可知，对任意实数X，有

P(X

因为R是在〔0，1〕区间上的均匀分布的连续随机变量，故P(R

(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｝

此过程每调用一次能产生两个相互正交的标准正态分布随机数。

————————————————

概率算法简介

很多算法的每一个计算步骤都是固定的，而在下面我们要讨论的概率算法，允许算法在执行的过程中随机选择下一个计算步骤。许多情况下，当算法在执行过程中面临一个选择时，随机性选择常比最优选择省时。因此概率算法可在很大程度上降低算法的复杂度。

概率算法的一个基本特征是对所求解问题的同一实例用同一概率算法求解两次可能得到完全不同的效果。这两次求解问题所需的时间甚至所得到的结果可能会有相当大的差别。一般情况下，可将概率算法大致分为四类:数值概率算法，蒙特卡罗（Monte Carlo）算法，拉斯维加斯（Las Vegas）算法和舍伍德（Sherwood）算法。

数值概率算法常用于数值问题的求解。这类算法所得到的往往是近似解。而且近似解的精度随计算时间的增加不断提高。在许多情况下，要计算出问题的精确解是不可能或没有必要的，因此用数值概率算法可得到相当满意的解。

蒙特卡罗算法用于求问题的准确解。对于许多问题来说，近似解毫无意义。例如，一个判定问题其解为”是”或”否”，二者必居其一，不存在任何近似解答。又如，我们要求一个整数的因子时所给出的解答必须是准确的，一个整数的近似因子没有任何意义。用蒙特卡罗算法能求得问题的一个解，但这个解未必是正确的。求得正确解的概率依赖于算法所用的时间。算法所用的时间越多，得到正确解的概率就越高。蒙特卡罗算法的主要缺点就在于此。一般情况下，无法有效判断得到的解是否肯定正确。

拉斯维加斯算法不会得到不正确的解，一旦用拉斯维加斯算法找到一个解，那么这个解肯定是正确的。但是有时候用拉斯维加斯算法可能找不到解。与蒙特卡罗算法类似。拉斯维加斯算法得到正确解的概率随着它用的计算时间的增加而提高。对于所求解问题的任一实例，用同一拉斯维加斯算法反复对该实例求解足够多次，可使求解失效的概率任意小。

舍伍德算法总能求得问题的一个解，且所求得的解总是正确的。当一个确定性算法在最坏情况下的计算复杂性与其在平均情况下的计算复杂性有较大差别时，可以在这个确定算法中引入随机性将它改造成一个舍伍德算法，消除或减少问题的好坏实例间的这种差别。舍伍德算法精髓不是避免算法的最坏情况行为，而是设法消除这种最坏行为与特定实例之间的关联性。

本文简要的介绍一下数值概率算法和舍伍德算法。

首先来谈谈随机数。随机数在概率算法设计中扮演着十分重要的角色。在现实计算机上无法产生真正的随机数，因此在概率算法中使用的随机数都是一定程度上随机的，即伪随机数。产生随机数最常用的方法是线性同余法。由线性同余法产生的随机序列a1,a2,„,an满足：

a0=d

an=(ban-1+c)mod m n=1,2„„.

其中,b>=0, c>=0, d>=m。d称为该随机序列的种子。

下面我们建立一个随机数类RadomNumber，该类包含一个由用户初始化的种子randSeed。给定种子之后，既可产生与之相应的随机数序列。randseed是一个无符号长整型数，既可由用户指定也可由系统时间自动产生。

const unsigned long maxshort=65536L;

const unsigned long multiplier=1194211693L;

const unsigned long adder=12345L;

class RandomNumber

{

private:

//当前种子

unsigned long randseed;

public:

//构造函数，缺省值0表示由系统自动产生种子

RandomNumber(unsigned long s=0);

//产生0-n-1之间的随机整数

unsigned short Random(unsigned long n);

//产生[0，1)之间的随机实数

double fRandom(void);

};

RandomNumber::RandomNumber(unsigned long s)

{

if(s==0)

randseed=time(0);

else

randseed=s;

}

unsigned short RandomNumber::Random(unsigned long n)

{

randseed=multiplier*randseed+adder;

return (unsigned short)((randseed>>16)%n);

}

double RandomNumber::fRandom(void)

{

return Random(maxshort)/double(maxshort);

}

函数Random在每次计算时，用线性同余式计算新的种子。它的高16位的随机性较好，将randseed右移16位得到一个0-65535之间的随机整数然后再将此随机整数映射到0-n-1范围内。

对于函数fRandom,先用Random(maxshort)产生一个0-(maxshort-1之间的整型随机序列)，将每个整型随机数除以maxshort，就得到[0，1]区间中的随机实数。

下面来看看数值概率算法的两个例子：

1.用随机投点法计算π

设有一半径为r的圆及其外切四边形，如图所示。向该正方形随机投掷n个点。设落入圆内的点在正方形上均匀分布，因而所投入点落入圆内的概率为πr^2/4r^2，所以当n足够大时，k与n之比就逼近这一概率，即π/4。由此可得使用随机投点法计算π值的数值概率算法。具体实现时，只需要在第一次象限计算即可。

double Darts(int n)

{

static RandomNumber dart;

int k=0;

for(int i=1;i<=n;i++){

double x=m();

double y=m();

if((x*x+y*y)<1)

k++;

}

return 4*k/double(n);

}

再简单举个舍伍德算法的例子。

我们在分析一个算法在平均情况下的计算复杂性时，通常假定算法的输入数据服从某一特定的概率分布。例如，在输入数据是均匀分布时，快速排序算法所需的平均时间是O(n

logn)。但是如果其输入已经基本上排好序时，所用时间就大大增加了。此时，可采用舍伍德算法消除算法所需计算时间与输入实例间的这种联系。

在这里，我们用舍伍德型选择算法随机的选择一个数组元素作为划分标准。这样既能保证算法的线性时间平均性能又避免了计算拟中位数的麻烦。非递归的舍伍德型算法可描述如下：

template

Type select(Type a[], int l, int r, int k)

{

static RandomNumber rnd;

while(true){

if(l>=r)

return a[l];

int i=l, j=l=(r-l+1);

Swap(a[i], a[j]);

j=r+1;

Type pivot=a[l];

while(true)

{

while(a[++i]

while(a[–j]>pivot);

if(i>=j)

break;

Swap(a[i], a[j]);

}

if(j-l+1==k)

return pivot;

a[l]=a[j];

a[j]=pivot;

if(j-l+1

{

k=k-j+l-1;

l=j+1;

}

else

r=j-1;

}

}

template

Type Select(Type a[], int n, int k)

{

if(k<1||k>n)

throw OutOfBounds();

return select(a, 0, n-1, k);

}

平时我们一般开始考虑的是一个有着很好平均性能的选择算法，但在最坏情况下对某些实例算法效率较低。这时候我们用概率算法，将上述算法改造成一个舍伍德型算法，使得该算法对任何实例均有效． 不过在有些情况下，所给的确定性算法无法直接改造成舍伍德型算法。这时候就可以借助随机预处理技术，不改变原有的确定性算法，仅对其输入进行随机洗牌，同样可以得到舍伍德算法的效果。还是刚才的例子，换一种方法实现：

template

void Shuffle(Type a[], int n)

{

static RandomNumber rnd;

for(int i=1;i

int j=(n-i)+i;

Swap(a[i],a[j]);

}

}