2012数据结构提交第十一章随机数

第十一章 随机数许多应用中都要用到随机数，例如密码学、仿真、抽样、数值分析、计算机程序设计、决策、数据挖掘、进化计算、美学、娱乐等。本章主要介绍一些常用的随机数产生方法、检验方法及随机数的应用。主要内容生成随机数随机数检验随机排列与随机组合应用生成随机数随机数被应用在许多不同的研究领域中程序设计：比较解决同一问题的不同算法的性能例如对于排序算法的测试，在输入规模为1000的情况下执行10000次排序，除一些特定的输入数据（已排序或逆序的数据）外，采用一个随机数生成程序来产生这些测试数据是一个很好的选择随机算法中：随机算法对于相同的输入会执行不同的运行过程从而得到不同的输出，可以用随机数来确定算法的下一步运行过程。例如在棋类（围棋、五子棋等）程序中，计算机可以在不同的开局中选择一种，对手下棋后，计算机要在一系列可能的走法中选择一种，这些都可以用随机数来进行选择。生成随机数为追求真正的随机序列，人们曾采用很多种原始的物理方法用于生成一定范围内满足精度（位数）的均匀分布数值序列。可以用放射性同位素的衰变来产生随机数.1942年，RAND公司以随机脉冲源，用电子旋转轮产生随机数表；1975

至1978年，Frigerio做了一个高分率的计数器，每小时能产生6000个31位的随机数。但使用物理方法产生随机数的缺点在于：速度慢、效率低、需占用大量存储空间、设备的维修费用昂贵且不可重现等。生成随机数在计算机问世后不久，人们转而探求用计算机产生随机数的有效方法：利用数学方法，根据特定的迭代公式在计算机上产生随机数。对此J.V.Neumann曾说过：任何考虑用算术方法产生随机数字的人都是在做错事。因为任何数据都要靠产生它们的算法来得到，不可能随机。一般来说，由计算机程序产生的随机序列通常被称作伪随机序列，是由递推公式和初值推出的。但在一定范围的应用中，只要伪随机数能通过一系列统计检验，就可以把它们当成“真”随机数使用。生成随机数1946年，J.V.Neumann提出了第一个利用计算机的普通算术操作产生随机数的方法，被称为“平方取中法”取一个N位数做为初始值，从其平方（2N位的数，如果不足可以在左边补０）中抽取中间的N位作为序列中的下一个数，依次类推产生随机序列。例如，如果N=10，以为初始值，把它平方后得到29，于是下一个数就是

3413663732，接下来的几个数分别为1000751721，5040070844，314112538等等。尽管平方取中法描述起来很简单，但它在实际应用中却并不是一个产生随机数的好方法，原因在于产生的序列容易出现重复元素的短循环；现已很少被采用。11.1

生成随机数均匀分布随机数–

11.1.1.1

线性同余法–

11.1.1.2

其它方法其它分布随机数–

11.1.2.1

正态分布–

11.1.2.2

指数分布–

11.1.2.3

泊松分布均匀分布随机数均匀分布的随机序列是指在一个指定范围内，所有数出现的概率相等。最重要、最基本的随机序列是在[0,1]上均匀分布的随机序列。一般地，其他分布的随机序列可以借助于均匀分布的随机序列来实现。因此，我们首先考虑生成均匀分布随机序列的方法：由这些方法得到的随机数处于0和某个整数m-1之间，然后通过简单的计算，就可以变换为[0,1]上均匀分布的随机数。均匀分布随机数线性同余法–1951年，D.H.Lehmer首次提出用线性同余方法生成随机数，使得随机数生成方法的研究向前迈进了一大步。线性同余方法是由一个初始值X0，三个数m，a，c及以下递推关系式产生序列X1，X2，…，这个序列被称为线性同余序列Xi+1

=

(aXi+c)

mod

m模数m，乘数a，增量c和初始值X0为整数且分别满足m>0,0≤a<m，0≤c<m，0≤X0<m，0≤i<m.由此递推式产生的随机整数介于0~m-1之间，初始值X0称为种子。若想得到0和1之间的随机序列，只需令Un

=Xn/m即可。均匀分布随机数线性同余法–

例：m

=10，a

=c

=5，X0=1，得到的序列为：10，0，5，8，1，10，0，5，…在序列中第二次产生同一个数会导致重复的序列，上例说明了同余序列会进入一个循环的事实。序列中没有重复的最大长度称为该序列的周期。一个好的随机序列应当有相当长的周期，周期长度为m-1显然是最好的。均匀分布随机数线性同余法增量c的取值会影响生成随机数的时间和序列周期长度：c

=

0时随机数的生成过程比

c

„

0时要稍快些，但c

=0会缩短序列的周期长度。Thomson和Rotenberg分别在1958年和1960年指出：取c

„0会得到更长周期的随机序列。通常分别称c

=0和c„0的线性同余法为乘同余法和混合同余法。均匀分布随机数线性同余法–

