第四章 随机数产生原理.ppt

1、第四章 随机数产生原理,一、引言 二、伪随机数产生原理 三、0，1均匀分布随机数的算法 四、其他分布随机数的产生 五、正态分布随机数的产生 六、MATLAB统计库中的随机数发生器 七、随机数的检验 八、案例3 九、习题,一、引言,以随机数产生为基础的Monte-Carlo方法已成为现代科研的重要手段之一。其意义早以超出了概率论与数理统计的范畴。广泛应用于计算方法、随机数规划、管理科学、物理化学中的高分子结构的研究等领域。我们来看一些例子。,1、数值计算的研究,数值计算的研究可以说是一切Monte-Carlo应用的基础，在计算数学领域我们遇到了很多的复杂计算，一个典型的例子是计算积分。对于一维、

2、二维的问题早已获得解决。但是当遇到高维积分问题时，我们传统的数值方法都由于计算量太大而陷于了困境。但是高维积分问题又偏偏是物理、高分子化学、运筹学和最优化问题迫切而必须解决的问题。我们来看一个例子。,这里,这是一个众所周知的积分公式，我们当然也可以把它一般的看为是一个高维积分，如果从传统的数值计算方法来看待，则高维问题是随着维数的增加计算量成指数增加的，计算很快就失去控制。但是如果我们换一个角度来看待这个问题，从概率论的角度，它实际是：,即是f（x）的均值，对于均值我们有一个很好的估计，即,【例4.1.1】 用Monte-Carlo 对 积分,解：将积分区域和值域看成是一个边长为一的正方形。利

3、用均匀分布随机数将点撒在正方形中，计算小于函数的个数并除全部点数。这就是积分的近似值。,% 利用Monte-Carlo方法计算定积分 x=rand(1,1000); x_2=x.2; JF=mean(x_2) % 作Monte-Carlo积分示意图 for i=1:1000xx=rand(1,100);yy=rand(1,100); end x1=linspace(0,1,50); y1=x1.2; plot(xx,yy,.,x1,y1,linewidth,2) axis equal h = legend(x-y,x2); JF = 0.3346,面积计算结果为： s = 0.3482,【例4

4、.1.2】 利用Monte-Carlo方法计算定积分。,解：抽两组随机数，求每组元素的平方代入给定的函数，然后求平均值即得积分的近似值。 % Monte-Carlo方法积分二重积分并与数值方法的结果进行比较 Q = dblquad(sin(x.2+y.2),0,1,0,1) % 数值求积分命令 x=rand(2,100000); % 产生两组随机数 Sin_xy= sin(x(1,:).2+x(2,:).2); % 代入函数 JF_Sin_xy=mean(Sin_xy) % 用Monte-Carlo方法求积分值,计算结果为： Q = 0.5613 JF_Sin_xy = 0.5612 当抽样数

5、很大时结果很接近。我们可以从Monte-Carlo方法中看出，如果维数增加实际计算难度并不增加，因此是处理高维问题的有效方法。,这里 x 是积分定义域中的均匀分布随机数，这是革命性的一个视角。从这个视角，我们把繁难的积分计算变成了简单的算术平均，而大量的抽样计算又是计算机的拿手好戏，更重要的是当维数增加时并不增加计算难度，从而用 Monte-Carlo 方法研究高维积分问题已是当今计算数学界的热门课题。,2、管理科学的系统仿真研究,管理科学中的系统仿真研究，如服务系统、库存系统等。其共性就是研究的对象是随机数，如顾客到达时间一般是一个服从指数分布的随机数，而服务时间也可以看成是服从某种分布的随

6、机数，当一个系统是多队多服务体系时，问题就变的相当复杂了。我们很难用数学的解析式来表达。这时Monte-Carlo方法也是有利的武器。对于这个领域的已有各种比较成熟的专用软件如GPSS、SIMULATION等可以使用。,3. 物理化学中的分子领域,50年代科学家已经在高分子领域使用Monte-Carlo方法了。这一领域所研究的问题更加复杂，计算量非常之大。高分子材料是由几乎是无穷的高分子链组成，而每一个链的长度又是10的好几次方。而链的构象又是千差万别，而且是随机游动的。如何在中微观上几乎是无规律的现象中去判断其宏观的性质？用数学的解析式来说明这样的现象是苍白无力的。Monte-Carlo方法

