均匀试验设计的理论、方法和应用)))历史回顾
文章编号:1002)1566(2004)03)0069)1269
均匀试验设计的理论、方法和应用
)))历史回顾
方开泰
(香港浸会大学,数学系)
摘要:本文回顾计算机仿真试验设计的主要两种方法:拉丁超立体抽样和均匀设计,在过去二十五
年的发展,特别是均匀设计的发展,包括均匀设计的优良性研究、新的均匀性测度、均匀设计表的
构造,以及均匀性在因子设计中的应用。
关键词:均匀设计;拉丁超立体抽样;因子设计;正交性;均匀性
中图分类号:O212文献标识码:A
一、历史回顾
廿世纪七十年代,在系统工程、高科技发展的推动下,计算机仿真(仿真)试验(computerexper-iments)的需求十分强烈,迫切要求高质量的试验设计。于是计算机仿真试验设计(Designofcomput-erexperimtnts)在那时成为一个最有挑战的课题。在北美洲,三位学者(McKay,M.D.,Beckman,R.J.andConover,W.J.(1979))在/Technometrics0提出了/拉丁超立方体抽样0(LatinHypercubeSam-pling)(简称LHS)的方法,并立即得到广泛的应用,一批学者对其理论和方法作了系统地研究和发展,形成了一个独立的分枝。差不多在同一时间,在中国,方开泰和王元院士提出了/均匀设计0(UniformDesign)(简称UD)。文章最初在1978年发表在中国科学院数学研究所的内部通讯,后来中、英文稿分别发表在5应用数学学报6和5科学通报6。那时,中国正处于文化大革命刚结束,百废待兴的时代,学术上与世界几乎隔绝。有趣的是,LHS和UD有异曲同工之处。表现于:
(A)两种方法均将试验点均匀地散布于输入参数空间,故在文献中广泛使用术语/充满空间的设计0(spacefillingdesign)LHS给出的试验点带有随机性,故称为抽样;而UD是通过均匀设计表来安排试验,不带有随机性。
(B)两种方法的最初理论均来自/总均值模型0(OverallMeanModel)LHS希望试验点对输出变量的总均值提供一个无偏估值,且方差较小,而UD是希望试验点能给出输出变量总均值离实际总均值的偏差最小。
(C)两种设计均基于U-型设计。
(D)两种设计能应用于多种多样的模型,且对模型的变化有稳健性。
经过了廿多年的发展,两种不同思路的方法是分道扬镳,还是相互补充,相互融合呢?本文想作一些回顾和讨论。为此,我们需要介绍两种方法的思路、模型、方法和应用。
二、总均值模型
设输入变量x1,,,xs与输出变量有一个确定性的关系
y=f(x1,,,xs),x=(x1,,,xs)ICs.
这里假定试验区域为单位立方体Cs=[0,1]s,变量y在Cs上的总均值为sEy)C1,,,x1,,(2.1).2)
70 中文核心期刊 数理统计与管理 23卷3期 2004年5月若在Cs上取了n个试验点,x1,,,xn,y在这n个试验点上的均值为
y (Dn)=f(xi),niE=1
此处Dn={x1,,,xn}代表这n个点的一个设计。n(2.3)
LSH方法是用抽样的方法来选取Dn使相应的估计 y(Dn)是无偏的,即E(y (Dn))=E(y),且方差Var(y (Dn))尽可能地小。McKay,BeckmanandConover(1979)指出LSH比简单随机抽样要好,即前者所获得的总均值比后者有较小的方差。若设计点集Drandom中的点x1,,,xn独立同分布,遵从Cs上的均匀分布,相应样本均值y (Drandom)是E(y)的无偏估计,其方差为Var(f(x))/n,其中x在Cs上均匀分布。如果试验点同分布,但相互之间有相关性,得
Var(DLHS)=Var(f(x))/n+(n-1)Cov(f(x1),f(x2))/n,(2.4)
右边第一项是随机抽样时样本均值的方差,故Var(DLHS)
(1)将试验区域(不失一般性可设为单位正方形),每边均分8份,故单位正方形被分成64=82个小正方形。
(2)随机地取(1,2,,,8)的两个置换,例如(2,5,1,7,4,8,3,6)和(5,8,3,6,1,4,7,8)将它们排成一个矩
2 5 8
1 3
阵,得7 6
4 1
8 4
3 7
6 行、每列有一个小正方形。
(3)在已选中的8个小正方形上分别选点,使小正方形中的每个点被选中的机会相等,于是得图1(b)的结果。这8个点就是一个LHS。。由矩阵的每一行(2,5)、(5,8)、,、(6,2)决定了8个小正方形,如图
1(a),我们看到,每
(a) (b)
图1 LHS(n=8,s=2)
理论上可以证明,LHS的样本均值的方差为
均匀试验设计的理论、方法和应用)))历史回顾
Var(y (DLHS))=(f(x))-+o,c>0,nnn71(2.5)
故LHS比随机抽样的方案好,因为其样本值有较小的方差。
均匀设计也是考虑了/总均值模型0,由数论方法中著名的Koksma-Hlawka不等式,可得,
E(y)-y (Dn)FV(f)D(Dn),(2.6)
式中V(f)是函数f在Cs上的总变差,函数f平稳,V(f)则小,函数f波动大,V(f)则大;D(Dn)为点集在Dn与Cs上的偏差,偏差是度量点集均匀性的一种测度。显然,我们希望D(Dn)越小越好,即点在Cs上散布越均匀越好。仍以n=8,s=2为例
,UD构造的步骤如下:
[i]将试验区域每边均分8份,共得82=64个小正方形。
[ii]考虑在64个小正方形中选取8个小正方形,使得每行每列恰有
一个小正方形。这种取法共有N=8°@8°=403202U
1,625,702,400种,对每一种取法,将试验放在小正方形之中心。然
后比较这幺多设计的偏差(即均匀性),使偏差达到最小的一个即为
均匀设计。
LHS及UD的构造均需要U-型设计的概念。
定义一:有n次试验,s个输入变量的任一个U-型设计可表成一个n
@s的矩阵,第j(j=1,,,s)列的元素由1,2,,,qj组成,且这qj个
数出现的次数相同。该设计记为U(n;q1@,@qs)。若某些qj相图2 均匀设计(n=8,s=2)
m),r1,,等时,该设计为U(n;qr11@,@qr,rm为正整数,且r1+,+rm=s。当q1=,=qs=qm
时,记为U(n;qs)。
在文献中,U-型设计又称为均衡设计(Balanceddesign)或格子点设计(Latticedesign)。
三、LHS的发展
LHS在西方十分流行,有关的文章很多,综合评述的文章不少,如Sack,Welch,MitchellandWynn(1989),Betes,Buck,RiccomangoandWynn(1966),KoehlerandOwen(1996)。熟悉随机抽样的人都知道,随机抽样的表现不稳定,故LHS提供的设计,有的很好,也有的很差。例如图3(a)的设计很好,图3(b)的设计很差。于是需要进一步来改进。其改进的思路列举如下:
(1)中心化的LHS(CenteredLHSorMidpointLHS):研究发展LHS中的步骤(3)并不重要,不如免去在小立方体上抽样,而将试验点放小立方体的中心。以后讲到的LHS均为中心化LHS。
(2)对称LHS(SymmetricLatinHypercube):人们发现,比较好的LHS有某种对称性。图3给出两个LHS(n=6,s=2)。显示(a)的设计比(b)均匀,它们若用矩阵来表示则为
1 1 2 6
D6-1=3 2
4 5
5 1
6 ,and D6-2=2 43 34 25 16 .
注意,D6-1的第一行与第六行之和,第二行与第五行之和,第三行与第四行之和均为(7,7),而D1,
72 中文核心期刊 数理统计与管理 23卷3期 2004年5月Ye,LiandSudjianto(2003)
。
(a) (b)
图3 对称和非对称的LHS
(3)列正交的LHS:如果将设计表的水平安排成零对称,两个列的内积如为零,则这两个称为正交。若一个LHS的任两列都列正交,则有一系列好的性质。表1就是一个列正交的LHS。详见Ye(1998)。
表1 列正交LHSNo
1
2
3
4
5
6
7
8
9112340-4-3-2-12-21-430-34-123-4-3210-1-23443-4-120-214-3(4)将正交设计或均匀设计随机化:这一想法是大大缩小LHS抽样的范围,即只在好的设计中抽样。所谓OA-basedLHS是将正交设计随机化,而(t,m,s)一网的随机化(randomlypermuted(t,m,s)-net)是将均匀设计随机化。这一方法可使均值y (Dn)的方差降至(f(x))-+03/2,d>0.nnn(5)最优LHS(OptimalLHS):改进LHS的另一办法是外加一个准则,用此准则来筛选LHS,求得在此准则下最优的设计。在文献中,流行的准则有IMSE(总均方差)、熵(Entropy)、极小极大(Minimax)、极大极小(Maximin)距离,以及在此基础上定义的Ép
准则等。这一类的方法,是选定一个准则(如上述中的一个),然后在所有U-型设计中找一个使该准则达到最优者,这就是我们追求的OptimalLatinHypercubeDesign。详细讨论可参见KoehlerandOwen(1996)。显然,如果选均匀性作准则,则相应的最优设计就是均匀设计。这时LHS和UD两种不同的思路达到异途同归。美国福特公司的工程经理AgusSudjianto博士正是用这种方法将LHS和UD组合在一起放在他们的软件包中。而且,他们发现,如果均匀性准则取中心化偏差,在上述诸准则中,以均匀设计最好计算,这也是他们近几年来使用均匀设计的原因之一。
四、均匀设计的发展
由于均匀性准则是个几何准则,与统计推断似乎无关,故均匀设计的理论发展比LHS的理论要艰难得多。在过去的十年中,均匀设计的理论与方法有一个飞跃的发展,表现在:
(1)均匀设计的优良性
/总均值模型0虽然在发展LHS及UD起了关键性的作用,但是一个试验的结果,仅仅估计总均值是远远不够的。绝大部份试验的目的,是通过试验数据来建立一个近似模型(有解析表达、好计算 y=f^(1,xs^(),.1)
均匀试验设计的理论、方法和应用)))历史回顾
用它来逼近模型(2.1)。如果,对预先给定的D>0,能使
f(x1,,,xs)-f^(x1,,,xs)
该近似模型就十分理想。在文献中,常常取 f^(x)=B1f1(x)+,+Bmfm(x),73(4.2)其中f1,,,fm为已知函数,B1,,,Bm是欲估的参数。如果能找到这样的f1,,,fm,什么样的设计是最好的呢?谢民育和方开泰(XieandFang(2000))运用统计决策理论的思想和方法,证明了均匀性测度是允许的,minimax的而且在一定的条件下是最优的。
当(4.2)不能很好近似真模型(2.1)时,稳健设计的思想是考虑
y=B1f1(x)+,+Bmfm(x)+h(x)(4.3)
这里h(x)表示拟合不足的部份。稳健设计的思想是假定h(x)来自一个函数类H中,来找一个设计达到最好。众所周知,均方误MSE(Dn,h)= [f(x)-f^(x)]2dx是公认的用来衡量f^(x)近似f(x)质量的标准。由于h(x)的任意性,在文献中常用/平均均方误0和/最大均方误0作为准则,前者是MSE(Dn,h)在hIH中的平均,后者是max[MSE(Dn,h)]。Hickernell(1999)证明,均匀设计在这两个准则下均是最优的。
Wein(1991)也考虑模型(4.3),在试验有重复时,对近似模型(4.1)可作/lackoffittest0,并以(4.3)作备择假设。在一定的函数类H下,他定义了两个准则,并证明均匀设计在两个准则下都是最优的。上述理论的建立使均匀设计建立在更加坚实的基础上,表现均匀设计对建立近似模型f^(x)是最好的。
最近,HickernellandLiu(2002)指出/当用线性模型来拟合数据时,混杂现象不利于模型系数的估计以及确定回归中的项是否显著。如果模型已知,最优设计能给出回归系数的有效估计;稳健试验设计则留心模型的不确定性给回归系数的估计带来的不精确性。虽然很难要求一个设计既是有效的,又是稳健的。但我们将指出均匀设计可限制混杂现象所带来的影响,能提供一个合理的有效及稳健的估计。0
(2)均匀性的度量
王元和方开泰二十五年前在提出均匀设计时,用/星偏差0(stardiscrepancy)作为均匀性的度量,因为星偏差广泛用于数论方法(或伪蒙特卡罗方法)之中,星偏差正是分布似合检验中的柯尔莫哥洛夫)))斯来尔洛夫统计量。由于星偏差不便于计算,我们采用一种近似的算法。Hickernell(1998)指出,星偏差及Lp-偏差对坐标系的旋转缺乏不变性。于是他改用泛函分析的工具)))希尔伯特再生核空间,来定义均匀性测度。
设F(x)为点集N上的均匀分布,Fn(x)是点集Dn={x1,,,xn}的经验分布函数。他们之间的距离可通过泛函分析中的模+F(x)-Fn(x)+来定义。如果用再生核空间,设核函数为对称和非负定的K(x,t),则文献中各种偏差均可表为
D(Dn,K)= K(x,t)d[F(x)-Fn(x)]d[F(t)-Fn(t)].NN
通过取不同的核函数K(x,t),及不同的试验点区域N,可获得不同的偏差。当N为单位立方体时,Hickernell(1998)提出了好几种偏差,其中以中心化L2-偏差(centeredL2-discrepancy),可卷L2-偏差(wrap-aroundL2-discrepancy)具有好的性质:
(a)对坐标系旋转有不变性;
(b)便于计算;
(c
74 中文核心期刊 数理统计与管理 23卷3期 2004年5月
s
当N,,时,可获得散差偏差(discretediscrepancy)。这三种新偏差的提2q2q2q
出,特别是离散偏差的提出,使均匀设计与因子设计的研究相互渗透,相互促进。详细讨论可见Fang,GeLiu(2002b),HickernellandLiu(2002)。
(3)均匀设计表的构造
均匀设计表的构造,有如下的方法:
[1]好格子点法(goodlatticepointmethod):由王元院士和方开泰首先使用,这是数论方法中最经典的方法,由Korobov于1959年提出,至今仍是一种实用而可行的方法。后来,许多学者对此法有很多改进。
[2]拉丁方法:由Fang,Pan,Shiu(1999)提出。
[3]正交表扩充法:由方燕(1995)提出。
[4]折迭法(collapsingmethod):是将两个均匀设计表,用Kronecker乘积的方法折迭在一起。设U和V分别为n@s及m@t的均匀设计表,令 DU,V=(1máUsVá1n)
式中1m为m@1且元素均为1的间量,á表示kronecker乘法。方开泰和覃红(FangandQin(2003))证明用这种方法构造的设计有较好的均匀性。例如 U=
1 2 1 2 1 1 3
则DU,V=1 2 2 2.2 1 2 2
1 2 3 1
2 1 3 [5]切割法(cuttingmethod):用好格子法中的指数生成间量法来生成均匀设计表Up(ps)的计算量很小,当p为素数时,相应的表Up(ps)有很好的均匀性。马长兴和方开泰(MaandFang(2003))提出了所谓的/均割法0,从p个点中取出n(
[6]RBIBD法:可分解的平衡不完全区组设计(resolvablebalancedincompleteblockdesigns)在过去的几十年中有系统的理论和丰富的设计。我和我的合作者们葛根年、刘民千、殷嘉兴、覃红以及陆璇和他的合作者建立了RBIBD和均匀设计之间的联系,通过这个联系可将许多已获得的RBIBD转化为均匀设计表Un(qn)或一些混合水平的表。这些表是在离散偏差下达到最优。有关的研究可类似地用于构造超饱和设计,所谓超饱和设计是指试验的数目少于欲估的参数的数目(主效应、交互效应、误差方差等)。有关的讨论很多,可参见:LuandMeng(2000),LuandSun(2001),Fang,GeandLiu(2002a;b),Fang,GeandLiu(2003),Fang,Ge,LiuandQin(2002,2003),Fang,LinandLiu(2003)等。
[7]翻转法:翻转法(foldover)常用于二水平的试验,例如一个L4(23),通过翻转可以变为L8(23),使得得后者的分辨率(resolution)可以提高一个等级。均匀设计也可通过翻转法来扩大,以及作序贯均匀设计。Ye(1998)给出了一种方法,但更一般地情形尚有待进一步研究。
[8]数值优化法:求均匀设计是一个典型的优化问题。我和P.Winker(WinkerandFang(1998))用(,and1 22 11 3,V=2 23 1
均匀试验设计的理论、方法和应用)))历史回顾75
(2002))用中心化偏差(FangandMa(2001))用可卷偏差,也获得了许多均匀设计表。最近Fang,LuandWinker(2003)对二水平及三水平的均匀表其中心化L2偏差的下界给了一个估计,从而加速了计算的过程。
(4)均匀性在因子设计中的应用
均匀性表面上是一个几何的准则,但它在统计推断中仍有丰富的内涵。在过去的八年中,这些丰富的内涵被逐一地揭示出来。下面将介绍其中的一部份:
[1]均匀性与字长:在正规的因子设计中,字长用来刻划一个因子设计有关主效应和交互效用的混杂情形。若一个s因素二水平设计D,其试验数为n=2s-p的正规试验,它的s个字长组合一个向量W(D)=(A1(D),A2(D),,,As(D))。两个重要的准则/最大分辨力0和/最小低阶混杂0都是W(D)的函数。它们的定义可参见方开泰、马长兴(2001)N2.7。这两个准则似乎与均匀性毫无相关,有趣的是我们(FangandMukerjee(2000))对正规的2s-p的因子设计D,给出了它的中心化L2-偏差和字长间的解析关系
[CD(D)]=212s-232s+9s1+
i=1EsAi(D).9由这个关系不难看出均匀性和最小低阶个混杂准则很相近。随后我们(MaandFang(2001))将上述关系推广至非正规因子设计和可卷L2-偏差。在因子设计的研究中,长期以来人们停留于正规设计,因为非正规因子设计那时没有字长的定义。而均匀性既可用于正规设计,也可用于非正规设计,显示了它的优势。直至2001年,MaandFang(2001),XuandWu(2001)同时将字长由正规因子设计推广至非正规因子设计,但字长的计算比均匀性(如中心化L2-偏差)的计算复杂。均匀性的上述简单性和普遍性显示了它的优点。但均匀性的数值大小有时不易直接给出统计推断的意义。如正交表L9(34)和均匀正交表UL9(34)的前三列的平方中心化L2-偏差分别为0.050059和0.0493645,它们之间均匀性的差别说明了什么?如果仅仅是看它们的数值,只能讲均匀正交表的均匀性比通常的L9(34)好,不能得到更多的性质。我和覃红(FangandQin(2002))提出了一个向量(MI1(D),,,MIs(D)),称为/均匀性指针0(uniformitypattern)。该向量具有字长型W(D)的优点,又有均匀性的通用性。用这个工具,可以方便地比较各种设计。
[2]均匀性与正交性:常见的大多数正交表是强度为2的正交数组,这些表要求一维和二维的投影均匀性,而均匀设计表是要求一维和s维的投影均匀性。其间必有许多联系,也有许多差异。(a)我和Winker在1998年发现,常用的许多正交表是在中心化L2-偏差下的均匀设计,可以通过优化计算来获得,并提出一个猜想:/一切正交表均是一定均匀性测度下的均匀设计0。后来,我们发现这个猜想是不对的,但在一些特殊情形正是对的,Ma,FangandLin(2003)给予了详细的讨论。(b)Fang,MaandMukerjee(2001)对一个部分因子qs-p设计D,定义了一个正交性的向量(B1(D),,,Bs(D)),其中B1(D)表示设计D与强度为1的正交性偏差,B2(D)表示设计D与强度为2的正交性偏差,,。我们在文中给出了中心化L2-偏差,可卷L2-偏差表为(B1(D),,,Bs(D))的解析函数。
(c)均匀性的正交表:如果一个设计既是均匀又是正交设计D,则称为均匀正交设计。我和马长兴在5正交与均匀试验设计6一书中,对L9(34)和UL9(34)进行了比较,指出,UL9(34)能减少混杂。
[3]均匀性与正交表的同构性:用d(n,q,s)表示安排s个q水平的因素并用n次试验的设计。两个d(n,q,s)的设计d1和d2称为同构,如果通过试验次序的交换,因素排列的置换,以及一个或多个因素水平的置换,能将一个设计变为另一个设计,则称为这两个设计同构。要识别两个d(n,q,,°sq,
76 中文核心期刊 数理统计与管理 23卷3期 2004年5月献中快速的算法不多。Ma,FangandLin(2001)利用设计点的均匀性和投影均匀性给出了一个算法,可以有效地鉴别不同构的设计。将这一思想发挥,FangandGe(2003)给出了一个有效算法来识别不等价的哈打马矩阵,并首次发现,36阶的哈达马矩阵至少有392个不等价。
[4]均匀性用于构造超饱和设计:若一个设计的试验点数少于欲估计的参数(主效应、交互效应等)个数,该设计称为超饱和设计。超饱和设计常用来筛选因素。为了构造超饱和设计,文献中提出了许多准则,如E(s2),Ave(V2),Ave(f2),E(fNOD)等,Fang,LinandLiu(2003)指出,在二水平时E(fNOD)与E(s2)等价,在三水平时E(fNOD)与Ave(V2)等价。Fang,Ge,Lin(2002a)指出,在二水平时离散偏差与E(s2)等价,在三水平时,离散偏差与Ave(V2)等价,也与E(d2)(参见LuandSun(2001))等价。而离散偏差可以用于多水平和混合水平的情形,故在上述文献中,利用可分解的平衡不完全区组设计与均匀设计的关系,可以构造出一大批超饱和设计。
五、建模方法(Modelling)
在计算机仿真试验中,欲通过试验数据寻找一个近似于真模型(2.1)的近似模型
y=f^(x1,,,xs)(5.1)
在文献中对f^(x1,,,xs)有好几种术语,如/Approximatemodel0,/metamodel0,,/modelofthemodel0。方开泰、马长兴(2001)在书中详细讨论了用多项式回归模型、人工神经网络、多元样条函数等方法来建模。在实践中,线性模型 y=B1f1(x)+,+Bmfm(x)(5.2)是大多数学者采用的近似模型,式中f1,,,fm已知,B1,,,Bm未知。由于可能组成的f1,,,fm数量很多,回归分析中各种筛选变量的技术十分需要,除了众所周知的/向前法0、/后退法0、/逐步回归法0及/最优子集法0外,由范钊青和李润泽(FangandLi(2001))提出的一种新的惩罚函数(pena-lizedfunction)方法具有很好的性质。李润泽(Li(2002))利用这一方法进行了两个均匀设计案例的研究。
Kriging方法为许多学者所推崇,它是在模型(5.2)的基础上,加了一项平稳高斯过程z(x) y=jEBfj(x1,,,xs)+z(x1,,,xs)j=1m(5.3)
使得在试验点上真值与估计值相等,故这一方法本质上是一个插值的方法。高斯程z(x)中有一些
22未知参数H1,,,Ht,R,我们要同时估计H1,,,Ht,R,B1,,Bm。在软件包MATLAB中已有Kriging
的软件(MATLABBOX)。LiandSudjianto(2003)给出了利用kriging方法建模的均匀设计案例研究,这是美国特汽车公司的一个实验。
Morris,MitchellandYlvisaker(1993)在kriging模型的框架下提出了一个贝叶斯算法,该方法要求y对每个输入变量xj的方向导数,故许多案例是不可行的。
考虑到计算机试验不存在随机误差,许多统计方法直接应用有困难,为了建模,建立一套/Sen-sitivityAnalysis0,其思想是将原模型分解
y=f(x1,,,xs)=f0+
i=1Efi(xi)+Efij(xi,xj)+1=i
希望这种分解可有助于将函数的变化能类似于方差分析中的平方和分解。不失一般性,设试验区域为单位立方体。为此,Sobol(2001)引进了ANOVA分解,他要求函数已中心化,即
1 0ft,t(xt,,,xt)dxt=0, j=1,,,k.1k1kj(5.5),
均匀试验设计的理论、方法和应用)))历史回顾
sf0= [0,1]f(x1,,,xs)dx1,dxs,
s-1fi(xi)= [0,1]f(x1,,,xs)
s-2fij(xi,xj)= [0,1]f(x1,,,xs)77kXiFdxk-f0,kXi,jFdxk-f0-fi(xi)-fj(xj).
这里f0相当于总平均,fi(xi)相当于变量xi的主效应,fij(xi,xj)相当于变量xi和xj的交互效
2ss2应,,。然后,他又引进/方差0的概念 D= [0,1]f(x1,,,xs)dx1,dxs-f0,
s
并证明方差有平方和分解 D=kEEDt1,tk=1t
其中 Dt1,tk= f2t,tdxt1k1,dxtk,
Dkk=1,,,sD
Dij的大小,考核xi和xj的交互作用,用 Sij=.D
Dt1,tk一般地,考核xt1,,xtk的交互作用,用 St1,tk=D.
上述的Sk,Sij,St1,tk称为敏感指针(sensitivityindices),用它可以来分析因素主效应和交互作用的考核变量的主效应,用 Sk=
重要程度,并删去不重要的主效应和交互效应,从而用剩余的项组成一个近似模型。在上述的分析中,要计算许多高维积分。在文献中,发展了一套算法称为FAST(Fourieramplitudesensitivitytest)。该方法,是通过一个Fourier变换,将高维积分变为一维积分,详细讨论可参见SaltelliandBolado(1998)。
由于建构的复杂性,尽管学者们发展了许多方法,但仍然远远不能满足于实际的需要,这是一个很有挑战性的研究方向。
六、应用
均匀设计的应用日益广泛,成功的案例与日俱增,读者不难从各种文献库中发现这些案例。近年来,均匀设计走向国际,有关均匀设计和均匀性的文章在国际刊物上已发表了几十篇,包括国际上顶尖的一些杂志,如/Biometrika0、/Technometrics0、/MathematicsComputation0,SIAM的刊物等。物别令人鼓舞的是,美国福特汽车公司正将均匀设计作为他们推行6Sigma,以及研制新型引擎的常规方法,福特汽车公司的工程经理A.Sudjianto,邀请方开泰去夏和今夏去福特汽车公司讲学,并合作课题,他在邀请信中讲道7Inthepastfewyears,wehavetremendousinusingUniformDesignforcomputerexperiments.Thetechniquehasbecomeacriticalenablerforustoexecute/DesighforSixSig-ma0tosupportnewproductdevelopment,inparticular,automotiveenginedesign.Today,computerexper-imentsusinguniformdesignhavebecomestandardpracticesatFordMotorCompanytosupportearlystageofproductdesignbeforehardwareisavailable.Wewouldliketosharewithyouoursuccessfulrealworldindustrialexperiencesinapplyingthemethodologythatyoudeveloped.Additionally,yourvisitwillbeveryvaluableforustogainmoreinsightaboutthemethodologyaswellastolearnthelatestdevelopmentinthearea.8。福特汽车公司是美国最大的汽车公司,均匀设计在福特汽车公司的运用成功,将会产生巨
78 中文核心期刊 数理统计与管理 23卷3期 2004年5月
[参考文献]
[1] Bates,R.A.,Buck,R.J.,Riccomagno,E.andWynn,H.P.Experimentaldesignandbservationforlargesystems[J].
J.R.Statist.Soc.B.,1996(58).77--94.
[2] Box,G.E.P.andDraper,N.R.Abasisfortheselectionofaresponsesurfacedesign[J].J.Amer.Statist.Assoc.1959
(54).622-654.
[3] Fan,J.andLi,R.VariableSelectionviaNonconcavePenalizedLikelihoodandItsOracleProperties[J].J.Amer.
Statist.Assoc,2001(96).1348-1360.
[4] Fang,K.T.Theuniformdesign:applicationofnumber-theoreticmethodsinexperimentaldesign[J].ActaMath.Ap-
pl.Sinica,1980.363-372.
[5] Fang,K.T.Uniformdesigns[M],EncyclopediaofStatistics,2ndEdition,Wiley,NewYork,toappear.
[6] Fang,K.T.,Ge,G.N.andLiu,M.Q.Uniformsupersaturateddesignanditsconstruction[J].ScienceinChina,Der.
A,2002a(45),1080-1088.
[7] Fang,K.T.,Ge,G.N.andLiu,M.Q.ConstructiononE(fNOD)-optimalsupersaturateddesignsviaroomsquares[J].
inCalcuttaStatisticalAssociationBulletinVol52,A.ChaudhuriandM.GhoshEds.,2002b,71-84.
[8] Fang,K.T.,Ge,G.NandLiu,M.Q.Constructionofoptimalsupersaturateddesignsbythepackingmethod[J].Sc-i
enceinChina,toappear.
[9] Fang,K.t.,Ge,G.N,Liu,M.Q.andQin,H.Constructiononminimumgeneralizedaberrationdesigns[J]Metrika.
2002(57),37-50.
[10] Fang,K.T.andLin,D.K.J.Uniformdesignsandtheirapplicationinindustry,inHandbookonStatistics22:Statistics
inIndustry[M].EdsbyKhattree,R.andRao,C.R.,Elsevier,North-Holland,131-170.
[11] Fang,K.T.,Lin,D.K.J.andLiu,M.Q.Optimalmixed-levelsupersaturateddesign[M].Metrika,toappear.
[12] Fang,K.T.,Lin,D.K.J.andQin,H.Anoteonoptimalfoldoverdesign[J].Statist.&Prob.Letters,2003(62),245-
250.
[13] Fang,K.T.,Lin,D.K.J.,Winker,P.andZhang,Y.Uniformdesign:TheoryandApplications[J].Technometrics,
2000(42),237-248.
[14] Fang,K.T.andMa,C.X.Wrap-aroundL2-discrepancyofrandomsampling,Latinhypercubeanduniformdesigns
[M].J.Complexity,2001,608-624.
[15] Fang,K.T.,Ma,C.X.andWinker,P.CenteredL2-discrepancyofrandomsamplingandLatinhypercubedesign,and
constructionofuniformdesign[J].MathComputation,2002(71)275-296.
[16] Fang,K.T.andMukerjee,R.Aconnectionbetweenuniformityandaberrationinregularfractionsoftwo-levelfactor-i
als[J].Biometrika,2000(87)193-198.
[17] Fang,K.T.andQin,H.Anoteonconstructionofnearlyuniformdesignswithlargenumberofruns[J].Statist.&
Prob.Letters,2002(61)215-224.
[18] Fang,K.T.andWang,Y.Number-theoreticMethodsinStatistics[M].ChapmanandHall,London.1994.
[19] Fang,K.T.andWinker,P.UniformityandOrthogonality,TechnicalReport[M].MATH-175,HongKongBaptistU-
niversity.1998.
[ [R].M,1998a(,299
均匀试验设计的理论、方法和应用)))历史回顾79
[21] Hickernell,F.J.Latticerules:howwelldotheymeasureup?inRandomandQuas-iRandomPointSets[J].EdsP.
HellekalekandG.Larcher,Springer-Verlag,1998b,106-166.
[22] Hickernell,F.J.Goodness-of-FitStatistics,DiscrepanciesandRobustDesigns[J].Statist.Probab.Lett,1999(44),
73--78.
[23] Johnson,M.E.,L.M.Moore,andD.Ylvisaker.Minimaxandmaxmindistancedesigns[J].J.Statist.Plan.andIn-
fer.,1990(26),131-148.
[24] Koehler,J.R.andA.B.Owen.Computerexperiments,inHandbookofStatistics[J].Vol.13,EdsbyGhosh,S.and
Rao,C.R.,ElsevierScienceB.V.,Amsterdam,1996,261-308.
[25] Li,R.Modelselectionforanalysisofuniformdesignandcomputerexperiment[J].InternationalJournalofreliabil-i
ty,QualityandSafetyEngineering,2002(9),305-315.
[26] Liang,Y.Z.,Fang,K.T.andXu,Q.S.Uniformdesignanditsapplicationsinchemistryandchemicalengineering
[J].ChemometricsandIntelligentlaboratorySystems,2001(58),43-57.
[27] Lu,X.andMeng,Y.Anewmethodintheconstructionoftwo-levelsupersaturateddesigns[J].JournalofStatistical
PlanningandInference,2000(86),229-238.
[28] Lu,X.andSun,Y.Supersaturateddesignwithmorethantwolevels[J].ChineseAnnalsofMathematics(SeriesB),
2001(22),183-194.
[29] Ma,C.X.andFang,K.T.Anoteongeneralizedaberrationinfactorialdesigns[J].Metrika,2001(53),85-93.
[30] Ma,C.X.andFang,K.T.Anewapproachtoconstructionofnearlyuniformdesigns[M].InternationalJournalof
MaterialsandProductTechnology,toappear.2003.
[31] Ma,C.X.Fang,K.T.andLin,D.K.J.Onisomorphismoffactorialdesigns[J].J.Complexity,2001(17),86--97.
[32] Ma,C.X.,Fang,K.T.andLin,D.K.J.Anoteonuniformityandorthogonality[J].J.Statist.Plan.Infer.,2003
(113),323-334.
[33] Morris,M.D.,Mitchell,T.J.andYlvisaker,D.Bayesiandesignandanalysisofcomputerexperiments:useofderiva-
tivesinsurfaceprediction[J].Technometrics,1993(35),243-255.
[34] McKay,M.D.,Beckman,R.J.andConover,W.J.Acomparisonofthreemethodsforselectingvaluesofinputvar-i
ablesintheanalysisofoutputfromacomputercode[J].Technometrics,1979(21),239-245.
[35] Owen,A.B.Orthogonalarraysforcomputerexperiments,integrationandvisualization[J].Statist.Sinica,1992(2),
439--452.
[36] Owen,A.B.Latticesamplingrevisited:MonteCarlovarianceofmeansoverrandomizedorthogonalarrays[J].Ann.
Statist,1994(22),930-945.
[37] Owen,A.B.Randomlypermuted(t,m,s)-netsand(t,x)-sequences,InMonteCarloandQuas-iMonteCarloMeth-
odsinScientificComputing,[J].Niederreiter,H.andShiue,P.J.D.,Eds,Springer,NewYork,1995,299-317.
[38] Sacks,J.,Welch,W.J.,Mitchell,T.J.andWynn,H.P.Designandanalysisofcomputerexperiments(withdiscus-
sion)[J].Statist.Sinica,1989(4),409-435.
[39] Saltell,A.andBolado,R.AnalternativewaytocomputeFourieramplitudesensitivitytest[J].Comp.Statist.&Data
Analysis,1998(26),445-460.
[40] Shewry,M.C.andWynn,H.P.Maximumentropysampling[J].J.Appl.Statist.,1987(14),165-170.
[41] Sobol,I.M.GlobalsensitivityindicesfornonlinearmathematicalmodelsandtheirMonteCarloestimates[J].Math.,
80 中文核心期刊 数理统计与管理 23卷3期 2004年5月
[42] Tang,B.Orthogonalarray-basedLatinhypercubes[J].J.Amer.Statist.Asso.,1993(88),1392-1397.
[43] Wang,Y.andFang,K.T.Anoteonuniformdistributionandexperimentaldesign[J].KeXueTongBao,1981(26),
485-489.
[44] Wang,Y.andFang,K.T.Numbertheoreticmethodinappliedstatistics(Ò)[J].ChineseAnnalsofMath.Ser.B,
1990(11),384-394.
[45] Wiens,D.P.Designsforapproximatelylinearregression:twooptimalitypropertiesofuniformdesigns[J].Statist.
Probab.Lett.,1991(12),217-221.
[46] Winker,P.andFang,K.T.OptimalU-typedesign,inMonteCarloandQuas-iMonteCarloMethods1996[M]eds.by
Niederreiter,H.,Zinterhof,P.andHellekalek,P.,Springer,1998,436-448.
[47] Xie,M.Y.andFang,K.T.Admissibilityandminimaxityoftheuniformdesigninnonparametricregressionmodel
[J].J.Statist.Plan.Inference,2000(83),101-111.
[48] Ye,Q.OrthogonalcolumnLatinhypercubesandtheirapplicationincomputerexperiments[J].J.Ammer.Statist.As-
soc.,1998(93),1430-1439.
[49] Ye,Q.andLi,W.andSudjianto,A.AlgorithmicconstructionofoptimalsymmetricLatinhypercubedesigns[J].J.
Statist.Plan.Inference,2003(90),145-159.
[50] Yue,R.X.andHickernell,F.J.Robustdesignsforfittinglinearmodelswithmisspecification[J].Statist.Sinica,
1999(9),1053-1069.
[51] Xu,Q.S.,Liang,Y.Z.andFang,K.T.Theeffectsofdifferentexperimentaldesignsonparameterestimationinthe
kineticsofareversiblechemicalreaction[J].ChemometricsandintelligentLaboratorySystems,2000(52),155-166.
均匀试验设计的理论、方法和应用)))历史回顾
文章编号:1002)1566(2004)03)0069)1269
均匀试验设计的理论、方法和应用
)))历史回顾
方开泰
(香港浸会大学,数学系)
摘要:本文回顾计算机仿真试验设计的主要两种方法:拉丁超立体抽样和均匀设计,在过去二十五
年的发展,特别是均匀设计的发展,包括均匀设计的优良性研究、新的均匀性测度、均匀设计表的
构造,以及均匀性在因子设计中的应用。
关键词:均匀设计;拉丁超立体抽样;因子设计;正交性;均匀性
中图分类号:O212文献标识码:A
一、历史回顾
廿世纪七十年代,在系统工程、高科技发展的推动下,计算机仿真(仿真)试验(computerexper-iments)的需求十分强烈,迫切要求高质量的试验设计。于是计算机仿真试验设计(Designofcomput-erexperimtnts)在那时成为一个最有挑战的课题。在北美洲,三位学者(McKay,M.D.,Beckman,R.J.andConover,W.J.(1979))在/Technometrics0提出了/拉丁超立方体抽样0(LatinHypercubeSam-pling)(简称LHS)的方法,并立即得到广泛的应用,一批学者对其理论和方法作了系统地研究和发展,形成了一个独立的分枝。差不多在同一时间,在中国,方开泰和王元院士提出了/均匀设计0(UniformDesign)(简称UD)。文章最初在1978年发表在中国科学院数学研究所的内部通讯,后来中、英文稿分别发表在5应用数学学报6和5科学通报6。那时,中国正处于文化大革命刚结束,百废待兴的时代,学术上与世界几乎隔绝。有趣的是,LHS和UD有异曲同工之处。表现于:
(A)两种方法均将试验点均匀地散布于输入参数空间,故在文献中广泛使用术语/充满空间的设计0(spacefillingdesign)LHS给出的试验点带有随机性,故称为抽样;而UD是通过均匀设计表来安排试验,不带有随机性。
(B)两种方法的最初理论均来自/总均值模型0(OverallMeanModel)LHS希望试验点对输出变量的总均值提供一个无偏估值,且方差较小,而UD是希望试验点能给出输出变量总均值离实际总均值的偏差最小。
(C)两种设计均基于U-型设计。
(D)两种设计能应用于多种多样的模型,且对模型的变化有稳健性。
经过了廿多年的发展,两种不同思路的方法是分道扬镳,还是相互补充,相互融合呢?本文想作一些回顾和讨论。为此,我们需要介绍两种方法的思路、模型、方法和应用。
二、总均值模型
设输入变量x1,,,xs与输出变量有一个确定性的关系
y=f(x1,,,xs),x=(x1,,,xs)ICs.
这里假定试验区域为单位立方体Cs=[0,1]s,变量y在Cs上的总均值为sEy)C1,,,x1,,(2.1).2)
70 中文核心期刊 数理统计与管理 23卷3期 2004年5月若在Cs上取了n个试验点,x1,,,xn,y在这n个试验点上的均值为
y (Dn)=f(xi),niE=1
此处Dn={x1,,,xn}代表这n个点的一个设计。n(2.3)
LSH方法是用抽样的方法来选取Dn使相应的估计 y(Dn)是无偏的,即E(y (Dn))=E(y),且方差Var(y (Dn))尽可能地小。McKay,BeckmanandConover(1979)指出LSH比简单随机抽样要好,即前者所获得的总均值比后者有较小的方差。若设计点集Drandom中的点x1,,,xn独立同分布,遵从Cs上的均匀分布,相应样本均值y (Drandom)是E(y)的无偏估计,其方差为Var(f(x))/n,其中x在Cs上均匀分布。如果试验点同分布,但相互之间有相关性,得
Var(DLHS)=Var(f(x))/n+(n-1)Cov(f(x1),f(x2))/n,(2.4)
右边第一项是随机抽样时样本均值的方差,故Var(DLHS)
(1)将试验区域(不失一般性可设为单位正方形),每边均分8份,故单位正方形被分成64=82个小正方形。
(2)随机地取(1,2,,,8)的两个置换,例如(2,5,1,7,4,8,3,6)和(5,8,3,6,1,4,7,8)将它们排成一个矩
2 5 8
1 3
阵,得7 6
4 1
8 4
3 7
6 行、每列有一个小正方形。
(3)在已选中的8个小正方形上分别选点,使小正方形中的每个点被选中的机会相等,于是得图1(b)的结果。这8个点就是一个LHS。。由矩阵的每一行(2,5)、(5,8)、,、(6,2)决定了8个小正方形,如图
1(a),我们看到,每
(a) (b)
图1 LHS(n=8,s=2)
理论上可以证明,LHS的样本均值的方差为
均匀试验设计的理论、方法和应用)))历史回顾
Var(y (DLHS))=(f(x))-+o,c>0,nnn71(2.5)
故LHS比随机抽样的方案好,因为其样本值有较小的方差。
均匀设计也是考虑了/总均值模型0,由数论方法中著名的Koksma-Hlawka不等式,可得,
E(y)-y (Dn)FV(f)D(Dn),(2.6)
式中V(f)是函数f在Cs上的总变差,函数f平稳,V(f)则小,函数f波动大,V(f)则大;D(Dn)为点集在Dn与Cs上的偏差,偏差是度量点集均匀性的一种测度。显然,我们希望D(Dn)越小越好,即点在Cs上散布越均匀越好。仍以n=8,s=2为例
,UD构造的步骤如下:
[i]将试验区域每边均分8份,共得82=64个小正方形。
[ii]考虑在64个小正方形中选取8个小正方形,使得每行每列恰有
一个小正方形。这种取法共有N=8°@8°=403202U
1,625,702,400种,对每一种取法,将试验放在小正方形之中心。然
后比较这幺多设计的偏差(即均匀性),使偏差达到最小的一个即为
均匀设计。
LHS及UD的构造均需要U-型设计的概念。
定义一:有n次试验,s个输入变量的任一个U-型设计可表成一个n
@s的矩阵,第j(j=1,,,s)列的元素由1,2,,,qj组成,且这qj个
数出现的次数相同。该设计记为U(n;q1@,@qs)。若某些qj相图2 均匀设计(n=8,s=2)
m),r1,,等时,该设计为U(n;qr11@,@qr,rm为正整数,且r1+,+rm=s。当q1=,=qs=qm
时,记为U(n;qs)。
在文献中,U-型设计又称为均衡设计(Balanceddesign)或格子点设计(Latticedesign)。
三、LHS的发展
LHS在西方十分流行,有关的文章很多,综合评述的文章不少,如Sack,Welch,MitchellandWynn(1989),Betes,Buck,RiccomangoandWynn(1966),KoehlerandOwen(1996)。熟悉随机抽样的人都知道,随机抽样的表现不稳定,故LHS提供的设计,有的很好,也有的很差。例如图3(a)的设计很好,图3(b)的设计很差。于是需要进一步来改进。其改进的思路列举如下:
(1)中心化的LHS(CenteredLHSorMidpointLHS):研究发展LHS中的步骤(3)并不重要,不如免去在小立方体上抽样,而将试验点放小立方体的中心。以后讲到的LHS均为中心化LHS。
(2)对称LHS(SymmetricLatinHypercube):人们发现,比较好的LHS有某种对称性。图3给出两个LHS(n=6,s=2)。显示(a)的设计比(b)均匀,它们若用矩阵来表示则为
1 1 2 6
D6-1=3 2
4 5
5 1
6 ,and D6-2=2 43 34 25 16 .
注意,D6-1的第一行与第六行之和,第二行与第五行之和,第三行与第四行之和均为(7,7),而D1,
72 中文核心期刊 数理统计与管理 23卷3期 2004年5月Ye,LiandSudjianto(2003)
。
(a) (b)
图3 对称和非对称的LHS
(3)列正交的LHS:如果将设计表的水平安排成零对称,两个列的内积如为零,则这两个称为正交。若一个LHS的任两列都列正交,则有一系列好的性质。表1就是一个列正交的LHS。详见Ye(1998)。
表1 列正交LHSNo
1
2
3
4
5
6
7
8
9112340-4-3-2-12-21-430-34-123-4-3210-1-23443-4-120-214-3(4)将正交设计或均匀设计随机化:这一想法是大大缩小LHS抽样的范围,即只在好的设计中抽样。所谓OA-basedLHS是将正交设计随机化,而(t,m,s)一网的随机化(randomlypermuted(t,m,s)-net)是将均匀设计随机化。这一方法可使均值y (Dn)的方差降至(f(x))-+03/2,d>0.nnn(5)最优LHS(OptimalLHS):改进LHS的另一办法是外加一个准则,用此准则来筛选LHS,求得在此准则下最优的设计。在文献中,流行的准则有IMSE(总均方差)、熵(Entropy)、极小极大(Minimax)、极大极小(Maximin)距离,以及在此基础上定义的Ép
准则等。这一类的方法,是选定一个准则(如上述中的一个),然后在所有U-型设计中找一个使该准则达到最优者,这就是我们追求的OptimalLatinHypercubeDesign。详细讨论可参见KoehlerandOwen(1996)。显然,如果选均匀性作准则,则相应的最优设计就是均匀设计。这时LHS和UD两种不同的思路达到异途同归。美国福特公司的工程经理AgusSudjianto博士正是用这种方法将LHS和UD组合在一起放在他们的软件包中。而且,他们发现,如果均匀性准则取中心化偏差,在上述诸准则中,以均匀设计最好计算,这也是他们近几年来使用均匀设计的原因之一。
四、均匀设计的发展
由于均匀性准则是个几何准则,与统计推断似乎无关,故均匀设计的理论发展比LHS的理论要艰难得多。在过去的十年中,均匀设计的理论与方法有一个飞跃的发展,表现在:
(1)均匀设计的优良性
/总均值模型0虽然在发展LHS及UD起了关键性的作用,但是一个试验的结果,仅仅估计总均值是远远不够的。绝大部份试验的目的,是通过试验数据来建立一个近似模型(有解析表达、好计算 y=f^(1,xs^(),.1)
均匀试验设计的理论、方法和应用)))历史回顾
用它来逼近模型(2.1)。如果,对预先给定的D>0,能使
f(x1,,,xs)-f^(x1,,,xs)
该近似模型就十分理想。在文献中,常常取 f^(x)=B1f1(x)+,+Bmfm(x),73(4.2)其中f1,,,fm为已知函数,B1,,,Bm是欲估的参数。如果能找到这样的f1,,,fm,什么样的设计是最好的呢?谢民育和方开泰(XieandFang(2000))运用统计决策理论的思想和方法,证明了均匀性测度是允许的,minimax的而且在一定的条件下是最优的。
当(4.2)不能很好近似真模型(2.1)时,稳健设计的思想是考虑
y=B1f1(x)+,+Bmfm(x)+h(x)(4.3)
这里h(x)表示拟合不足的部份。稳健设计的思想是假定h(x)来自一个函数类H中,来找一个设计达到最好。众所周知,均方误MSE(Dn,h)= [f(x)-f^(x)]2dx是公认的用来衡量f^(x)近似f(x)质量的标准。由于h(x)的任意性,在文献中常用/平均均方误0和/最大均方误0作为准则,前者是MSE(Dn,h)在hIH中的平均,后者是max[MSE(Dn,h)]。Hickernell(1999)证明,均匀设计在这两个准则下均是最优的。
Wein(1991)也考虑模型(4.3),在试验有重复时,对近似模型(4.1)可作/lackoffittest0,并以(4.3)作备择假设。在一定的函数类H下,他定义了两个准则,并证明均匀设计在两个准则下都是最优的。上述理论的建立使均匀设计建立在更加坚实的基础上,表现均匀设计对建立近似模型f^(x)是最好的。
最近,HickernellandLiu(2002)指出/当用线性模型来拟合数据时,混杂现象不利于模型系数的估计以及确定回归中的项是否显著。如果模型已知,最优设计能给出回归系数的有效估计;稳健试验设计则留心模型的不确定性给回归系数的估计带来的不精确性。虽然很难要求一个设计既是有效的,又是稳健的。但我们将指出均匀设计可限制混杂现象所带来的影响,能提供一个合理的有效及稳健的估计。0
(2)均匀性的度量
王元和方开泰二十五年前在提出均匀设计时,用/星偏差0(stardiscrepancy)作为均匀性的度量,因为星偏差广泛用于数论方法(或伪蒙特卡罗方法)之中,星偏差正是分布似合检验中的柯尔莫哥洛夫)))斯来尔洛夫统计量。由于星偏差不便于计算,我们采用一种近似的算法。Hickernell(1998)指出,星偏差及Lp-偏差对坐标系的旋转缺乏不变性。于是他改用泛函分析的工具)))希尔伯特再生核空间,来定义均匀性测度。
设F(x)为点集N上的均匀分布,Fn(x)是点集Dn={x1,,,xn}的经验分布函数。他们之间的距离可通过泛函分析中的模+F(x)-Fn(x)+来定义。如果用再生核空间,设核函数为对称和非负定的K(x,t),则文献中各种偏差均可表为
D(Dn,K)= K(x,t)d[F(x)-Fn(x)]d[F(t)-Fn(t)].NN
通过取不同的核函数K(x,t),及不同的试验点区域N,可获得不同的偏差。当N为单位立方体时,Hickernell(1998)提出了好几种偏差,其中以中心化L2-偏差(centeredL2-discrepancy),可卷L2-偏差(wrap-aroundL2-discrepancy)具有好的性质:
(a)对坐标系旋转有不变性;
(b)便于计算;
(c
74 中文核心期刊 数理统计与管理 23卷3期 2004年5月
s
当N,,时,可获得散差偏差(discretediscrepancy)。这三种新偏差的提2q2q2q
出,特别是离散偏差的提出,使均匀设计与因子设计的研究相互渗透,相互促进。详细讨论可见Fang,GeLiu(2002b),HickernellandLiu(2002)。
(3)均匀设计表的构造
均匀设计表的构造,有如下的方法:
[1]好格子点法(goodlatticepointmethod):由王元院士和方开泰首先使用,这是数论方法中最经典的方法,由Korobov于1959年提出,至今仍是一种实用而可行的方法。后来,许多学者对此法有很多改进。
[2]拉丁方法:由Fang,Pan,Shiu(1999)提出。
[3]正交表扩充法:由方燕(1995)提出。
[4]折迭法(collapsingmethod):是将两个均匀设计表,用Kronecker乘积的方法折迭在一起。设U和V分别为n@s及m@t的均匀设计表,令 DU,V=(1máUsVá1n)
式中1m为m@1且元素均为1的间量,á表示kronecker乘法。方开泰和覃红(FangandQin(2003))证明用这种方法构造的设计有较好的均匀性。例如 U=
1 2 1 2 1 1 3
则DU,V=1 2 2 2.2 1 2 2
1 2 3 1
2 1 3 [5]切割法(cuttingmethod):用好格子法中的指数生成间量法来生成均匀设计表Up(ps)的计算量很小,当p为素数时,相应的表Up(ps)有很好的均匀性。马长兴和方开泰(MaandFang(2003))提出了所谓的/均割法0,从p个点中取出n(
[6]RBIBD法:可分解的平衡不完全区组设计(resolvablebalancedincompleteblockdesigns)在过去的几十年中有系统的理论和丰富的设计。我和我的合作者们葛根年、刘民千、殷嘉兴、覃红以及陆璇和他的合作者建立了RBIBD和均匀设计之间的联系,通过这个联系可将许多已获得的RBIBD转化为均匀设计表Un(qn)或一些混合水平的表。这些表是在离散偏差下达到最优。有关的研究可类似地用于构造超饱和设计,所谓超饱和设计是指试验的数目少于欲估的参数的数目(主效应、交互效应、误差方差等)。有关的讨论很多,可参见:LuandMeng(2000),LuandSun(2001),Fang,GeandLiu(2002a;b),Fang,GeandLiu(2003),Fang,Ge,LiuandQin(2002,2003),Fang,LinandLiu(2003)等。
[7]翻转法:翻转法(foldover)常用于二水平的试验,例如一个L4(23),通过翻转可以变为L8(23),使得得后者的分辨率(resolution)可以提高一个等级。均匀设计也可通过翻转法来扩大,以及作序贯均匀设计。Ye(1998)给出了一种方法,但更一般地情形尚有待进一步研究。
[8]数值优化法:求均匀设计是一个典型的优化问题。我和P.Winker(WinkerandFang(1998))用(,and1 22 11 3,V=2 23 1
均匀试验设计的理论、方法和应用)))历史回顾75
(2002))用中心化偏差(FangandMa(2001))用可卷偏差,也获得了许多均匀设计表。最近Fang,LuandWinker(2003)对二水平及三水平的均匀表其中心化L2偏差的下界给了一个估计,从而加速了计算的过程。
(4)均匀性在因子设计中的应用
均匀性表面上是一个几何的准则,但它在统计推断中仍有丰富的内涵。在过去的八年中,这些丰富的内涵被逐一地揭示出来。下面将介绍其中的一部份:
[1]均匀性与字长:在正规的因子设计中,字长用来刻划一个因子设计有关主效应和交互效用的混杂情形。若一个s因素二水平设计D,其试验数为n=2s-p的正规试验,它的s个字长组合一个向量W(D)=(A1(D),A2(D),,,As(D))。两个重要的准则/最大分辨力0和/最小低阶混杂0都是W(D)的函数。它们的定义可参见方开泰、马长兴(2001)N2.7。这两个准则似乎与均匀性毫无相关,有趣的是我们(FangandMukerjee(2000))对正规的2s-p的因子设计D,给出了它的中心化L2-偏差和字长间的解析关系
[CD(D)]=212s-232s+9s1+
i=1EsAi(D).9由这个关系不难看出均匀性和最小低阶个混杂准则很相近。随后我们(MaandFang(2001))将上述关系推广至非正规因子设计和可卷L2-偏差。在因子设计的研究中,长期以来人们停留于正规设计,因为非正规因子设计那时没有字长的定义。而均匀性既可用于正规设计,也可用于非正规设计,显示了它的优势。直至2001年,MaandFang(2001),XuandWu(2001)同时将字长由正规因子设计推广至非正规因子设计,但字长的计算比均匀性(如中心化L2-偏差)的计算复杂。均匀性的上述简单性和普遍性显示了它的优点。但均匀性的数值大小有时不易直接给出统计推断的意义。如正交表L9(34)和均匀正交表UL9(34)的前三列的平方中心化L2-偏差分别为0.050059和0.0493645,它们之间均匀性的差别说明了什么?如果仅仅是看它们的数值,只能讲均匀正交表的均匀性比通常的L9(34)好,不能得到更多的性质。我和覃红(FangandQin(2002))提出了一个向量(MI1(D),,,MIs(D)),称为/均匀性指针0(uniformitypattern)。该向量具有字长型W(D)的优点,又有均匀性的通用性。用这个工具,可以方便地比较各种设计。
[2]均匀性与正交性:常见的大多数正交表是强度为2的正交数组,这些表要求一维和二维的投影均匀性,而均匀设计表是要求一维和s维的投影均匀性。其间必有许多联系,也有许多差异。(a)我和Winker在1998年发现,常用的许多正交表是在中心化L2-偏差下的均匀设计,可以通过优化计算来获得,并提出一个猜想:/一切正交表均是一定均匀性测度下的均匀设计0。后来,我们发现这个猜想是不对的,但在一些特殊情形正是对的,Ma,FangandLin(2003)给予了详细的讨论。(b)Fang,MaandMukerjee(2001)对一个部分因子qs-p设计D,定义了一个正交性的向量(B1(D),,,Bs(D)),其中B1(D)表示设计D与强度为1的正交性偏差,B2(D)表示设计D与强度为2的正交性偏差,,。我们在文中给出了中心化L2-偏差,可卷L2-偏差表为(B1(D),,,Bs(D))的解析函数。
(c)均匀性的正交表:如果一个设计既是均匀又是正交设计D,则称为均匀正交设计。我和马长兴在5正交与均匀试验设计6一书中,对L9(34)和UL9(34)进行了比较,指出,UL9(34)能减少混杂。
[3]均匀性与正交表的同构性:用d(n,q,s)表示安排s个q水平的因素并用n次试验的设计。两个d(n,q,s)的设计d1和d2称为同构,如果通过试验次序的交换,因素排列的置换,以及一个或多个因素水平的置换,能将一个设计变为另一个设计,则称为这两个设计同构。要识别两个d(n,q,,°sq,
76 中文核心期刊 数理统计与管理 23卷3期 2004年5月献中快速的算法不多。Ma,FangandLin(2001)利用设计点的均匀性和投影均匀性给出了一个算法,可以有效地鉴别不同构的设计。将这一思想发挥,FangandGe(2003)给出了一个有效算法来识别不等价的哈打马矩阵,并首次发现,36阶的哈达马矩阵至少有392个不等价。
[4]均匀性用于构造超饱和设计:若一个设计的试验点数少于欲估计的参数(主效应、交互效应等)个数,该设计称为超饱和设计。超饱和设计常用来筛选因素。为了构造超饱和设计,文献中提出了许多准则,如E(s2),Ave(V2),Ave(f2),E(fNOD)等,Fang,LinandLiu(2003)指出,在二水平时E(fNOD)与E(s2)等价,在三水平时E(fNOD)与Ave(V2)等价。Fang,Ge,Lin(2002a)指出,在二水平时离散偏差与E(s2)等价,在三水平时,离散偏差与Ave(V2)等价,也与E(d2)(参见LuandSun(2001))等价。而离散偏差可以用于多水平和混合水平的情形,故在上述文献中,利用可分解的平衡不完全区组设计与均匀设计的关系,可以构造出一大批超饱和设计。
五、建模方法(Modelling)
在计算机仿真试验中,欲通过试验数据寻找一个近似于真模型(2.1)的近似模型
y=f^(x1,,,xs)(5.1)
在文献中对f^(x1,,,xs)有好几种术语,如/Approximatemodel0,/metamodel0,,/modelofthemodel0。方开泰、马长兴(2001)在书中详细讨论了用多项式回归模型、人工神经网络、多元样条函数等方法来建模。在实践中,线性模型 y=B1f1(x)+,+Bmfm(x)(5.2)是大多数学者采用的近似模型,式中f1,,,fm已知,B1,,,Bm未知。由于可能组成的f1,,,fm数量很多,回归分析中各种筛选变量的技术十分需要,除了众所周知的/向前法0、/后退法0、/逐步回归法0及/最优子集法0外,由范钊青和李润泽(FangandLi(2001))提出的一种新的惩罚函数(pena-lizedfunction)方法具有很好的性质。李润泽(Li(2002))利用这一方法进行了两个均匀设计案例的研究。
Kriging方法为许多学者所推崇,它是在模型(5.2)的基础上,加了一项平稳高斯过程z(x) y=jEBfj(x1,,,xs)+z(x1,,,xs)j=1m(5.3)
使得在试验点上真值与估计值相等,故这一方法本质上是一个插值的方法。高斯程z(x)中有一些
22未知参数H1,,,Ht,R,我们要同时估计H1,,,Ht,R,B1,,Bm。在软件包MATLAB中已有Kriging
的软件(MATLABBOX)。LiandSudjianto(2003)给出了利用kriging方法建模的均匀设计案例研究,这是美国特汽车公司的一个实验。
Morris,MitchellandYlvisaker(1993)在kriging模型的框架下提出了一个贝叶斯算法,该方法要求y对每个输入变量xj的方向导数,故许多案例是不可行的。
考虑到计算机试验不存在随机误差,许多统计方法直接应用有困难,为了建模,建立一套/Sen-sitivityAnalysis0,其思想是将原模型分解
y=f(x1,,,xs)=f0+
i=1Efi(xi)+Efij(xi,xj)+1=i
希望这种分解可有助于将函数的变化能类似于方差分析中的平方和分解。不失一般性,设试验区域为单位立方体。为此,Sobol(2001)引进了ANOVA分解,他要求函数已中心化,即
1 0ft,t(xt,,,xt)dxt=0, j=1,,,k.1k1kj(5.5),
均匀试验设计的理论、方法和应用)))历史回顾
sf0= [0,1]f(x1,,,xs)dx1,dxs,
s-1fi(xi)= [0,1]f(x1,,,xs)
s-2fij(xi,xj)= [0,1]f(x1,,,xs)77kXiFdxk-f0,kXi,jFdxk-f0-fi(xi)-fj(xj).
这里f0相当于总平均,fi(xi)相当于变量xi的主效应,fij(xi,xj)相当于变量xi和xj的交互效
2ss2应,,。然后,他又引进/方差0的概念 D= [0,1]f(x1,,,xs)dx1,dxs-f0,
s
并证明方差有平方和分解 D=kEEDt1,tk=1t
其中 Dt1,tk= f2t,tdxt1k1,dxtk,
Dkk=1,,,sD
Dij的大小,考核xi和xj的交互作用,用 Sij=.D
Dt1,tk一般地,考核xt1,,xtk的交互作用,用 St1,tk=D.
上述的Sk,Sij,St1,tk称为敏感指针(sensitivityindices),用它可以来分析因素主效应和交互作用的考核变量的主效应,用 Sk=
重要程度,并删去不重要的主效应和交互效应,从而用剩余的项组成一个近似模型。在上述的分析中,要计算许多高维积分。在文献中,发展了一套算法称为FAST(Fourieramplitudesensitivitytest)。该方法,是通过一个Fourier变换,将高维积分变为一维积分,详细讨论可参见SaltelliandBolado(1998)。
由于建构的复杂性,尽管学者们发展了许多方法,但仍然远远不能满足于实际的需要,这是一个很有挑战性的研究方向。
六、应用
均匀设计的应用日益广泛,成功的案例与日俱增,读者不难从各种文献库中发现这些案例。近年来,均匀设计走向国际,有关均匀设计和均匀性的文章在国际刊物上已发表了几十篇,包括国际上顶尖的一些杂志,如/Biometrika0、/Technometrics0、/MathematicsComputation0,SIAM的刊物等。物别令人鼓舞的是,美国福特汽车公司正将均匀设计作为他们推行6Sigma,以及研制新型引擎的常规方法,福特汽车公司的工程经理A.Sudjianto,邀请方开泰去夏和今夏去福特汽车公司讲学,并合作课题,他在邀请信中讲道7Inthepastfewyears,wehavetremendousinusingUniformDesignforcomputerexperiments.Thetechniquehasbecomeacriticalenablerforustoexecute/DesighforSixSig-ma0tosupportnewproductdevelopment,inparticular,automotiveenginedesign.Today,computerexper-imentsusinguniformdesignhavebecomestandardpracticesatFordMotorCompanytosupportearlystageofproductdesignbeforehardwareisavailable.Wewouldliketosharewithyouoursuccessfulrealworldindustrialexperiencesinapplyingthemethodologythatyoudeveloped.Additionally,yourvisitwillbeveryvaluableforustogainmoreinsightaboutthemethodologyaswellastolearnthelatestdevelopmentinthearea.8。福特汽车公司是美国最大的汽车公司,均匀设计在福特汽车公司的运用成功,将会产生巨
78 中文核心期刊 数理统计与管理 23卷3期 2004年5月
[参考文献]
[1] Bates,R.A.,Buck,R.J.,Riccomagno,E.andWynn,H.P.Experimentaldesignandbservationforlargesystems[J].
J.R.Statist.Soc.B.,1996(58).77--94.
[2] Box,G.E.P.andDraper,N.R.Abasisfortheselectionofaresponsesurfacedesign[J].J.Amer.Statist.Assoc.1959
(54).622-654.
[3] Fan,J.andLi,R.VariableSelectionviaNonconcavePenalizedLikelihoodandItsOracleProperties[J].J.Amer.
Statist.Assoc,2001(96).1348-1360.
[4] Fang,K.T.Theuniformdesign:applicationofnumber-theoreticmethodsinexperimentaldesign[J].ActaMath.Ap-
pl.Sinica,1980.363-372.
[5] Fang,K.T.Uniformdesigns[M],EncyclopediaofStatistics,2ndEdition,Wiley,NewYork,toappear.
[6] Fang,K.T.,Ge,G.N.andLiu,M.Q.Uniformsupersaturateddesignanditsconstruction[J].ScienceinChina,Der.
A,2002a(45),1080-1088.
[7] Fang,K.T.,Ge,G.N.andLiu,M.Q.ConstructiononE(fNOD)-optimalsupersaturateddesignsviaroomsquares[J].
inCalcuttaStatisticalAssociationBulletinVol52,A.ChaudhuriandM.GhoshEds.,2002b,71-84.
[8] Fang,K.T.,Ge,G.NandLiu,M.Q.Constructionofoptimalsupersaturateddesignsbythepackingmethod[J].Sc-i
enceinChina,toappear.
[9] Fang,K.t.,Ge,G.N,Liu,M.Q.andQin,H.Constructiononminimumgeneralizedaberrationdesigns[J]Metrika.
2002(57),37-50.
[10] Fang,K.T.andLin,D.K.J.Uniformdesignsandtheirapplicationinindustry,inHandbookonStatistics22:Statistics
inIndustry[M].EdsbyKhattree,R.andRao,C.R.,Elsevier,North-Holland,131-170.
[11] Fang,K.T.,Lin,D.K.J.andLiu,M.Q.Optimalmixed-levelsupersaturateddesign[M].Metrika,toappear.
[12] Fang,K.T.,Lin,D.K.J.andQin,H.Anoteonoptimalfoldoverdesign[J].Statist.&Prob.Letters,2003(62),245-
250.
[13] Fang,K.T.,Lin,D.K.J.,Winker,P.andZhang,Y.Uniformdesign:TheoryandApplications[J].Technometrics,
2000(42),237-248.
[14] Fang,K.T.andMa,C.X.Wrap-aroundL2-discrepancyofrandomsampling,Latinhypercubeanduniformdesigns
[M].J.Complexity,2001,608-624.
[15] Fang,K.T.,Ma,C.X.andWinker,P.CenteredL2-discrepancyofrandomsamplingandLatinhypercubedesign,and
constructionofuniformdesign[J].MathComputation,2002(71)275-296.
[16] Fang,K.T.andMukerjee,R.Aconnectionbetweenuniformityandaberrationinregularfractionsoftwo-levelfactor-i
als[J].Biometrika,2000(87)193-198.
[17] Fang,K.T.andQin,H.Anoteonconstructionofnearlyuniformdesignswithlargenumberofruns[J].Statist.&
Prob.Letters,2002(61)215-224.
[18] Fang,K.T.andWang,Y.Number-theoreticMethodsinStatistics[M].ChapmanandHall,London.1994.
[19] Fang,K.T.andWinker,P.UniformityandOrthogonality,TechnicalReport[M].MATH-175,HongKongBaptistU-
niversity.1998.
[ [R].M,1998a(,299
均匀试验设计的理论、方法和应用)))历史回顾79
[21] Hickernell,F.J.Latticerules:howwelldotheymeasureup?inRandomandQuas-iRandomPointSets[J].EdsP.
HellekalekandG.Larcher,Springer-Verlag,1998b,106-166.
[22] Hickernell,F.J.Goodness-of-FitStatistics,DiscrepanciesandRobustDesigns[J].Statist.Probab.Lett,1999(44),
73--78.
[23] Johnson,M.E.,L.M.Moore,andD.Ylvisaker.Minimaxandmaxmindistancedesigns[J].J.Statist.Plan.andIn-
fer.,1990(26),131-148.
[24] Koehler,J.R.andA.B.Owen.Computerexperiments,inHandbookofStatistics[J].Vol.13,EdsbyGhosh,S.and
Rao,C.R.,ElsevierScienceB.V.,Amsterdam,1996,261-308.
[25] Li,R.Modelselectionforanalysisofuniformdesignandcomputerexperiment[J].InternationalJournalofreliabil-i
ty,QualityandSafetyEngineering,2002(9),305-315.
[26] Liang,Y.Z.,Fang,K.T.andXu,Q.S.Uniformdesignanditsapplicationsinchemistryandchemicalengineering
[J].ChemometricsandIntelligentlaboratorySystems,2001(58),43-57.
[27] Lu,X.andMeng,Y.Anewmethodintheconstructionoftwo-levelsupersaturateddesigns[J].JournalofStatistical
PlanningandInference,2000(86),229-238.
[28] Lu,X.andSun,Y.Supersaturateddesignwithmorethantwolevels[J].ChineseAnnalsofMathematics(SeriesB),
2001(22),183-194.
[29] Ma,C.X.andFang,K.T.Anoteongeneralizedaberrationinfactorialdesigns[J].Metrika,2001(53),85-93.
[30] Ma,C.X.andFang,K.T.Anewapproachtoconstructionofnearlyuniformdesigns[M].InternationalJournalof
MaterialsandProductTechnology,toappear.2003.
[31] Ma,C.X.Fang,K.T.andLin,D.K.J.Onisomorphismoffactorialdesigns[J].J.Complexity,2001(17),86--97.
[32] Ma,C.X.,Fang,K.T.andLin,D.K.J.Anoteonuniformityandorthogonality[J].J.Statist.Plan.Infer.,2003
(113),323-334.
[33] Morris,M.D.,Mitchell,T.J.andYlvisaker,D.Bayesiandesignandanalysisofcomputerexperiments:useofderiva-
tivesinsurfaceprediction[J].Technometrics,1993(35),243-255.
[34] McKay,M.D.,Beckman,R.J.andConover,W.J.Acomparisonofthreemethodsforselectingvaluesofinputvar-i
ablesintheanalysisofoutputfromacomputercode[J].Technometrics,1979(21),239-245.
[35] Owen,A.B.Orthogonalarraysforcomputerexperiments,integrationandvisualization[J].Statist.Sinica,1992(2),
439--452.
[36] Owen,A.B.Latticesamplingrevisited:MonteCarlovarianceofmeansoverrandomizedorthogonalarrays[J].Ann.
Statist,1994(22),930-945.
[37] Owen,A.B.Randomlypermuted(t,m,s)-netsand(t,x)-sequences,InMonteCarloandQuas-iMonteCarloMeth-
odsinScientificComputing,[J].Niederreiter,H.andShiue,P.J.D.,Eds,Springer,NewYork,1995,299-317.
[38] Sacks,J.,Welch,W.J.,Mitchell,T.J.andWynn,H.P.Designandanalysisofcomputerexperiments(withdiscus-
sion)[J].Statist.Sinica,1989(4),409-435.
[39] Saltell,A.andBolado,R.AnalternativewaytocomputeFourieramplitudesensitivitytest[J].Comp.Statist.&Data
Analysis,1998(26),445-460.
[40] Shewry,M.C.andWynn,H.P.Maximumentropysampling[J].J.Appl.Statist.,1987(14),165-170.
[41] Sobol,I.M.GlobalsensitivityindicesfornonlinearmathematicalmodelsandtheirMonteCarloestimates[J].Math.,
80 中文核心期刊 数理统计与管理 23卷3期 2004年5月
[42] Tang,B.Orthogonalarray-basedLatinhypercubes[J].J.Amer.Statist.Asso.,1993(88),1392-1397.
[43] Wang,Y.andFang,K.T.Anoteonuniformdistributionandexperimentaldesign[J].KeXueTongBao,1981(26),
485-489.
[44] Wang,Y.andFang,K.T.Numbertheoreticmethodinappliedstatistics(Ò)[J].ChineseAnnalsofMath.Ser.B,
1990(11),384-394.
[45] Wiens,D.P.Designsforapproximatelylinearregression:twooptimalitypropertiesofuniformdesigns[J].Statist.
Probab.Lett.,1991(12),217-221.
[46] Winker,P.andFang,K.T.OptimalU-typedesign,inMonteCarloandQuas-iMonteCarloMethods1996[M]eds.by
Niederreiter,H.,Zinterhof,P.andHellekalek,P.,Springer,1998,436-448.
[47] Xie,M.Y.andFang,K.T.Admissibilityandminimaxityoftheuniformdesigninnonparametricregressionmodel
[J].J.Statist.Plan.Inference,2000(83),101-111.
[48] Ye,Q.OrthogonalcolumnLatinhypercubesandtheirapplicationincomputerexperiments[J].J.Ammer.Statist.As-
soc.,1998(93),1430-1439.
[49] Ye,Q.andLi,W.andSudjianto,A.AlgorithmicconstructionofoptimalsymmetricLatinhypercubedesigns[J].J.
Statist.Plan.Inference,2003(90),145-159.
[50] Yue,R.X.andHickernell,F.J.Robustdesignsforfittinglinearmodelswithmisspecification[J].Statist.Sinica,
1999(9),1053-1069.
[51] Xu,Q.S.,Liang,Y.Z.andFang,K.T.Theeffectsofdifferentexperimentaldesignsonparameterestimationinthe
kineticsofareversiblechemicalreaction[J].ChemometricsandintelligentLaboratorySystems,2000(52),155-166.