第三节均匀设计表的构造和运用本节介绍均匀设计表的构造和使用表的来源,其中均匀性度量──偏差将起关键作用,我们将介绍偏差的定义,并给出正交设计与均匀设计各自偏差的比较,从中可以了解为什么均匀设计可以比正交设计节省试验次数,本节还介绍拟水平在均匀设计中的使用和有关表的构造,熟悉本节内容对于正确理解和使用均匀设计有很大帮助。3.1均匀设计表的构造定义1每一个均匀设计表是一个方阵,设方阵有n行m列,每一行是{1,2,...,n}的一个置换(即1,2,…,n的重新排列),表的第一行是{1,2,…,n}的一个子集,但不一定是真子集。显然,第一章表4-6列举的U6*(64),U7(74)和U7*(74)都符合上述定义。符合定义1的均匀设计表数量太多,本节仅介绍用好格子点法(goodlatticepoint)构造的均匀设计表,其方法如下:1)给定试验数n,寻找比n小的整数h,且使n和h的最大公约数为1。符合这些条件的正整数组成一个向量h=(h1,…,hm)。2)均匀设计表的第j列下法生成iijihu[modn](3.1)这里[modn]表示同余运算,若jhi超过n,则用它减去n的一个适当倍数,使差落在[1,n]之中。Uij可以递推来生成uhjj1nhuhuujijjijji,1nhunhujijjij若若1,,1ni(3.2)例如,当n=9时,符合条件1)的h有1,2,4,5,7,8;而h=3或h=6时不符合条件1),因为最大公约数(3,9)=3,(6,9)=3,均大于1.所以9U最多只可能有6列,又如当h34时,用公式(3.2)来生成该列时其结果依次如下:uuuuuuuuu132333435363738393444884123934774112924664101914559,,(mod),(mod),(mod),其结果列于表16的第三列。表16U969()123456112457822481573363636448721555127846636363775184288754219999999用上述步骤生成的均匀设计表记作Unnm(),向量h称为该表的生成向量,有时为了强调h的作用,可将Unnm()记成Uhn().给定n,相应的h可以象上例那样方便地求得,从而m也就确定.所以m是n的一个函数,这个函数曾由大数学家欧拉研究过,称为欧拉函数,记为E(n).这个函数告诉我们均匀设计表最多可能有多少列.下面的结果来自数论:i)当n为素数时,E(n-1)=n-1所谓素数就是一个正整数,它与其所有比它小的正整数的最大公约数均为1.如2,3,4,5,11,13,…均为素数。ii)当n为素数幂时,即n可表成n=pl,这里p为素数l,l为正整数,这时Ennp()()11(3.3)例如n=9可表为n32,于是E()()991136即U9至多可以有6列。iii)若n不属于上述两种情形,这时n一定可以表为不同素数的方幂积,即npppllsls1212(3.4)这里pps1,,为不同的素数,lls1,,为正整数,这时Ennp()()11…()11ps(3.5)例如n=12可表为n=232,于是E()()()12121121134即U12最多只可能有4列。上述三种情形中,以素数情形为最好,我们最多可以获得n-1列,而非素数情形,在上述表的结构中永远不可能有n-1列,例如n=6=2311,此时E()()()661121132,这说明,当n=6时,用上述办法生成的均匀设计表只有2列,即最多只能安排两个因素,这是太少了,为此,王元,方开泰(1981)建议,可将U767()表的最后一行去掉来构造U6,为了区别于由(3.2)生成的均匀设计表,我们记它为U666*(),在U的右上角加一个“*”号,表U666*()列于表17,对照表16我们看到U表和U*表之间的关系和各自特点:i)所有的Un*表是由Un1表中划去最后一行而获得;ii)Un表的最后一行全部由水平n组成,Un*表的最后一行则不然。若每个因素的水平都是由低到高排列,Un表中最后一号表17)6(6*6UNo.123456112345622461353362514441526355316426654321试验将是所有最高水平相组合,在有些试验中,例如在化工试验中,所有最高水平组合在一起可能使反应过分剧烈,甚至爆炸。反之,若每个因素的水平都是由高到低排列,则Un表中最后一号试验将是所有低水平的组合,有时也会出现反常现象,甚至化学反应不能进行。Un*表则没有类似现象,比较容易安排试验。iii)若n为偶数,Un*表比Un表有更多的列。如上面讨论过的U6表只有2列,而Un*表可以有6列。iv)若n为奇数,则Un*表列数通常少于Un表。v)Un*表比Un表有更好的均匀性,应优先采用Un*表,其细节将在下节讨论。vi)若将Un或Un*的元素组成一个矩阵的秩最多分别为En()12及En()112。3.2均匀性准则和使用表的产生在1.6节曾指出均匀设计在使用时由于选择的列不同,试验的效果也大不相同,于是建议读者按使用表的推荐去选列,那么使用表又是如何产生的呢?设我们要从均匀设计表Unnm()中选出s列,则可能的选择有()sm种可能,我们要从中选择一个最好的,这里必须对“好”和“坏”有明确的含义,表Unnm()是由它的生成向量hhhm(,,)1所唯一确定的,选择s列,本质上就是从h中选择s个hhm1,,,由这s个数生成的均匀设计表为Uhhnii(,,),这是一个n×s矩阵。它的每一行是s维空间Rs中的一个点,故n行对应Rs中的n个点,若这n个点在试验范围内均匀,则试验效果好,否则试验效果不好。因此,比较两个均匀设计表Uhhniis(,,)1和Uhhnjjs(,,)1的好坏等价于比较由它们所对应的两组点集的均匀性。于是我们必须给出均匀性度量。度量均匀性准则很多,其中偏差(discrepancy)是使用历史最久,为公众所广泛接受的准则,我们先给出它的定义。设Unnm()是一个均匀设计表,若把它的每一行看成m维空间的一个点,则Unnm()给出了n个试验点,这些点的坐标由{1,2,…,n}组成,用线性变换将{1,…,n}均匀地变到(0,1)之间如下:iinin21212,,,,若用qki表示Unnm()中的元素,则上面的变换等价于令n,,1k),x,,x(xn,1k,m,,1i,n21qki2xkm1kkki(3.6)于是n个试验点变换成[,]01mmC中的n个点:xxn1,,.考虑原n个试验点的均匀性,等价于考核xxn1,,在Cm的均匀性。定义2设xxn1,,为Cm中的n个点,任一向量xxxCmm(,,)1,记vxxxm()1为矩形[0,x]的体积,nx为xxn1中落入[0,x]的点数,则xvnnsup)x,,x(DxCxn1m(3.7)称为点集{,,}xxn1在Cm中的偏差(discrepancy)。为什么偏差可以用于度量点集散布的均匀性呢?若n个点xxn1,,在Cm中散布均匀,则nnx/表示有多少比例的点落在矩形[0,x]中,它应当和该矩形的体积v(x)相差不会太远。如果用统计学的语言来解释偏差,令n1kkn}xx{In1)x(F(3.8)表示的{,,}xxn1经验分布函数,式中I{.}为示性函数,令F(x)为Cm上均匀分布的分布函数,于是(3.7)定义的偏差可表为xFxFsup)x,,x(DnRxn1m(3.9)偏差实际上就是在分布拟合检验中的Kolmogorov-Smirnov统计量,它给出了经验和理论分布之间的偏差。在Cm中任给n个点xxn1,,,如何计算它们的偏差对均匀设计表的构造十分重要。长期以来,一直没有人编出一个实用的算法。当方开泰在1978年提出均匀设计时,只好把偏差展开成级数,取其首项,给出近似偏差的准则。此方法方便计算,但有时有大的偏差,而且只适用于格子点法构造的均匀设计,不能计算正交设计等其它方法所产生试验点的偏差,最近Bundschuh和Zhu(朱尧辰)[17]给出了计算偏差的算法,当因素数不太多时,他们的算法可以精确地求出任何点集的偏差,他们已用MATLAB编出有关的程序,本讲议中的计算,都是用该程序获得的。设我们要从均匀设计表Unnm()中选出s列,使其相应的均匀设计有最小的偏差.当m和s较大时,由m列中取出s列的数目有()sm之多,要比较这么多组点集的均匀性工作量很大.于是需要有简化计算和近似求解的方法.详细讨论可参看方开泰[2],方开泰、郑胡灵[12]等.这里仅仅介绍利用整数的同余幂来产生hhiis1,,的办法。令a为小于n的整数,且a,a2(modn),…,at(modn)互不相同,at+1=1(modn),则称a对n的次数为t,例如222423211234,,,(mod5)则2对5的次数为3.又如333935343112345,,,,(mod9)表示3对9的次数为4.一般若a对n的次数大于或等于s-1,且(a,n)=1,则可用(,,,)aaas01(modn)(3.10)作为生成向量,故a称为均匀设计的生成元.然后在一切可能的a(最多n-1个)中去比较相应试验点的均匀性,工作量则大大减少.理论和实践证明,这种方法获得的均匀设计使用表仍能保证设计的均匀性.于是,给定n和s,只要求得最优的a,便可获得生成向量,从而获得相应的均匀设计表。表18对奇数n(5≤n≤31,n=37)给出了Un表的生成元及其相应均匀设计的偏差.同时对偶数n(6≤n≤30)给出了Un*表的生成元和相应的偏差.类似地,对奇数n,我们也获得Un*表的生成向量和相应均匀设计表的偏差(表19).表18和19的结果取自FangandLi[14].综合两个表的结果,我们有如下的说明。i)对奇数n,Un*表比Un表有更好的均匀性,例如n=15,s=4时,U15(154)的偏差为D=0.2772,而U15415*()的偏差为D=0.1511,后者比前者相对降低了0277201511027724549%....表19中p%一列给出了所有情形偏差降低的百分比.为了直观起见,我们将表18和表19的偏差点成图11.我们按s=2,3,4,5分成四个图.图中“+”表示奇数n的Un表的偏差,“*”表示偶数表18Un和Un*的生成元和相应设计的偏差n23456752(.3100)2(.4570)63(.1875)3(.2656)3(.2990)73(.2398)3(.3721)3(.4760)84(.1445)4(.2000)2(.2709)94(.1944)4(.3102)2(.4066)107(.1125)7(.1681)5(.2236)5(.2414)7(.2994)117(.1634)7(.2649)7(.3528)7(.4286)7(.4942)125(.1163)6(.1838)6(.2233)4(.2272)6(.2670)6(.2768)135(.1405)6(.2308)6(.3107)6(.3814)6(.4439)6(.4992)1411(.0957)7(.1455)7(.2091)1511(.1233)7(.2043)7(.2772)1610(.0908)5(.1262)5(.1705)5(.2070)10(.2518)2(.2769)1711(.1099)10(.1832)10(.2501)10(.3111)10(.3667)10(.4174)188(.0779)9(.1394)9(.1754)4(.2047)3(.2245)9(.2247)198(.0990)8(.1660)14(.2277)14(.2845)14(.3368)14(.3850)2013(.0947)5(.1363)10(.1915)10(.2012)10(.2010)2113(.094