三、伪随机数产生方法

3.1平方取中法

平方取中法是冯·纽曼(John van Neumann)40年代中期提出的。这个方法首先从某个初始的种子数开始,求出这个数的平方。取这个平方数的中间几位作为随机数序列中的第2个数;再求出第2个数的平方,又取这个平方数的中间几位作为随机数序列中的第3个数;不断按这个方式继续此算法,即可得到相应的伪随机序列。

1:利用平方取中法产生4位数的随机数序列,序列的种子数取为x0= 3187,通过计算得到

x0= 3187

(3187)2=10156969     x1=1569

(1569)2=02461761     x2=4617

(4617)2=21316689     x3=3166

(3166)2=10023556     x4=0235

(0235)2=00055225     x5=0552

(0552)2=00304704     x6=3047

(3047)2=09284209     x7=2842

将此过程继续下去,还可以得到076959139635833242228542,…,

这个方法有几个缺点。首先,利用这个方法产生的伪随机数序列的重复周期通常较短。第二,对于较长的伪随机数序列,利用这种方法可能无法通过随机性的统计检验。第三,当在任何时候生成之后,其后产生的数都将为0。如果这种现象在一个较复杂的仿真研究过程中出现,它将会使仿真分析人员误入歧途。这种现象如下例所示。

2:利用平方取中法产生两位数的随机数序列,种子数取为x0=44。通过计算得到:

x0=44

(44)2=1936     x1=93

(93)2=8649     x2=64

(64)2=4096     x3=09

(09)2=0081     x4=08

(08)2=0064     x5=06

(06)2=0036     x6=03

(03)2=0009     x7=00

(00)2=0000     x8=00

这样,随机数就无法继续产生了。

利用平方取中法的另一个问题是这个方法可能产生退化,即总是得到相同的xi值。

3:设在产生四位数的随机数过程中,得到了一个xi值为4500,即

                xi=4500

从而,                  (4500)2=20250000  xi+1=2500

            (2500)2=06250000  xi+2=2500

相继产生的数值总为2500

由于这些原因,平方取中法已经被许多新的能提供更可靠的随机数序列的算法所取代。

3.2线性同余法

线性同余法在1951年由菜默尔(Lehmer)首先提出。目前大多数随机数发生器都采用这种方法。在这个算法中,随机数序列中的数由如下的递推关系产生   

                        n0             (1)

初始值x0称为种子,常数a称为乘子,常数c称为增量,而常数m称为模数。

对于(1)式,当c=0时,该算法称为乘同余法;当c0时,该算法称为混合同余法。从例4中可以看出利用这种方法产生的序列的重复性,一般来讲任何由此方法产生的序列都存在重复性。在大多数情况下,合理地选择常数a,c,x0m,可以使重复周期充分的长。

(1)式中,显然可以得出

0xn m-1                           2

为了得到[0,1]区间上分布的随机数,可以令

                                                      3

Rn为满足要求的随机数。

常数m,a,c的值对所产生的随机数序列的周期长度有很大的影响。

4:设a=5c=3m=15,取x0=7,利用线性同余法产生随机数序列。

x0=7

x1=(5*7+3)mod16=6    R1=0.375

x2=(5*6+3)mod16=1    R2=0.063

x3=(5*1)+3mod16=8    R3=0.500

x4=(5*8+3)mod16=11   R4=0.688

x5=(5*11+3)mod16=10  R5=0.625

x6=(5*10+3)mod16=5   R6=0.313

x7=(5*5+3)mod16=12   R7=0.750

x8=(5*12+3)mod16=15  R8=0.938

x9=(5*15+3)mod16=14  R9=0.875

x10=(5*14+3)mod16=9    R10=0.563

x11=(5*9+3)mod16=0     R11=0.000

x12=(5*0+3)mod16=3     R12=0.188

x13=(5*3+3)mod16=2     R13=0.125

x14=(5*2+3)mod16=13    R14=0.813

x15=(5*13+3)mod16=4    R15=0.250

x16=(5*4+3)mod16=7     R16=0.438

x17=(5*7+3)mod16=6     R17=0.375

x18=(5*6+3)mod16=1     R18=0.063

x19=(5*1+3)mod16=8     R19=0.500

n=16时出现循环。

对于常数的选择讨论如下:

  m的选择

由于重复周期的长度总是小于m,因此需要将m取大的数值,更进一步,所选用的m的值应能简化同余关系的解,在计算机中数字都是用二进制表达的,因此已经证明m取值为(2k)是很好的。

  ac的选取

当且仅当下列条件满足时,一个由线性同余法产生的随机数序列的最大可能重复周期为m

l         cm互质,即同时能被cm整除的正整数只有1

l         如果m能被4整除,则(a-1)也能被4整除,即a=1+4k

特别地,当选择a=216+5=65541a=216+3=65539时可以得到满意的结果。至子c的选择,只要满足cm互为质数的条件即可。

  x0的选取

   如果随机数序列的周期为m,因为能产生完全的序列,即在一个周期内可以取到0至(m-1)的所有值,因此x0的选取是不重要的。但仍然要小心,例如取x0=0时会产生退化的序列。

对于乘同余法,由于c=O,无论怎样选择m,都无法满足cm互质的条件,因而不可能得到满周期。若选择m=2k,则所产生的随机数序列的周期p2k-2,即在0m-1之间的整数至多只有四分之一可能成为xn的值,而且这四分之一的整数在0(m-1)之间是如何分布的尚难确定。这与种子数x0的选取有关。若取乘子为a=8L3a=8L+5形式的整数,种子x0取为奇数.则可以达到最长的周期p=2k-2

5:使用乘同余法,对a=13m=2^6=64,x0=1234,求产生器的周期。可以看出当种子为13时,序列的周期都为16。而当种子为24时,序列的周期分别为84

i

X0=1

X0=2

X0=3

X0=4

0

1

2

3

4

1

13

26

39

52

2

41

18

59

36

3

21

42

63

20

4

17

34

51

4

5

29

58

23

 

6

57

50

43

 

7

37

10

47

 

8

33

2

35

 

9

45

 

7

 

10

9

 

27

 

11

53

 

31

 

12

49

 

19

 

13

61

 

55

 

14

25

 

11

 

15

5

 

15

 

16

1

 

3

 

从(1)式中得到的序列χ0,χ1…实质上完全不是随机的,因为若设     

                     4

则有

                         =

                        

                         =

即一旦m,a,c,x0确定,xi就完全被唯一地确定下来了。而且在[0,1]区间上得到的数值Rn最多只能取到m个值。

 

3.3加同余法

加同余法需要n个数的序列作为它的种子,这n个数的序列可以应用其它的方法产生,应用加同余法可以使这个序列不断扩大,加同余法的算法是

                                     5

这种方法的主要好处是速度快,因为它不需要作乘法运算。它可以得到大于m的周期但这种方法产生随机数的过程不象混合同余法那样清楚,因此由此方法产生的随机数序列需要经过仔细的确认。

3.4二次平方同余法

二次平方同余法适用于m2的幂次的情况,这种方法的递推关系式为

            n0                     (6)

种子数x0必须满足条件