7、是一个很好的工具，它使得科学家用Monte-Carlo方法去探索高分子运动的规律。一个典型的例子是：对于高分子多链体的研究这是一个很复杂的问题，直到标度理论和重正化群理论方法的引入，才使得单链构象统计问题得到了较好的解决。,例：用计算机模拟高分子链,链的末端距,末端距：空间一链的末端与始端的距离为末端距，由于我们将始端放在坐标原点，所以末端距的 计算公式为：末端距=（X2+Y2+Z2)1/2 这里X，Y，Z为链的末端点的坐标。显然末端距随链的不同而不同，即为随机变量。,三根链的起点 (0,0,0),蓝链的末端距,绿链的末端距,红链的末端距,二、伪随机数产生原理,前面Monte-Carlo方法的

8、例子是以高质量的随机数为基础的。通过完全的随机抽样或调查可以产生随机序列。当我们用Monte-Carlo方法研究一个实际问题时，我们需要快速地获得大量的随机数。用计算机产生这样的随机数是非常方便的，用数学方法在计算机上产生的随机数称为伪随机数。,基本定理： 如果随机变量X的分布函数F（x）连续，则R = F（x）是 0，1上的均匀分布的随机变量。 证：因为分布函数F（x）是在（0，1）上取值单调递增的连续函数，所以当X在（，x）内取值时，随机变量R则在（0，F(x)）上取值，且对应于（0，1）上的一个R值，至少有一个x满足，见图4.2.1,以 表示随机变量 R 的分布函数，则有,(4.2.1)

9、,=,证毕,图4.2.1,(4.2.2),基本定理给出了任一随机变量和均匀分布R之间的关系。而有些随机变量可以通过分布函数的逆变换来获得，因此如果我们可以产生高质量的均匀分布，我们就可以通过变换获得高质量的其他分布。见公式 (4.2.3),(4.2.3),例4.2.1 求指数分布的随机数。令,从而我们用服从0，1上的随机数R，通过上面的公式就可以得到指数分布的随机数了。,例4.2.1 产生1000个均匀分布随机数，利用变换产生=6的指数分布并进行拟合优度检验。,clc,clear x=linspace(0,20,100); R=rand(1,1000); % 产生1000个（0，1）均匀随机数

10、 ex=-6*log(1-R); % 变换为指数分布随机数 ex=ex; m=mean(ex) v=var(ex) subplot(1,2,1),cdfplot(ex) subplot(1,2,2),hist(ex) kstest(ex, ex expcdf(ex, 6) % 拟合优度检验,结果为： H=0, 接受原假设，变换后的确为=6的指数分布,三、 （0，1）均匀分布随机数的产生,1、算法要求 （1） 产生的数值序列要具有均匀总体简单子样的一些概率统计特性，通常包括分布的均匀性，抽样的随机性，试验的独立性和前后的一致性。 （2） 产生的随机数要有足够长的周期。 （3） 产生随机数速度快，

11、占用内存小。,为了达到快速的要求，一般采用递推公式,（4.3.1),目前最常用的方法是上述方法的一个特例：,混合同余法,(4.3.2),其中a，b，M 以及初值 y 都是正整数，容易看出 x 满足0x1。其中 mod M 运算定义为：任一整数 y 可唯一表示为公式,则,乘同余法,当 b = 0 时，有,(4.3.4),加同余法,以下形式的同余法称为加同余法,(3.4.5),例4.3.1 历史上比较有名的称为“菲波那西”数列为加同余法的特例。,(4.3.6),当 M=8 时，取初值得“菲波那西”数列。0，1，1，2，3，5，8，13，21，34，55，89，144，253 对上述数列取模得：0，

