随机数的产生和检验

随机数的产生和检验

随机数的产生与检验

摘要

本文通过对常用的随机数的产生方法简单的分析和理论上的验证,对比研究随机数的产生机理以及产生的随机数的好坏,并以此为依据提出自己的一些改进方法，以便对随机模拟更好的利用。

关键词

随机数、随机数的产生随机数的检验

一、引言

随机数的产生方法的研究已经有较长的历史.至今仍有统计学者继续研究随机数的产生的方法和理论.随机数的产生,最早的方法称为手工方法.即采用抽签、掷骰子、抽牌、摇号或者从搅乱的罐子中取带数字的球等方法，许多彩票的发行仍采用这种方法。

随着计算机和模拟方法的应用，计算机来产生随机数成为新的课题。利用计算机产生随机数有两种方法，在计算机内输入随机数表和把具有随机性质的物理过程变换为随机数，如粒子的辐射性，裂变等等。后者得到的随机数均匀性和随机性都很好，而且取之不尽的，但是缺点也明显，对计算的结果不能重复检验，这种物理随机数的产生需要大量的人力物力去检查和维修，成本过高。而数学方法产生的随机数得到了广泛的应用，虽然产生的随机数为伪随机的，正是因为它的占用内存少、速度快、可重复性的优点。

论文

随机数的应用范围很广，对于随机数的均匀性，随机性，独立性的检验也是不可缺少的，只有通过了检验的随机数才有更大的利用空间。本文通过对几种常见的随机数的产生方法进行比较分析，总结其优缺点，并提出一些改进方法。

二、产生随机数的几种常用方法

2.11线性同余法（LCG）

初值 线性同余法通过满足公式(2.1)产生随机序列，主要参数为a, c, M。只有选择合适的参数才能得到随机数的周期接近或达到M。我们把a=137,M=256,c=187用公式(2.1)产生的伪随机数产生方法称为方法T1（见附录1）（周期为256）。类似的，我们把a=1103515245/65536,M=32768(Linux下M=2147483647)，c=12345/65536用公式(2.1)产生的伪随机数称为方法MO（见附录2），它就是我们通常所使用的标准库函数rand。

2.12素数模乘同余法

素数模乘同余法通过满足公式(2.2)来产生随机序列。我们把a=23,M=108+1，满足公式(2.2)的伪随机数产生方法称为方法T2。（见附录3）

2 / 15

2.13线性同余组合发生器（扩大周期法）

利用两个随机数产生器相结合扩大周期的方法。算法如下：

已知两个LCG。

1)用第一个LCG产生k个随机数，一般取k=128.这k个数被顺序地放在

矢量 中。置 。

2)用第二个LCG产生一个随机整数j,要求 ；

3)令 ,然后在用第一个LCG产生一个随机数y， 置.

4)重复第二和第三步，得到随即数列,即为组合同余发生器产生的数列。

若第一个LCG的模为M，令 ，则为均匀随机数列。

这里我们把上述产生伪随机数的方法称为方法D（见附录四）。Gebhardt 1976年证明了这种组合发生器具有随机性强,周期增大的性质,当两个混合式发生器的周期都为 时,组合同余发生器的周期达到了 。

2.14对以上方法总结

MO是默认的使用得最多的产生伪随机数方法，一般在标准C库中

就是使用这种方法，既然它可以被作为库函数使用，必有其独特之处。我们把它作为基准系统来看待。

T1方法可以明显看出周期相比MO缩短了很多，从随机序列行为上

明显有别与MO。

T2方法使用了另外一种实现方法，但其周期比较大。

D方法是速度最慢的，它使用两个产生伪随机数方法相结合的方法，

这两个产生伪随机数方法序列可由上面的产生伪随机数方法产生，集合了上

论文

面几种产生伪随机数方法的特点。

2.2常用的统计性质检验

上面我们介绍了一些常用产生伪随机数的方法，并在宏观上明显地指出它们是有差异的，那么在微观上如何描述一个随机序列的性质呢？Knuth在其经典著作中，详细描述了14种统计检验、12种经验检验以及多种理论检验及谱检验等，并且推断，产生伪随机数生成器的好坏取决于特定的应用。本节将简单介绍几种那些最常用的统计检验。

统计检验的一般方法：

首先假设总体具有某种统计特性，然后由样本值检验这个假设是否可信，这种方法成为假设检验，或成为统计检验，具体的步骤如下：

○1提出假设：总体分布为 ， ；

○2选取适当的统计量 ,其中 为样本求出在 成立时的分布；

○3给定显著性水平，确定检验法，即给出否定域 ： 使

W=1；

○4由观测样本值来计算T值；

○5作统计推断，当T W时否定；当 时，相容。

4 / 15

2.21参数检验

均匀随机数的参数检验是检验由某个发生器产生的随机数序列的均值、方差或者各阶矩等与均匀分布的理论值是否有显著差异。

设 是某个发生器长生的随机数，：样本的特征量与均匀总体的特征量没有有显著差异。首先对特征量作统计检验。在是均匀总体的简单随机样本的假设下，统计量