模数的选择由于序列的周期不可能大于m，为了使序列的周期尽可能的大，我们希望m的取值也稍微大些对m取值的另一个原则是计算速度快。对于32位的二进制计算机，令w

=232.若取m=w，计算上很方便，但序列Xi的低位与高位相比，其随机性弱得多；而当m=w-1时，情况就不同了，序列Xi的低位和高位一样地随机，但此时计算时间上比m=w时差些。还有一种方案是令m为小于w的最大素数，例如取m

=231-1，是一个被研究者推荐使用的取值。均匀分布随机数线性同余法模数的选择乘数的选择合理地选择乘数a，能够得到具有长周期的序列。

研究表明：当m是不同素数的积时，只有a

= 1才产生全周期，但当m可被某个素数的高次幂整除时，a有相当大的选择范围。定理11.1

由m，a，c和X0所定义的线性同余序列有周期长度m当且仅当c与m互素；对于m的每个素因子p，a-1是p的倍数；如果m是4的倍数，则a-1也是4的倍数。引理11.1设p是一个素数，l是一个正整数，其中pl>2，如果x

≡

1

(modulo

pl)，

x

≢

1

(modulo

pl+1

)，则xp≡1

(modulo

pl+1

)，

xp

≢1

(modulo

pl+2

).证明：对于不是p

的倍数的某个整数q，有x

=1+qpl.由二项式公式有3pllq

ppp-1

lp-1l+12

2l

p-1

(

)l1

p

-1

x

=1+p-1

(

)

+qp

ppl

p

p

21

1

1

qp

+=1+qp

1+p

pp

p

p

p

qp

++q

p

++

q

p圆括弧中的量是一个整数，而且圆括弧中除去第一项之外，各项的系数都可被p

整除，因此各项都是p

的倍数，因此对于1<k<p，1

k

qk-1

p(k-1)lp

p

可为p(k-1)l

所整除。括弧中最后一项为qp-1p(p-1)l-1，当pl>2

时，(p-1)l>1，因此它可为p

所整除，ppq

pp-1l-12

2l-1p-1

(

)l-1

2

3

x

=1+qpl+1+qpl+1qp

+

p

p

p

q

p

++因此，xp

≡1+

qpl+1

(modulo

pl+2)，证毕。▐证明：对t

用归纳法，只须证明，如果m1

与m2

互素，则由(X0，a，c，m1m2)所确定的线性同余序列的周期长度l是序列(X0

mod

m1，a

mod

m1，c

mod

m1，m1)和(X0

mod

m2，a

mod

m2，c

mod

m2，m2)的周期长度l1

和l2

的最小公倍数。这三个序列的元素分别记为Xn，Yn，Zn，则有Yn

=Xn

mod

m1，Zn

=Xn

mod

m2，由于m1

与m2

互素，则Xn

=

Xk

当且仅当

Yn

=

Yk

和

Zn

=

Zk

(11-3)设l’是l1

和l2

的最小公倍数，我们希望证明l’=l

.由于对于所有适当大的n，Xn

=

Xn+λ，有Yn

=Yn+λ（因此l是l1

的倍数），且Zn

=Zn+λ（因此l是l2

的倍数），所以必有l

‡l’.其次，我们知道，对于所有适当大的n，Yn

=Yn+λ’，Zn

=Zn+λ’；因此由式(11-3)，Xn

=Xn+λ’这说明l

£

l’；所以l

=l’.▐1

tll1

t=

p

p0，则由(X

，a，c，m)所确定的线性同余序0j

jj

jlllljjjjjmod

p

,a

mod

p

,c

mod

p

,p

)的周期长度l

的最小引理11.2

设把m