12、1，1，2，3，5，0，5，5，7，1，1 (4.3.7) 再除以模 M 我们可得到 (0,1) 之间的序列 。,我们知道对于一个来自均匀分布的随机序列，应该满足独立性、均匀性等统计特性，但伪随机数往往有一些缺陷，例如 (4.3.7) 序列到一定长度后，又开始重复下面的序列，这称为周期性是一种明显的规律，与随机性矛盾。通常我们只能选用一个周期内的序列作为我们的伪随机数。因此研究一种算法，使得其产生的随机序列的周期尽可能长，我们可以通过调节（4.3.1）的参数来实现。 因此如何来获得一个周期比较长的序列，就成了我们研究的一个内容：有关伪随机数序列的周期有如下的一些结论：,定理 4.3.1 混合同

13、余法产生序列达最大周期 M 的充要条件：（1） b 与 M 互素（2） 对于 M 的任意素因子 p，有（3） 如果 4 是 M 的因子，则显然乘同余法产生的序列达不到周期 M（不满足（1）。当 取 （k为任意整数）时，因为 M 只有一个素因子2， 且4是 M 的因子，则由条件（2）、（3）有 ， 从而混合同余发生器达到最大周期的算法为：,（3.4.8),其中c，d为任意整数。混合同余发生器是否达到最大周期M与初始值无关。,对于乘同余发生器，由同余运算的定义，知其由如下性质,(1) 如果,则有：,（2）如果,则,其中（c，M）是c，M的最大公约数。,利用这些性质可得到以下定理。 定理4.3.2

14、对乘同余发生器，若 ，则使成立得最小正整数 V 就是此发生器得周期。,在数论中称 V 为 a 关于 M 的阶数，对于乘同余发生器，若种子与 M 互素，则其周期就是关于 M 的阶数。这样一来，寻找达到最大周期的同余发生器的问题就转化为数论方面寻求 M达到最大阶数 a 的问题了。Knuth 对这一问题的研究作了总结。,从算法上我们知道公式是递推的，因此一般的随机发生器程序都要预先赋初值，这种初值为种（Seed），有些统计软件如SPSS要求用户给出Seed. 以均匀分布（0，1）随机变量R变换成的随机变量。以r,u,分别表示 (0,1) 均匀分布,指数分布，N(0,1) 标准正态分布。其他常见的分布

15、如卡方分布、F分布等的抽样方法见表4.3.1。,四、其他分布随机数的产生,1、 直接抽样法由基本定理我们知道，对于有些随机变量可以建立与R的函数关系，因此只需对R进行抽样，利用函数的映射关系我们就可以方便地得到该随机变量的抽样了。如前面的指数分布随机数。 2、变换抽样产生随机变量的变换抽样方法，是讨论均匀分布的不同函数分布，为随机变量抽样提供一些简单可行的算法。在概率论中，从不同的角度出发，对随机变量函数进行了讨论，以下列出一些结果。,设随机变量 X 具有密度函数 ，是对随 机变量 X 的变换，且 的逆函数存在，记为 有一阶连续导数，则随机变量的密度函数,，,（4.4.1),例4.4.1 R，

16、1-R均为（0，1）上的均匀分布随机变量，设随机向量（X，Y）具有二维联合密度。对于随机变量 X，Y 进行函数变换,其中，函数 ， 的逆变换存在，记为,且存在一阶偏导数，设 J 为 Jacobi 变换,则随机变量的二维联合密度为,(4.4.2),例4.4.2 用变换抽样产生标准正态分布的随机变量 U 随机变量U 的密度函数为:,取随机数,，则,相互独立、服从N（0，1）分布。解上面的两个方程，得逆变换公式,由(4.4.2)，随机变量 的密度函数,从而我们知道是两个独立的标准正态分布。 下面还有几个常用的二维函数的变换结果。1） 随机变量和 的密度函数,（4.4.3),2) 随机变量的积 的密度

17、函数,(4.4.4),3) 随机变量的商 X/Y,(4.4.5),3、值序抽样 值序统计量 ，通常称为值序量。是随机向量的一个函数，取其一个实现并排序得其中的第 l 个值为函数值。显然这是一个统计量,若随机变量的各个分量独立且同分布，则值序统计量 的密度函数 和分布函数分别为：,这里，F(x), f(x) 分别为随机变量在分布函数和密度函数。,特别当随机变量 为0，1上的均匀分布时， 得密度函数为：,这是 分布的密度函数,因此我们可以很容易地产生特殊的 分布的随机数。,例4.4.3 产生服从 分布的随机数 若随机变量 的密度函数为：,利用（0，1）均匀分布随机数我们可以如下产生特殊 随机数,（