渐进服从N（0，1）。

对于在上一节中提到的四种方法，我们用Matlab输出每种方法的随机数并做出相关的检验如下表（以U1为例）：

表1

（此表中数据按照公式2.3为手工计算）

论文

通过表中看出上述四种方法全部通过了检验，通过比较u 1或者P值可以给这四种方法按其优点从大到小排序：MO>T1>D>T2。

其中MO是我们通常所使用的标准库函数rand，他是对线性乘同余法中取到了合适的a和M，就可达到很好的效果。另外对于组合发生器，他的有点在于大大增大了伪随机数的周期。

2.22均匀性检验

随机数的均匀性检验又称为频率检验，它用来检验由某个发生器产生的随机数序列是否均匀地分布在[0,1]区间上。也就是是检验经验频率与理论频率的差异是否显著。

设 是某个发生器长生的随机数，假设： 为均匀总体的简单样本。将[0,1)区间分为m个小区间，以表示第i个小区间，设｛ ｝落入第i个小区间的数目为 。

根据均匀性假设， 落入每个小区间的概率为，第i个小区的理论频率为

统计量

（ ）

渐进服从分布。给定显著性水平，差分布表得到临界值后，即可对经验频率与理论频率的差异性做出显著性检验。

表2

6 / 15

（此表中数据计算方法见附录5）在均匀性检验中n=1000,对于T1我们已经知道了其周期为256,故采用n=256来计算T1的统计量，不可思议的很小（n=100，V=8.1）。通过Matlab 计算的结果V，可以比较V的大小来看出四种方法的优劣，其中V越小说明均匀性越好。上述四种方法全部通过了检验，通过比较，T1，即线性同余的方法产生的均匀性比较好，这也是取特殊的a,c,M 使得产生的随机数比较好。

2.23独立性检验

两个随机变量的相关系数反映他们之间线性相关程度，若两个随机变量独立，则他们的相关系数比为零（反之不一定），故可以利用相关系数来检验随机数的独立性。

设 是一组待检验的随机数，假设：相关系数 。考虑样本的j阶自相关系数

（ ）

当n-j充分大，且 成立时， 渐进服从N（0，1）分布；在实际检验中常取。利用统计量 可以进行相关性检验。

论文

表3

（此表中数据计算方法见附录6）在上表中只是简单通过j阶的相关系数检验来说明这四种方法的独立性好坏。当然此方法有些简单，但它直观的说明了这四种方法产生的随机数都有很好的独立性。

三、结束语

通过对上面常用的四种方法的分析和检验，特别是三种检验中，不同的方法产生的随机数对于随机性、均匀性、独立性等等的统计性质有一定的差别，这种方法或许它的随机性很好，但是均匀性或者独立性却不太好；或许某种方法它的独立性很好，但是随机性却不算好。所以呢，我们在使用这些方法产生随机数的时候，我们可以更具某个具体问题它针对的是那个方面，对于某一方面比如随机性要求很高的话，我们就选取产生随机数随机性强的来产生，而其他的方面我们可以选取不同的方法来产生数据。

另外一方面，对于同一种方法选取的初值，模长，增量，倍数以及产生随机

8 / 15

数的数目，这五个方面在选取适当的情况下，统计性质会得到很好的体现，我们在T1 和MO这两种方法的比较中很明显的看出来。所以呢，我们在应用这些方法产生随机数的时候也要注意上述五个方面的影响，在随机模拟的时候要充分发挥数学方法产生随机数的可重复性，选取不同的合适的初值，模长等使得我们的随机模拟的精度精益求精。

参考文献

本文的大部分公式，定理来源于北京大学出版社出版，高惠璇编著的《统计计算》。

附录：

附录1

论文

T1方法的程序：function r = T1(x0,n)

r = zeros(n,1);

x = zeros(n,1);

M = 256;

c=187;

a =137;

x(1) = a*x0 + c;

r(1) = x(1)/M;

for i=2:n

y = a*x(i-1)+c;

x(i) = mod(y,M);

r(i) = x(i)/M;

end

附录2

MO方法的程序：function r = MO(x0,n) M = 32768;

a =1103515245/65536; x(1) = 4*x0 + 1;

r(1) = x(1)/M;

10 / 15

for i=2:n

y = a*x(i-1);

x(i) = mod(y,M); r(i) = x(i)/M;

end

附录3

T2方法的程序：

function r = T2(x0,n)

format long;

M = power(10,8)+1;

a = 23;

r = zeros(n,1);

x = zeros(n+5,1);

x(1) = x0;

for i=2:n+5

y = a*x(i-1);

x(i) = mod(y, M);

end

r = x(6:(n+5))/M;

format short;

附录4

组合发生器D方法的程序：function x=zuhe(m)

论文

a=1103515245/65536; M=32768;

t(1)=12345/65536;

n=1;

for i=1:m

t(i+1)=mod(a*t(i),M);

r(i)=t(i)/M;

end

k=random('unif',1,m,1); j=floor(k);

x(1)=r(j);

for n=2:m