靈活的多層嵌套拉丁超立方體設計構(gòu)造

徐 琎，段曉君，王正明，晏 良

(1. 國防科技大學 文理學院， 湖南 長沙 410073； 2. 國防科技大學 前沿交叉學科學院， 湖南 長沙 410073)

幾乎所有的計算機試驗都存在多種精度的情況，如何將這多種試驗結(jié)合起來并充分利用所有數(shù)據(jù)信息是十分關鍵的問題。多種精度的試驗有許多的應用，如飛行器設計相關試驗等[6]。基于多種精度試驗的建模方法也層出不窮，如貝葉斯修正方法[7]、貝葉斯分層建模[8]和計算機試驗與實體試驗的聯(lián)合建模[9]等。在多種精度試驗的設計選點部分，文獻[10]提出了一種方法：嵌套拉丁方設計(Nested Latin Hypercube Designs, NLHD)。該設計保證了所有高精度的試驗點在低精度試驗中進行了同樣的試驗，這對后續(xù)試驗的建模有很好的幫助[11]，同時不同精度下的試驗點集均為LHD。然而，SLHD及其若干后續(xù)研究(例如文獻[12]和文獻[13])均限制低精度試驗的試驗次數(shù)必須是高精度的倍數(shù)。這就導致在具體的試驗中，試驗次數(shù)會有諸多限制。文獻[14]和文獻[15]介紹了結(jié)構(gòu)靈活的設計，分別為序貫設計(Sequential Designs, SD)和廣義嵌套拉丁超立方體設計(General Nested Latin Hypercube designs, GNLH)，但SD方法的均勻性不夠，其較低精度的試驗點在大多數(shù)情況下不是LHD，而GNLH只能構(gòu)造兩層的嵌套設計。因此，對于三種或三種以上精度的計算機試驗，現(xiàn)有工作無法給出結(jié)構(gòu)靈活的嵌套拉丁超立方體設計。

本文介紹一種結(jié)構(gòu)靈活的多層嵌套拉丁超立方體設計(Multi-layer Nested Latin hypercube Designs, MND)，其試驗設計點數(shù)可以靈活選取，同時每類精度試驗的試驗點均有很好的一維投影均勻性。

1 MND的結(jié)構(gòu)與性質(zhì)

1.1 MND的構(gòu)造方法

下面對MND構(gòu)造算法中使用到的符號進行說明。Zn表示整數(shù)集合{1,…,n}。對于向量v，v(i)表示其第i個元素。對于矩陣A，A(i,j)為第i行j列元素。

對于給定的整數(shù)u≥3,0m，Mm,n為

對于j=1,…,u-1，vj為lj維隨機向量，其中vj(1)服從Zlj+1/lj上的離散均勻分布；對于i=2,…,lj，vj(i)=Mlj+1/mj+1,lj+1/lj(ri,vj(i-1))，其中ri服從Z(lj+1/lj)-(lj+1/mj+1)+1上的離散均勻分布。

容易看出，該算法有效需滿足的條件是：對于i=1,…,u-1，mi+1是mi的倍數(shù)或mi+1≥li。由于實際應用中，高精度與低精度試驗的試驗點一般差距較大，故大多數(shù)的試驗可以滿足上述兩個條件。同時，對于j=1,…,u-1，根據(jù)lj維向量vj的定義可知，其元素均屬于Zlj+1/lj。再由τ1的元素均屬于Zm1可以得出τj+1的元素均屬于Zlj+1，因

算法1 構(gòu)造MND

此，算法1步驟8中元素vj(τj(k))對任意的j=1,…,u-1均存在。下面用一個例子來更好地說明該算法的流程。

例1對于u=3,m1=2,m2=5,m3=13,p=1，隨機生成τ1=(2,1)T，且v1來自矩陣M2,5。隨機構(gòu)造后得向量v1=(4,2)T，則

τ2(1)=5(τ1(1)-1)+v1(τ1(1))=7

同理，τ2(2)=4。由算法1步驟10隨機生成τ2剩余的3個元素，得τ2=(7,4,9,2,6)T。v2來自矩陣M10,13。隨機構(gòu)造后得向量v2=(1,13,10,8,6,4,2,2,12,9)T。則

τ3(1)=13(τ2(1)-1)+v2(τ2(1))=80

同理，得到

τ2=(80,47,116,26,69,124,31,95,15,90,105,60,3)T