18、1）,（2）,（3）,（4）,我们来验证（2），即是否服从a=1，b=5的贝塔分布，按公式抽5个均匀分布的随机数，取其中最小的为一个样本，共取1000个，然后用分布的Kolmogorov-Smirnov拟合优度检验，程序如下：,% 例4.4.3验证(2) % 抽5个均匀分布随机数，产生一个服从bata分布的随机数 for i=1:1000B(i)=min(rand(1,5); end B=B; % 转置为列向量 h=kstest(B, B betacdf(B,1,5) % 检验服从bata分布吗？,计算结果为： h=0 接受原假设，服从a=1，b=5的贝塔分布。学习者可以验证其余的抽样公式。,

19、五、正态分布随机数的产生,正态分布在数理统计中具有基础性的作用，因此产生高质量的正态分布有重要的 意义，这一节我们将介绍几种数值方法求正态分布，利用中心极限定理产生正态分布的随机数。,1） 利用中心极限定理产生随机数中心极限定理：,设,服从均值为，方差为2 的某一分布，令,（4.5.1),则当 n 充分大时， 渐近地服从标准正态分布N（0，1） 注意，这个定理没有指出随机变量x 是服从什么分布的。这 正是该定理的神奇之处。我们现在已经能产生0，1均匀分 布的随机数了，那么我们可以利用这个定理来产生标准正态 分布的随机数。,现在我们产生n个0，1均匀分布随机数，由公式（4.5.1）我们有：,为编

20、程方便，我们特别选 n = 12，则,这样我们很方便地就把标准正态分布随机数计算出来了。,例4.5.1 利用中心极限定理产生标准正态分布随机数并检验,% example 4.5.1 clc,clear for i=1:1000 R=rand(1,12); X(i)=sum(R)-6; end X=X; m=mean(X) v=var(X) subplot(1,2,1),cdfplot(X) subplot(1,2,2),histfit(X) h=kstest(X, X normcdf(X, 0,1),结果为： H=0, 接受原假设，变换后的确为标准正态分布。,2) Hasiting 有理逼近法

21、,这是一种计算速度快，也能满足一定精度的算法。我们可以构造分布函数反函数的近似逼近公式，来产生标准正态分布的随机数。其计算公式为：,（4.5.2),这里 ，系数为：,a0 = 2.515517 b1 = 1.432788a1 = 0.802853 b2 = 0.189269a2 = 0.010328 b3 = 0.001308,六、利用统计工具箱产生随机数,在MATLAB统计工具箱中为我们提供了大量的产生各种随机数发生器程序，我们只需要调用就可以产生我们想要的随机数。,随机数发生器,例4.5.2,例4.6.1 设正态二维分布的密度函数为：,作二维分布的散点图和直方图,% 例4.6.1 二为独立

22、正态分布的散点图和直方图 mu = 0 0; sigma = 1 0; 0 1; r = mvnrnd(mu,sigma,1000); subplot(1,2,1),plot(r(:,1),r(:,2),+) subplot(1,2,2),hist3(r,10 10),七、随机数的检验,我们已经基本搞清伪随机数的产生原理，由于并不是真正的随机数，很自然的问题是，它们是否具有真正随机数的那些统计性质如参数大小、独立性，均匀性等等。 设：随机数具有连续的分布函数F（X），则随机变量R=（X）是均匀分布（0，1）的随机变量，因此如果R通过统计检验随机变量 X 也可以通过。因此我们以下着重讨论均匀分布

23、R的检验问题，再简单地讨论正态随机数检验问题。,统计推断原理 统计量的定义： 设 为 随机变量序列，则随机序列的函数称为统计量。记为：,显然统计量 S 也是随机变量。既然是随机变量，它们就应该有其分布或称总体的规律，当然也有各种数字特征。例如均值、标准差、方差等等各阶矩。 我们的统计推断方式是： （1）H0：某假定成立 （2）在假定成立的条件下构造统计量S,（3）统计量构造完毕，我们也就知道了该统计量的全部统计规律。如它的分布函数，或密度函数各阶矩等。 （4）根据统计量的分布，在给定的显著性水平， 对统计量S 的一次抽样确定以 1-为概率的区域，该区域称为接受域 。如果该次抽样计算出统计量 S