分解成素因子m列的周期长度l，是诸线性同余序列(X公倍数，其中，1≤j≤t.引理11.3

假定1<a<pl，其中p

为素数，如果l是使得(aλ-1)/(a-1)”0(modulopl)的最小正整数，则a

”

1(mod

ulo

p),l

=

pla

”

1(mod

ulo

4),当且仅当当p

>2当p

=2证明：假设l＝pl，如果a

≢1(modulo

p)，则当且仅当an-1

”0(modulo

pl)时，(aλ-1)/(a-l

l于是

，

条

件

a

p

-1

”

0

(modulo

pl

)

意

味

着

a

p

”

1

(modulo

p)

；

但1)

”

0

(modulo

pl).la

p”a

(mod

ulo

p)，因此a≢1(modulo

p)导致矛盾。如果p

=2

和a

”3(modulo

4)，则有(a

2l

-1

-1(a

-1)”

0

(

mod

ulo

2l

)这些论证说明，每当l

=pl

时，一般地必有a

=1+qpf，其中pf>2

且q

不是p

的倍数。剩下要证明的是，这个条件对于导致l

=pl

是充分的。通过重复地应用引理11.1，可以发现对所有g

‡0，有”

1

(mod

ulo

p

f

+

g

),ga

p

≢1

(mod

ulo

p

f

+

g

+1

)ga

p因此，(a

pg-1

(a

-1)”

0

(mod

ulo

pg

)

,(a

pg-1

(a

-1)≢

0

(mod

ulo

pg

+1

)(11-5)特别是，(a

p

-1)(a

-1)”0(mod

ulo

pl

).现在，同余序列(0，a，1，pl)有X

=(an-1)/(a-1)lnmod

pl；因此它有长度为l

的周期，即当且仅当n

是l

的倍数时Xn

=0.

因此pl

是l的倍数仅当对于某个g，l

=pg

时才可能出现，而且式(11-5)意味着l

=pg，从而完成证明。▐经过研究者多年的深入研究，对于大多数的应用而言，31位的素数m

=231-1是常用的选择.1988年，Park和Miller提出：c

=0时，取m

=231-1和a=16807，可以得到一个最大的周期，它就是此类随机数生成器中最著名的最小标准随机数生成器。它经受了大量的统计测试，并且绝大多数都成功通过。线性同余法的ADL描述：算法LR(a，c，m，x.y)/*线性同余法*/LR1.

[乘法]i‹

a*x+c.LR2.

[模]y‹

i

MOD

m.

▐均匀分布随机数线性同余法其它方法–线性同余法是被研究得最透彻的方法，但它并不是生成随机数的唯一的方法。例11.2：在Xi+1

=(aXi+c)mod

m

中引入常数，使之成为：Xi+1

=

(48271Xi+1)

mod

(231-1)但此时如果取X0

=179424105，X1=

(48271*179424105+1)

mod

(231-1) =

179424105，得到的序列周期长度仅为1.例11.3：对Xi+1=(aXi+c)mod

m(1)做如下修改：Xi+1

=

((aXi)

mod

(m+1)+c)

mod

m

(2)式(2)产生的序列是否比式(1)产生的序列更随机呢？可以看到，式(2)的形式进入了Xi+1=f(Xi)的类型，其中f是随机选择的，本章开头提到的平方取中法就属于这种类型。它已经完全不符合线性同余法的理论基础，产生的序列很容易进入短循环。因此，若要对一个好的随机数生成器进行修改，必须保证其理论基础不被破坏，才可能得到“更好”的序列。二次同余法对线性同余法的另一种改进是二次同余方法，它实现了对线性同余法的真正改进：(2i

ii+1X

=

aX

+dX

+c

modm(11-8)如果对a，d，c

的取值加入限制条件，由(11-8)能够得到全周期序列。定理11.2

二次同余序列(11-8)有长度m的周期，当且仅当下列条件成立：c与m互素；对于所有m的奇素因子p，d和a-1都是p的倍数；如果m是4的倍数，则d是偶数，而且d

”a-1

(modulo

4)；如果m是2的倍数，则d

”a-1(modulo2)；如果m是9的倍数，则d

≢3c

(modulo

9).加法生成器有一种改进的方法是使得Xi+1依赖于Xi和Xi-1，而不仅仅依赖于Xi，这样，序列的周期长度最长可以达到m2.这是因为序列在出现(Xi+λ，Xi+λ+1)=(Xi，Xi+1)之前不会重复（l为整数）这类似于斐波那契序列，但简单的斐波那契序列的结果却并不令人满意.1958年，G.J.Mitchell和D.P.Moore基于该思想提出一个较好的加法生成程序，使Xi依赖于Xi-j和Xi-k，0<j,k≤i.加法生成程序由式(11-9)产生随机序列：Xi

=

(Xi-24+Xi-55)

mod

m

,

i≥55

(11-9)算法AR(X，Y，i，j，k)执行前，用一个随机数生成器产生55

个随机数X1,X2

,…,X55并将其存储在循环列表Y

中.Y[1]，Y[2]，…，Y[55]的值分别置为X55，X54，…，X1，j

和k的初始值分别为24

和55，t

是使2t≤w

的整数，逐次执行本算法产生X56，X57，…作为输出算法AR(X，Y，i，j，k，t)/*加法生成器*/AR1.[加法]/*若此时要输出Xi，Y[j]现在等于Xi-24，Y[k]等于Xi-55.*/Y[k]‹(Y[k]+Y[j])mod

2t.PRINT

(Y[k]).AR2.[前进]/*修改j

和k

的值，以便计算下一个Xi.*/j‹

j-1.

k‹

k-1.IF

j=0

THEN

j‹

55.ELSE

IF

k=0

THEN

k‹

55.

▐随机序列的组合方法前面介绍了二次同余法和加法方法，但研究者发现它们仍然不能满足一些实际应用的需求，所生成的序列不是充分随机的。具体思想是：对于由两种无关方法生成的两个随机序列<Xi>和<Yi>，先将X中的前k个数赋值到序列V中V[1],V[2],…V[k]，k为某个适当选择的数，通常在100左右，令序列Y中的数经过简单数学变换生成新数j，将序列X中的数X[k+1]赋值给V[j]，上述过程重复进行，新序列V即是最终得到的新随机序列。算法CR

就是基于这一思想提出的，逐个输出“更随机”的数。算法开始前，以<Xi>序列的头k

个值来初始化V.算法CR(X，Y，i，j，k.V)/*两个随机序列的随机组合算法。*//*基于两个随机序列，通过洗牌实现随机化。*/CR1.[生成X，Y]/*置p

和q

分别等于序列<Xi>和<Yi>的下一项。*/p‹

X[i].

q‹

Y[i].CR2.[产生j]/*m

是产生序列<Yi>的模，因此j

是一随机值，由Y

确定，0≤j<k.*/j‹

k*q/m.CR3.[交换]PRINT

(V[j]).V[j]‹

p.

▐上述方法也被称为洗牌方法。存在一个更好的洗牌方法，它类似于算法CR，但只要求输入一个随机序列，却可以给出具有更好性能的随机序列。首先用某种给定的方法生成随机序列<Xi>，使用辅助数组V[0],V[1],…,V[k-1]逐个输出“更随机的”序列。算法BR

开始前，用<Xi>序列的头k

个值来初始化V，y

置成<Xi>序列的第

k+1

个值，y

经过变换后得到随机变量j，输出V[j]为新的随机数，然后将V[j]更新为X[k+1].算法BR(X，j，k，y.V)/*基于一个随机序列的随机组合算法。*//*基于一个随机序列，通过洗牌实现随机化。*/BR1.[抽取j]/*

m

是生成序列<Xi>的模，j

是由y

确定的一个随机值，0≤j<k.*/j‹

k*y/m.BR2.[交换]y‹

V[j].k‹

k+1.PRINT(y).V[j]‹

X[k].

▐两种组合方法的比较即使把算法CR应用于象斐波那契那样的非随机序列，也能产生出令人满意的随机序列来。但如果<Xi>和<Yi>是密切相关的话，算法CR也许产生出不如原来序列那么随机的序列。但算法BR却不会出现这样的问题，因为算法BR不使任何序列减少随机性，而且由于它可能以非常小的额外代价增强随机性，因此推荐它和任何其它的随机数生成程序组合在一起使用。然而洗牌方法有一个固有的缺点：它们只改变所生成数的顺序，而不是数本身。其它分布随机数并不是所有的应用都需要均匀分布的随机数，某些仿真过程需要其它类型的随机数，例如在研究排队和工序问题时，独立事件（顾客或工件的到达）发生的间隔往往服从指数分布，需要用到指数分布的随机数。从原理上说，其它类型的随机数都可以由均匀分布的随机数来产生。其它分布随机数正态分布除均匀分布以外，一个最重要的连续分布是正态分布或称高斯分布，即所有数值点连续地分布在其平均值左右，而且离平均值远的数据与离平均值近的数据相比其出现概率要小得多，个分布在统计学里有着极其重要的地位常用的正态分布是具有均值为0和标准方差为1的正态分布，成为标准正态分布，记为N(0,1)

分布，其分布函数为：21x2dt2pe-t-¥F

(x)=其它分布随机数正态分布–1958年，Box和Muller给出了由两个均匀分布的随机变量生成两个正态分布的随机变量的算法，称为配极法设u1,u2

是区间(0,1)上均匀分布的随机变量（可以由随机整数X

除以模数m

得到）且相互独立。令x1

=

(-2log(u1))*cos(2pu2

)，x2

=

(-2log(u1))*sin(2pu2

)，那么x1,x2

服从N(0,1)

分布，且相互独立。其它分布随机数正态分布–

1964

年，Marsaglia和Bray提出了一种改进算法，避免使用三角函数。算法中的UR(a，c，m，x.i)表示生成均匀分布随机数的算法，可以是线性同余法、加法生成器或其它的算法，在本章中始终有这样的含义。算法NR(.nx，ny)/*正态分布随机数生成算法*//*计算两个独立的正态分布变量nx

和ny.*/NR1.[得到两个均匀分布的随机数]/*

v1，v2

是-1

和+1

之间的实数。*/UR(a，c，m，x.i).UR(a，c，m，y.j).i‹

i/m.

j‹

j/m.v‹

2*j-1.u‹

2*i-1.NR2.[计算s]s‹

u2+v2.NR3.

[s≥1?]IF

s≥1

THEN

GOTO

NR1.NR4.[计算nx,ny]IF

s=0

THEN

nx‹

ny‹

0.ELSE

（

nx

‹

u

(-2

ln

s

)

/

s

.

ny

‹

v

(-2

ln

s

)

/

s

.）▐由此算法得到的nx，ny

服从N(0,1)分布，且相互独立。其它分布随机数正态分布指数分布–指数分布的分布函数如下：–采用对数方法生成指数分布的随机数相对较容易：如果u是一个在区间(0,1)上均匀分布的随机数，则令x

=

-t*ln(u)x就是以t为平均值的服从指数分布的随机数F

(x

=1-

e-x

mx

‡

0其它分布随机数正态分布指数分布泊松分布–泊松分布（Poissondistribution）刻画了小概论事件发生的次数，如彩票中奖的事件等。泊松分布与指数分布有关，因此我们可以通过指数分布或均匀分布的随机数得到泊松分布的随机数。其它分布随机数正态分布指数分布泊松分布–首先生成具有均值1/a的独立指数分布的随机数u1,u2…，一旦u1+…+ur≥1则停止，令p=r-1，即得到一个服从泊松分布的随机数p.11.2

随机数检验序列的随机性和周期长度是衡量随机序列的两个重要标准，这一节我们来介绍一些检验序列随机性的方法。对于任何随机序列来说，令EX1，EX2…，是一些有效的检验，即使序列通过了检验EX1，…，EXn，也不能保证序列一定能够通过下一个检验EXn+1；每通过一种检验只能更增强我们对该序列的信心。一般来讲，如果对一个序列应用五六个不同类型的检验，并且该序列都通过了这些检验，则认为它是随机的。对伪随机数序列使用哪种检验方法，不仅取决于产生伪随机数的方法，而且取决于伪随机数的应用需求。存在两种类型的检验：经验检验是让计算机来处理序列中出现的数据，并计算一些统计量，以此来评价序列的随机性；理论检验使用形成该序列的递推关系式的数论方法来确定序列的特性。随机数检验一般检验方法–

c2检验c2分布是统计学中常用的抽样分布.c2检验也许是最著名、最基本的随机数检验方法，可以与其它检验方法一起使用。我们先来看一个例子，假定把两个硬币投掷100

次，每次投掷的结果和每个事件的概率及其期望的投掷数如下所示：两个正面一正一反两个反面观察结果Y285121事件概率ps1/41/21/4期望值

nps255025从上述结果中我们可以说投掷这两枚硬币时其结果是偏向于正面，但是这种偏向明显吗？计算每个事件的观察结果与期望值的方差和，是否能度量硬币的偏倚情况？对于上例而言，相对于“两个正面”的情况，在“一正一反”

情况下得到的方差应该更大，因为“一正一反”更有可能发生；或者说“两个正面”所对应的差距3比起“一正一反”所对应

的差距1意义更重要。因此应当将方差与事件的概率结合起来判断结果的偏倚性：V=

(28-25)2/25+(51-50)2/50+(21-25)2/25

=

1.02这就是在投硬币的实验中，三个观察量的c2统计。c2方法的一般应用假设我们对一个随机事件做n次独立的观察，即一次观察的结果对其它观察的结果没有影响每个观察得到的结果落入k个范畴（k种情况）之一设ps是每次观察落入范畴s的概率，1≤s≤k，并设Ys是真正落入范畴s的观察数，于是有统计–

由

Y1+

Y2+…+

Yk

=

n和

p1+p2+…+pk

=1，得(

)2ks

snpsY

-

npV

=

s

=11k

Y

2

V

=n

s

=1

ps

s

-

n对于通过观察而计算出的V值，如何检验它的合理性？我们可以在一个c2分布的百分率表中查看这些值。对于不同的ps和v，c2分布的值已经制成表格（几乎所有的数理统计书中都有c2分布的附表），自由度v=k-1，如果表中第v行第p列的值是x，则表明“如果n充分大，下表中的量V将以大约p的概率小于或等于x”.p=1%p=5%p=25%p=50%p=75%p=95%p=99%v=10.000160.003930.10150.45491.3233.8416.635v=20.020100.10260.57541.3862.7735.9919.210v=30.11480.351812132.3664.1087.81511.34v=40.29710.71071.9233.3575.3859.48813.28v=50.55431.14552.6754.3516.62611.0715.09在投硬币的例子中，有三个范畴，因此自由度v

=2看表中v=2这一行，其中p=1%这一列对应的值（0.02010）表示“对于所计算的值V以大约1%的概率小于或等于0.

0201”换句话说，如果重复实验100次，其中只有1次取V<0.0201，如果在多次试验中发现计算出的V均小于0.0201，那么就证实投硬币不是随机事件再来看上例的结果，V=1.02，落在25%和50%之间，它并不是过高的或过低的，因此认为这个观察结果是随机的。不管n的值是多少，也不管概率ps是多少，都使用同一张c2分布的百分率表，只有自由度v对结果有影响但这张表仅当n足够大时才正确。也就是说，应该做更多次的试验，查看V落入c2分布中的情况如何——通过扩大试验的规模，可以对观察结果的偏倚情况有更好的认识。那么，对于具体的问题而言，n究竟应该取多大呢？事实上，

n的适当选择是不大容易办到的，Knuth提出一个经验规则：使得每个期望值nps至少应该为5，当然如果n比这大得多，会得到更强有力的检验。在使用c2检验时应遵循以下几个步骤：应当做相当多次独立的观察，即n应该足够大;记录落入k个范畴中的每一个观察数，根据式(11-11)或(11-12)计算V.将V的值与c2分布的百分率表中的数进行比较。对于表中v=k-1行，如果V小于1%那列的值或大于99%那列的值，则拒绝这个数并认为它是不够随机的；如果V处于

1%和5%这两列值之间或95%和99%这两列值之间，则这些数是“可疑的”；如果V处于5%和10%这两列值之间或90%和95%这两列值之间，则这些数是“几乎可疑的”……如果V处于25%和75%这两列值之间，则认为这些数随机的。c2检验至少要对不同的数据组进行三次检验，如果三个结果中有两个是可疑的，则这些数被认为是不够随机的。随机数检验一般检验方法c2检验Kolmororov-Smirnov检验(KS检验)

c2检验适用于观察落入有限的k个范畴时的情况，然而当观察结果是无限多个值（例如0和1之间的随机实数）时，KS检验是一种好的选择。假设要描述一个随机量X

的分布，可以使用分布函数F(x)，F(x)

=

Pr(X£x)Pr(X£x)是(X£x)的概率。如果对X做n次独立的观察，并由此得到X1，X2，…，Xn，则有经验分布函数Fn(x)，(

)nF

x

=小于或等于x的X1,X2

,,Xn数目n(11-13)当n充分大时，好的随机数生成程序使得经验分布函数Fn(x)越来越好地近似于F(x)，而坏的随机数生成程序会导致Fn(x)和F(x)有很大的偏离。当F(x)没有跳跃时，可以使用KS

检验。该检验以F(x)与Fn(x)之间的差为基础，定义如下统计量：n

nn

n-¥

<x<+¥-¥

<x<+¥K+

=

n

max

(F

(x

-F

(xK-

=

n

max

(F

(x)-F

(x))(11-14)n+n其中，K

表示

F

(x)大于F(x)n-n时的最大偏离量，K

表示

F

(x)小于F(x)时的最大偏离量；n的作用是放大统计量K+和K-，使其与n

无关。n

np=1%p=5%p=25%p=50%p=75%p=95%p=99%n=10.010000.050000.25000.50000.75000.95000.9900n=20.014000.067490.29290.51760.70711.09801.2728n=30.016990.079190.31120.51470.75391.10171.3589n=40.019430.087890.32020.51100.76421.13041.3777n=50.021520.094710.32490.52450.76741.13921.4024n=60.023360.10020.32720.53190.77031.14631.4144n=70.025010.10480.32800.53640.77551.15371.4246n=80.026500.10860.32800.53920.77971.15861.4327n=90.027860.11190.32740.54110.78251.16241.4388n=100.029120.11470.32940.54260.78451.16581.4440n=110.030280.11720.33300.54390.78631.16881.4484n=120.031370.11930.33570.5453078801.17141.4521n=150.034240.12440.34120.55000.79261.17731.4606n=200.038070.12980.34610.55470.79751.18391.4698n=300.043540.13510.35090.56050.80361.19161.48012n

n+

-与c

检验类似，可以用

K

和

K

分布表中的相应值与实际得到的K

+和K

-的值进行对比，来判断它们是否n

n特别高或特别低。因为x

的个数是无限的，因此式(11-14)不易于在计算机上计算。由于F(x)nn+是递增的，且

F

(x)只在有限步中增长，因此可按如下步骤计算K

和nK-

：对X

做n

次独立的观察，并由此得到X1，X2，…，Xn；对X1，X2，…，Xn

按递增顺序排序，使得X1£X2£…£Xn；[3]用式(11-15)计算K+和K-.n

nnnnn1£

j£nK

=

n

max

j

-F

(X

)

j

j

-1K

=

n

max

F

(X

j

)-1£

j£n

+-(11-15)对于KS

检验来说，观察次数n

的选择也是很重要的：一方面，为保证随机变量X

的分布函数是正确的，应该使n

取很大的值；另一方面，为了检验出序列的局部非随机特性，要求取较小的n；在使用KS

检验时应遵循以下几个步骤：对某个连续分布函数F(x)做n

次独立的观察，给出独立的观察量X1，X2，…，Xn，F(x)必须是没有跳跃的；对X1，X2，…，Xn

按递增顺序排序，使得X1£X2£…£Xn；用式(11-15)计算K+和K-；n

n4.

将K+和K-的值与K+和K-分布表中的数进行比较,确定F(x)和n

n

n

nFn(x)之间偏离的不可能性程度。KS检验和c2检验是两种基本的统计检验，二者是为不同种类的应用设置的：KS检验应用于没有跳跃的分布函数F(x)中，而c2检验应用于跳跃的分布中；但如果把连续的F(x)的区域分成k部分，并忽略每个部分内部的变化，则也可能应用c2检验。虽然如此，但有待检验的是一个连续分布时，KS检验还是比c2检验具有明显的优势。某些情况下二者可以一起使用在对c2

检验总结时提到该检验至少要对不同的数据组进行三次检验，并考虑这些结果有多少是“可疑的”——相比之下有一种更好的方法，可以进行r

次c2

检验，并得到值V1，V2，……Vr，画出这r

个值的经r验分布函数，并把它同正确分布（表11.1）对比；或者通过计算统计量K+和K-，在表11.2

中进行对比以确认K+和K-是否“可疑”。r

r

r随机数检验一般检验方法经验检验方法等分布检验(频率检验)序列检验间隔检验序列相关检验子序列检验等分布检验(频率检验)既可以检验整数序列的随机性，也可检验0与1之间的均匀分布的实数序列的随机性若序列是0与1之间均匀分布的实数序列，对于F(x)=x，0≤

x

≤1，使用KS检验。若待检验的序列是均匀分布在0与d-1之间的整数序列Yj，d是一个方便计算的整数。对于每个整数r，计算Yj=r的次数，0≤r

<d，0≤j

≤n。令范畴k

=d，然后对每个范畴k和概率ps=1/d，应用c2检验，即计算c2统计值，并且在百分率表中进行比较，得出结论。序列检验希望序列中相继的数偶独立地均匀分布序列检验应用于整数序列Yj中相继的数偶，对长度为2n的序列进行n次观察对于0≤j<n，计算(Y2j,Y2j+1)=(q,r)出现的次数，这些计数是对0≤q,r≤d的每个整数偶(q,r)进行的，范畴数k=d2，将c2检验应用于范畴k和ps=1/d2.d的取值方式与在等分布检验中类似，但应该略小，因为一个有效的c2检验应当使n相对于k来说足够大（至少n≥5d2）间隔检验间隔检验用来考察在某个范围内随机数出现的“间隔”长度如果实数α和β满足0≤α<β≤1，考虑连续的子序列Uj，Uj+1

，…，Uj+r

，其中α≤Uj+r<β

，而其余的Up

(j≤p<j+r)不在范围[α,β)内，r+1个数的子序列表示长度为r的一个间隔，或称为Up落入范围[α,β)的间隔，这也是间隔检验得名的原因。算法G

的输入是[0,1)上的实数序列，该序列保存在数组U

中，计算长度为

0，1，…，t-1

的间隔个数，以及长度大于或等于t

的间隔的个数，保存在数组

COUNT

中，直到造出n

个间隔的表格来，作为对序列进行间隔检验的数据。算法G（U，α，β.COUNT）/*生成间隔算法*//*任意给定α和β

的值。*/G1.[初始化]j‹

0.

s‹

0.COUNT[r]‹

0.U[j]<β

THEN

GOTO

G4.FOR

r=0

TO

t

DOr‹

0.G2.[判断U[j]]j‹

j+1.IF

U[j]≥α

ANDG3.[增r]r‹

r+1.

GOTO

G2.G4.[记录间隔长度]/*长度为r

的一个间隔已经找到。*/IF

r≥t

THENCOUNT[t]‹COUNT[t]+1.ELSE

COUNT[r]‹

COUNT[r]+1.G5.[找到n

个间隔了？]s‹

s+1.IF

s<n

THEN

(r‹

0.

GOTO

G3.)▐序列相关检验22jjjjnUnVj

=0

j

=0

j

=0n-1

n

-1C

=

n-1

2

n-12

U

-V

-

j

=0

j

=0统计学中经常出现“相关系数”：若有n个量U0，U1，……Un-1和另外n个量V0，V1，……Vn-1，则它们之间的相关系数定义为n-1

n-1

n-1n(U

jVj

)

-(U

j

)(Vj

)j

=0

j

=0(11-16)我们的讨论排除U0

=U1

=……=Un-1

或V0

=V1=……=Vn-1的情况，因为此时式(11-16)的分母为0.

相关系数总处于(-1,+1)，当它为0

或非常接近0

时，表示Uj

和Vj（0≤j<1）是彼此独立的，当相关系数为–1时，表示完全线性相关。完全线性相关时，对所有的j，存在常数α

和β，使得Vj

=a

–bU

j

.对于随机序列U0，U1，……Un-1，定义下列统计数字为“序列相关系数”()2)222n

U0

1

1

2

n-2

n-1n-1

0

0

1

n-120

1

n-1

0

1

n-1C

=n(U

U

+U

U

++U

U

+U

U

)-(U

+U

++U+U

++U

)-(U

+U

++U(11-17)它测量Uj+1对Uj

的依赖程度。对于一个随机数序列，若式(11-17)中的C接近0，说明序列是充分随机。但实际上，由于UjUj+1

不完全独立于Uj+1Uj+2，所以不能期望C

恰巧为0；一个“好”的C

值将在mn

-2sn

和mn

+2sn

之间n-1m

=n

-12nn2，

s=(n

-1)2

(n

-

2)子序列检验在某些应用中可能一次需要多个随机数，如果一个程序需要r个随机数，则将原来的随机序列U1，U2…分为r个子序列。将应用于原序列<Ui>的检验方法用于这r个子序列进行检验。使用线性同余序列的经验表明，这些子序列的随机性不比原序列的随机性差，除非r和周期长度有共同的大因子。11.3

随机排列在某些应用中，我们会考虑计算N个对象的随机排列问题例如模拟纸牌游戏中，洗一副牌相当于对52张牌的随机排列，因此也把随机排列称为洗牌问题随机排列问题要求等概率地产生N个对象的随机排列例如在洗牌游戏中，52张牌的52!种排列的出现应是等概率的。随机排列一种直观使用均匀分布随机数解决随机排列问题的方法是，产生一个1与N!之间的随机整数，从包含所有可能的N!种排列的表中选择该随机整数对应的那种排列作为随机排列

这种方法在N很小时是很有效的，但当N很大时，N!比某些随机数的精度大得多随机排列我们可以用一种简单的方法来生成随机排列

要产生N!个排列中的任意一个排列，首先将N个元素放在一个数组中

然后生成一个1到N之间的随机整数

进而将该随机数所对应的数组元素与指定位置上的数组元素交换

上述过程重复进行，直到生成N-1个随机整数并进行N-1次交换后，一个随机排列就产生了。算法RDA(X，t)/*随机排列算法*//*设X[1]，X[2]，…，X[t]是要洗的t张牌，UR

算法返回[0,m)区间上符合均匀分布的随机数。*/RDA1.[初始化]j‹

t.RDA2.[生成u]UR(a，c，m，x.u).u‹

u/m.RDA3.[交换]X[k]«

X[j]k‹

j*u+1.RDA4.[减小j]j‹j-1.IF

j>1

THENGOTO

RDA2.

▐算法是否对所有的排列都是等概率的？如果算法中的UR是用线性同余方法实现的，则生成的排列数可能不会超过m个，因此最后的排列完全由所生成的第一个u的值所确定例如，若m=232，则当N‡13时，就有某些排列不会出现了，因为13!≈1.45×232为解决这个问题，我们可以用具有充分长周期的方法实现算法UR例如延搁的斐波那契生成程序，并且使用至少N!个不同的初始值来初始化这个生成程序。随机组合许多应用都要求从N个物品中随机、无偏地抽取出其中的n个例如，质量控制部门或统计计算中要从大量的物品中抽取其中的一部分我们可以将它看成从包含N个记录的文件中随机抽取出n个记录的问题这种抽样问题实质是进行一个随机组合的计算随机组合已经有许多方法用来解决这个问题直观的想法是以概率n/N来选择每个记录，这样能够保证最终抽样的记录数为n个，但结果却并不令人满意只需对这个过程做简单的修改就可以得到满意的结果假定已经在N个记录中考虑了t个记录，并且已经选择了t个记录中的p个，p<t<N，现在正在考虑第t+1个记录是否应该被选择那么对第t+1个记录选择的概率应为(n-p)/(N-t)算法SS(num，n.S)/*简单的随机抽样算法*//*从num

个记录X[1..num]中随机地选择n

个记录，0<n≤num；S

是抽样的结果集合，UR

是一个产生[0,m)上的均匀分布随机数的算法。*/SS1.[初始化]t‹

0.

p‹

0.

S‹

˘

.SS2.[生成u]UR(a，c，m，x.u).u‹

u/m.SS3.[检验]IF

(num-t)*u‡(n-p)THEN

GOTO

SS5.SS4.[选择]S‹

S¨

X[t].t‹

t+1.

p‹

p+1.IF

p<n

THEN

GOTO

SS2.ELSE

RETURN.SS5.[跳过该记录]t‹

t+1.

GOTO

SS2.▐随机组合这个算法产生的抽样完全是无偏的，即以相同的概率n/N选择所有的记录。如果事先不知道N的值为了使用上面的算法，可能要先将整个文件扫描一遍，以确定N的值，然后再执行该算法；一种更好的方法是在第一遍扫描时取出m‡n个原始记录，其中m比N小得多，然后执行算法时仅考虑这m个记录，但这样做的风险是最终的抽样结果不一定是一个真正随机的抽样。随机组合由于不知道输入文件何时结束，因此必须记录下输入记录中被选择的随机抽样即在读输入时，构造一个包含m个记录的“水库”（贮存器），这个水库中只包含前边的抽样中已出现过的记录。如果当前正在考虑的记录被选中，则将它放在该贮存器中。最初，头n个记录总是进入贮存器中，当输入第t个记录且t>n时，根据产生的随机数来判定记录t是否放入贮存器中。建立具有n个索引的一张表，它们指向贮存器中已经选定为随机抽样的那些记录。初始时它们指向前n个记录，当输入记录不断增加时，索引会不断变化。算法RS(S，n.R)/*水库抽样算法。给定n>0，从大小未知但大于n

的文件S

中随机选择n

个记录。“贮存器”M

存放抽样候选者的所有记录，变量m

表示贮存器的大小。索引表I

中的每个索引项指向“贮存器”M

中的一个记录.变量t

表示已经处理过的记录数*/RS1.[初始化]/*输入头n

个记录，并把它们复制到贮存器M

中。*/FOR

j=1

TO

n

DOI[j]‹