則嵌套向量d=(τ3-ε)/130。圖1顯示了嵌套向量d，即MND的一維投影均勻性。由圖可知，對任一層，其均可滿足在任一區(qū)間中有且僅有一個點。

圖1 MND的一維均勻性Fig.1 One-dimensional uniformity of MND

1.2 MND的統(tǒng)計性質(zhì)

進一步，對于MND的性質(zhì)有如下定理。

定理1對于給定的整數(shù)u≥3，0

1)對于任意的k=1,…,u及i,j=1,…,mu且i≠j，有

「τ(i)mk/l?≠「τ(j)mk/l?

2)對于i=1,…,mu，a=1,…l，有

P{τu(i)=a}=1/(mul)

證明：類似文獻[15]中命題1的證明，易知該定理成立。

□

關于MND設計D=(d(1),…,d(p))有如下兩個定理。

定理2對于i=1,…,u,j=1,…,p，{d(j)(1),…,d(j)(mi)}滿足在任一區(qū)間[0,1/mi),…,[(mi-1)/mi,1)中有且僅有一個設計點。

證明：由定理1 第1部分可知。

□

定理3對于i=1,…,mu，{d(1)(i),…,d(p)(i)}服從[0,1)p上的均勻分布。

證明：由定理1 第2部分可知。

□

2 算例

對于運行次數(shù)分別為m1,…,mu的多種精度試驗，有四種方法可對其進行處理，依次為：

1)獨立同分布抽樣(Independent and Identically Distributed sample, IID)：選取mu行的獨立同分布抽樣，利用其前mi行處理第i個試驗。

2)LHD：選取mu行的拉丁超立方體設計矩陣，利用其前mi行處理第i個試驗。

3)文獻[14]的方法SD。

4)多層嵌套設計MND。

1)M1：

2)M2：

其中，M1中的函數(shù)h1(x)由文獻[16]提出，M1中h2(x)和h3(x)來自文獻[17]。M2中h1(x)被稱為borehole函數(shù)[18]，在計算機試驗中經(jīng)常選其作為黑箱函數(shù)[13,16,19]。在M2中代替低精度函數(shù)的h2(x)和h3(x)由文獻[13]構(gòu)造。對于不同方法生成的設計矩陣D，其前m1試驗點對應模型h1(x)，前m2試驗點對應模型h2(x)，整個設計矩陣對應模型h3(x)。圖2和圖3分別為在M1和M2下不同方法估計μi的均方誤差。其中選定的參數(shù)為m1=2，m2=5，m3=13，每類試驗進行2000次，得出各方法的RMSE后進行對比。從兩圖中可以得到：①IID方法的均方誤差最大；②在高精度試驗中，SD和MND均比LHD表現(xiàn)好，其原因在于這兩者的子矩陣比LHD的子矩陣有更好的均勻性；③在估計μ1時， MND與SD有幾乎相同的均方誤差；④在估計μ2和μ3時，MND比SD有更小的均方誤差，原因在于MND比SD在整體上有更好的均勻性；⑤在估計μ3時，MND和LHD有相同的RMSE，這是因為MND和LHD的全矩陣有相同的均勻性。

(a) M1中模型h1(x)(a) Model h1(x) in M1

(b) M1中模型h2(x)(b) Model h2(x) in M1

(c) M1中模型h3(x)(c) Model h3(x) in M1圖2 M1下各方法RMSE比較Fig.2 RMSE comparison of the different schemes on M1

3 結(jié)論

本文構(gòu)造了靈活的多層嵌套拉丁超立方體設計，與一般的嵌套設計相比，該方法在參數(shù)選取上具有很強的靈活性。與SD相比，該方法各層在一維上均可達到更好的均勻性，這使得在估計模型均值時，其估計方差更小。兩個算例也體現(xiàn)了多層嵌套拉丁超立方體設計的性質(zhì)。雖然第1節(jié)中給出的算法也有其局限性，即對于i=1,…,u-1，mi+1是mi的倍數(shù)或mi+1≥li需滿足，但由于實際應用中，高精度與低精度試驗的試驗點一般差距較大，故大多數(shù)的試驗可以滿足上述兩個條件。

(a) M2中模型h1(x)(a) Model h1(x) in M2

(b) M2中模型h2(x)(b) Model h2(x) in M2

(c) M2中模型h3(x)(c) Model h3(x) in M2圖3 M2下各方法RMSE比較Fig.3 RMSE comparison of the different schemes on M2