24、 的值 s 落入该领域，我们就接受原假，否则推翻原假设。 这个就是小概率事件在一次实验实际不可能发生原理。 落入由 确定的区域是一个小概率事件，在一次实验中我们认 为是不可能发生的。,接受域,1、 统计检验中两类常用统计量的构造检验方法 （1）设随机变量 X 具有数学期望E（X）= ，和有限方差 D（X）= ，我们抽 N 次样本。X1，X2，XN，当N充 分大时，统计量,(4.6.1),以 N（0，1）为极限分布，取显著水平 = 0.05, 则拒绝为 （，1.96），（1.96，）。 当 1.96，则认为差异显著，拒绝假设 E（X）=,（2）将样本X1，X2，NN，按一定规则分为不相交的 K组

25、，记 i 组的观测频数为 ni（i=1,，k），若随机变量 X 落于弟 i 组的概率为 Pi，则得理论频数m i= N Pi，由ni，mi 构造统计量。,=,(4.6.2),渐近服从自由度为 的分布，简记 这里， l 是确定概率 P 时，由子样 X 中给出的约束条件数。为有效地进行 统计检验，一般要求（4.6.2）的样本数 N30；在（4.6.2）中k5，mi5。,当统计量的自由度 时， U = 渐近服从N（0，1）分布。 服从柯尔莫格洛夫斯米尔诺夫的统计量,取显著性水平= 0.05, 可拒绝原假设。,我们可以用这种方法进行独立性，及拟合优度检验。,2、 随机数的统计性质检验,设均匀分布的随机

26、样本为：（4.6.3） 是一组需要检验的0，1均匀分布随机数，利用(4.6.1)， (4.6.2) 我们构造统计量来对常用的参数进行统计检验。,1）参数检验 参数检验是检验参数估计值与理论值的差异是否显著的方法，设我们有0，1 均匀分布样本R1，R2，RN （4.6.4） 我们可以构造以下统计量（参数估计量）,即随机变量R的一阶距，二阶距和方差这些参数的估计量。根据随机变量R的理论分布，不难计算,现在我们可以应用中心极限定理及（4.5.1）来构造统计量了,渐近服从N（0，1）分布，从而可以推断样本估计值与理论参数的差异。我们还可以利用同样的方法去构造三阶矩、四阶矩，对随机变量的偏度和峰度进行参

27、数检验。,例4.7.1 对随机数发生器unifrnd产生的1000个随机数进行均值、方差、偏度和峰度等的参数的检验。 偏度计算方法为：,u3=mean(R-R_mean)/R_std).3)*0.408248*sqrt(n); 峰度计算方法为： uu=mean(R-0.5)/sqrt(1/12).4)-1.75 u4=uu*0.204124*sqrt(n);,function s1,s2,s3,s4=moment_test(R) % 对（0，1)均匀分布随机数进行矩检验 n=length(R); R_mean=mean(R);R_var=var(R);R_std=std(R); u1=sqrt

28、(12*n)*(R_mean-0.5); if abs(u1)1.96s1=pass; elses1=*; end,% 对方差进行检验 var(R) u2=sqrt(180*n)*(R_var-1/12) if abs(u2)1.96s2=pass; elses2=*; end % 对偏度进行检验 u3=mean(R-R_mean)/R_std).3)*0.408248*sqrt(n); if abs(u3)1.96s3=pass; elses3=*; end,% 对偏度进行检验 uu=mean(R-0.5)/sqrt(1/12).4)-1.75 u4=uu*0.204124*sqrt(n);

29、 if abs(u4)1.96s4=pass; elses4=*; end,调用该函数的主程序为：R=rand(1,1000); s1,s2,s3,s4=moment_test(R); S1 % 显示对均值的推断，pass为接受，*拒绝 S2 % 显示对方差的推断，pass为接受，*拒绝 S3 % 显示对偏度的推断，pass为接受，*拒绝 S4 % 显示对偏度的推断，pass为接受，*拒绝,计算结果为： s1 = pass s2 = pass s3 = pass s4 = pass,2) 均匀性检验或称分布的拟合优度检验 均匀性检验，又称频率检验，检验它的经验分布频率和理论频率的差异是否显著。

30、把0，1区间分为k个等区间，按 ri 取值的大小把（4.6.3）分为 k 组，设有ni个随机数属于i组，即共有个满足不等式（i1）/kri/k。根据均匀性假设，满足在每个小区间的概率为Pi=1/k，mi=N/k由（4.6.2）得统计量,渐近服从 分布，而统计量,(4.6.8),理论频数ni,mi,为累积频率检验，渐近服从柯莫洛戈洛夫斯米尔诺夫分布 取显著水平 = 0.05，当DN1.35 时拒绝均匀性假设。,（4.6.9),例4.7.2 对分布的拟合优度进行检验。 （1）渐近服从柯莫洛戈洛夫斯米尔诺夫检验,R=unifrnd(0,1,1000,1); h=kstest(R, R unifcdf

31、(R, 0, 1),结果为：h=0， 接受原假设，是来自（0，1）均匀分布随机数,R=unifrnd(0,1,1,1000); % 构造卡方统计量 k=12; n=length(R); n1=hist(R,k); % 计算每个区间的频数 kf_7 = k/n*(sum(n1-n/k).2); % 计算分位点即统计量 chi2_p=chi2cdf(k-1,kf_7); % 计算下侧概率 if chi2_p0.95chi2_str=pass; elsechi2_str=*; end chi2_str,（2）利用卡方分布进行检验,3）独立性检验 随机数独立性检验，其重点是检验（4.6.3）中前后各随

32、机数的统计相关是否异常。我们知道统计性能良好的随机数数，前后之间应该是统计独立的，若存在明显的相关性则视这批随机数为不合格。,这里： 表示i 阶滞后相关系数，为0时表示没有相关性。,（1）相关系数检验 相关系数取值为零是两个随机变量独立的必要条件，其值的大小给出它们之间线性相关强弱的测度，可以用来检验随机数的独立性。例如随机序列与其自身持滞后序列的相关性。 例如： 一阶滞后相关性,j 阶滞后相关性,利用(4.6.3)数据构造统计量,（4.6.10),对充分大的N（如Nj50），取零假设则统计量,（4.6.11),渐近服从N（0，1）分布。,例4.7.3对17阶滞后相关系数进行检验。,funct

33、ion sacf1,sacf2,sacf3,sacf4,sacf5,sacf6,sacf7=acf_test(R) % 独立性的自相关AFC检验 R_mean=mean(R);R_var=var(R);n=length(R); for i=1:7 rou(i)=sum(R(1:n-i).*R(i+1:n)-R_mean2)/R_var)*sqrt(1/(n-i); end rou if abs(rou(1)1.96sacf1=pass; elsesacf1=*; end if abs(rou(2)1.96sacf2=pass; elsesacf2=*; end,（2）联列表检验 在XY平面上，

34、将单位正方形分为个相等的小正方形，把随机数（4.6.3）按出现的先后顺序两两分组，如取,记落入小正方形 ( i, j) 内的观测数为 nij ( i, j = 1, 2, , k)，令,H0：数据（4.6.3)来自均匀分布随机数,则可以构造统计量,其极限分布是服从自由度为（k-1)2 的卡方分布，记为,(4.6.12),例4.7.4 对随机数进行联列表检验。,R=unifrnd(0,1,1000,2); N=2000 k=6 n=hist3(R,k k) % 产生每个小正方形落入的个数 ni=sum(n); nj=sum(n); nij=ni*nj; n_sum=sum(sum(n.2./ni

35、j)-1 chi2_2=N*n_sum chi2_p=chi2cdf(k-1)2,chi2_2); if chi2_p0.95chi2_str=pass; elsechi2_str=*; end chi2_str,计算结果：pass,4) 组合规律检验 该检验用于对数据（4.6.3）的组合规律进行检验，按随机数先后出现的顺序，根据一定的规则把这些随机数组合起来，检验他们的组合规律的样本性质与理论性质的差异，并进行判断。（1） 扑克（颜色）检验 将（4.6.3）的数据按顺序分为8个随机数一组，记组数为L，对每组的随机数取其小数点后第一位，并以八进制来表示。,例如：如下一组数1 2 3 4 5 6

36、 7 8 (0.001, 0.615, 0.516, 0.316, 0.815, 0.95, 0.11, 0.216) 小数点后面第一位用 8 进制表示，我们有（ 0, 6, 5, 3, 0, 1, 1, 2 ） 一般地，我们们记为显然有,当 中有 k 种相同数字时，称该序列为 k 色向量，例如: （1，1，3，4，4，6，7，0） 为6色 （5，5，5，5，5，5，5，5） 为1色 （0，2，5，4，7，6，1，3） 为8色将向量按颜色的不同分为 8 类 。从理论上我们有以下频率。,记 为落入第 k 组向量的个数， 为理论频数，则我们可以构造统计量如下：,(4.6.13),例4.7.5 对1

37、000个随机数进行扑克牌检验。 解题思路： 1）对1000个随机数按顺序每8个一组，共分125组。 2）分别对每组元素小数点后第一位进行模为8的运算。 3）计算每组的颜色数，并放入数组cc 4）构造自由度为4的卡方分布统计量进行检验。,% 扑克牌检验，2005.10.24 clear,clc R=rand(1,1000); % poke test n=length(R);kk=8; jj=1; % 求125组数据的颜色数，并放入矩阵cc中。 for ii=1:kk:n % 对每组数小数点后第一位取模为8运算 rr=10*R(ii:ii+7); pk=mod(fix(rr),8);,% 计算每组

38、的颜色数 pk=sort(pk); j=1; for i=1:7if pk(i)=pk(i+1);j=j+1;end end cc(jj)= j+1;jj=jj+1; % 将颜色数输入到数组cc中 end,% 构造自由度为4的卡方分布，并进行检验。 nn(1)=sum(n_pk(1:3); nn(2:4)=n_pk(4:6); nn(5)=sum(n_pk(7:8); m=0.02 0.17 0.4205 0.3195 0.0697*n; chi_4=sum(m-nn).2./m); % 构造统计量 p=chi2cdf(4,chi_4); if p0.95str_pk=pass; elsest

39、r_pk=*; end str_pk,(2) 连（run）检验 把随机数 (4.6.3)按一定规则进行分类，如分为两类，分别记为a, b, 得到形如 aa bbb aaaaaaaa bb 由两类元素a,b组成的序列，我们把位于异类元素之间的同类元素，如abbbba中的bbb类元素称为一个连，连中包含同类元素的个数称为连长，显然连长是一个随机变量。在随机数序列中，出现连长为 i 的连数记为 。总连数记为，构成进行检验的统计量。因此，随机数的连检验是按照随机数出现的先后顺序，重点检验它的连贯现象是否异常的一种方法。,l 正负连检验 把随机数序列 (4.5.3) 变换为 ri 1/2 按正负分为两类

40、，这时a = -1，b= 1 组成正负两类连，根据均匀性，独立性假设，出现 a, b 的概率都是0.5且有 E（l）=N/2+1， D（l ）=（N-1）/ 4， P1=k=2-k（k=1，），则可按（4.5.1），（4.5.2）构造统计量U，和卡方统计量进行检验。 l 升降连检验 把随机序列（4.5.3）按生序或降序规律分为两类，表示随机数的增减及其长度的变化规律，组成升降两类连，这时有 E(l)=(2N4)/3，D(1) = (16N-29)/90,从而我们又可以按(4.5.1)，(4.5.2)构造 U，和卡方统计量进行检验了。,例4.7.6 对1000个随机数进行升降连检验。 解题思路：

41、 1）对1000个随机数按升降连产生连长随机数。 2）分别构造标准正态分布统计量和卡方统计量。 3）对两个统计量进行统计检验。,function striud_u,striud_chi2=run_ud_test(R); % 升降连检验，包括正态和卡方检验 m1=0.625 0.275 0.079167 0.020833; n=length(R);R1=diff(R); AN=0 0 0 0; % 搜索并计算总连长 k=1;j=1;i=1;mk(1)=1; while i=n-2if R1(i)*R1(i+1)0mk(j+1)=i+1;j=j+1;endi=i+1;end,mkmax,mki=m

42、ax(mk); mk=mk(1:mki); mk=diff(mk); nn=hist(mk,5); n1=nn(1:3) sum(nn(4:5); Run_tol=sum(nn);m1=m1*Run_tol; % 计算正态和卡方检验检验 u=(3*Run_tol-2*n+4)/sqrt(1.6*n-2.9); chi_3=sum(n1-m1).2./m1); % 检验 if abs(u)1.96striud_u=pass; elsestriud_u=*; end,5、检验随机序列性质的其它一些统计量 除了上面介绍的方法外，事实上还有其他很多常用的统计量如，最小值MINMUN、即MINIMUN=

43、MINX1, X2, XN。最大值MAXIMUN，即MAXIMUN= MAXX1, X2, XN。极差RANG=MAXIMUN - MINIMUN等。我们可以用上面的方法构造统计量对它们进行检验。,if chi2cdf(3,chi_3) 0.95striud_chi2=pass; elsestriud_chi2=*; end,l 模（MODE） 这是一种中心趋势的度量，类似与均值。如对一个随机变量X抽了六个样6，10，10，4，4，10。 则这个样本的模是10，即出现次数最多的那个数，对一个大样本来讲，如果有N个数都出现同样的最多数，则取其中值为最小的那个。 l 中位数（MEDIAN） 这也是

44、中心趋势的一种度量，将样本X1，X2，XN从小到大排列。记为X（1），X（2），X（N），其最中间的那个数为中位数。,l 偏度SKEWNESS 前面的统计量是反映序列的中心离差趋，而偏度则是衡量X的密度是否偏向一边，即不对称的一种度量，在正态情况下，模，均值，中位数应近乎于重叠于一点，当这三点不重合时，则产生了偏度，当中位数在模的右边时，这时右尾长于左尾，这称为负偏，见图,偏度有时称三阶距，当偏度为0时，曲线对称，当偏度为正时，样本密集于均值左边，当偏度为负时，样本密集于均值右边。偏度公式：,SKWNESS =,(4.6.14),l 峰度KURTOSIS峰度是曲线凸起或平滑的一种度量，正态分布

45、的峰度为0，如果峰度为正表示曲线比真实曲线高出，而偏度为负时，较正态分布平滑，见图4.6.2,峰度有时称为四阶距。公式为：,KURTOSIS =,(4.6.15),八、案3例 随机数检验,用计算机放一组0，1均匀分布样本（N = 1000），进行十二类统计检验，共61个统计量。进行以下检验。1奇序列的14矩检验2偶序列的14矩验检3全序列的14矩验检4混合矩检验5延迟17的相关检验6分布均匀性检验7连检验8随机数的函数检验9顺序统计量检验 10组合规律 11模型计算检验 12子样间的KS检验,制作界面思路： 1）规划界面布局 2）拉出各种统计量的edit界面 3）拉出图形界面 4）拉出push

46、button按纽，注意所有的动作程序都在该按纽子程序下：例如：,function pushbutton1_Callback(hObject, eventdata, handles) % hObject handle to pushbutton1 (see GCBO) % eventdata reserved - to be defined in a future version of MATLAB % handles structure with handles and user data (see GUIDATA)n=str2num(get(handles.edit_num,String);

47、 R=rand(1,n); hist(R); s1,s2,s3,s4=moment_test(R); %u1=sqrt(12*n)*(mean(R)-0.5);set(handles.edit2,String,s1);set(handles.edit3,String,s2);set(handles.edit19,String,s3);set(handles.edit20,String,s4);,% 分布拟合优度检验，一、卡方。二、科尔莫哥若夫斯米尔若夫检验 % 卡方检验 k=8; chi2_str=chi2_test(R,k); set(handles.edit4,String,chi2_str);% 科尔莫哥若夫斯米尔若夫检验 s_ks=ks_test(R,k); set(handles.edit5,String,s_ks);% 组合规律检验 % poke test str_pk=pk_test(R) set(handles.edit13,String,str_pk);,