幻方在趨勢無關(guān)設(shè)計(jì)中的研究?

馬海南,陳雪平

(1-浙江工業(yè)職業(yè)技術(shù)學(xué)院人文社科部，紹興 312000;2-江蘇理工學(xué)院數(shù)理學(xué)院，常州 213001;3-東南大學(xué)數(shù)學(xué)系，南京 211189)

1 引言

試驗(yàn)設(shè)計(jì)自20世紀(jì)20年代問世至今，已在農(nóng)業(yè)、工業(yè)、生物醫(yī)藥、航空航天、經(jīng)濟(jì)管理等領(lǐng)域得到廣泛應(yīng)用，各種試驗(yàn)設(shè)計(jì)方法也得到了迅速發(fā)展，如區(qū)組設(shè)計(jì)、部分因析設(shè)計(jì)、正交設(shè)計(jì)、均勻設(shè)計(jì)、空間填充設(shè)計(jì)、飽和設(shè)計(jì)、超飽和設(shè)計(jì)、正交拉丁超立方體設(shè)計(jì)等等．一般的試驗(yàn)設(shè)計(jì)方法并不考慮試驗(yàn)順序，因?yàn)樵囼?yàn)順序常采用“隨機(jī)化”準(zhǔn)則．然而，在一些試驗(yàn)中，比如物理或者工程試驗(yàn)中，一次試驗(yàn)常常需要花費(fèi)一定的時(shí)間，而整個(gè)試驗(yàn)做完需要一個(gè)月或更久，那么，在這些需要長時(shí)間完成的試驗(yàn)中，試驗(yàn)單元或者試驗(yàn)環(huán)境有可能發(fā)生了趨勢性的變化．最后的試驗(yàn)結(jié)果易受到時(shí)間或空間效應(yīng)的干擾．于是，許多學(xué)者對設(shè)計(jì)順序進(jìn)行了研究.早期，Phillips[1-3]、Daniel和Wilcoxon[4]、Draper和Stoneman[5]等討論了帶有時(shí)間趨勢效應(yīng)的試驗(yàn)設(shè)計(jì)順序問題，他們將時(shí)間趨勢效應(yīng)T用一個(gè)p階多項(xiàng)式模型來代替，即

其中Ti為時(shí)間或空間趨勢效應(yīng)，i=1,2,…,n為相應(yīng)的試驗(yàn)順序，ai為相應(yīng)的模型參數(shù)．Cheng和Steinberg[6]研究了時(shí)間序列模型下的設(shè)計(jì)順序問題．

本文同樣采用多項(xiàng)式模型來表示趨勢干擾效應(yīng)．在這一個(gè)模型下，Draper和Stoneman[5]討論了與時(shí)間趨勢無關(guān)(trend-free)的二水平設(shè)計(jì)，并且還討論了最小水平變化(level change)設(shè)計(jì)．Cheng和Jacroux[7]給出了二水平時(shí)間趨勢無關(guān)設(shè)計(jì)的一種構(gòu)造方法，Jacroux[8]對帶有時(shí)間趨勢的混合水平設(shè)計(jì)進(jìn)行了研究，Wang[9]借助并列法對特殊的混合水平情形2p4q進(jìn)行了討論．近期，也對含有區(qū)組的帶有時(shí)間趨勢的設(shè)計(jì)順序進(jìn)行了研究[10]．然而，文獻(xiàn)中對于多水平和混合水平情形的研究仍然較少，對高水平和混合水平因析設(shè)計(jì)的研究則更少．

幻方，也稱為魔方(magic square)，我國數(shù)學(xué)家楊輝稱其為“縱橫圖”．其與群論，組合數(shù)學(xué)密切相關(guān)，在純數(shù)學(xué)和應(yīng)用數(shù)學(xué)的研究中都有著重要的意義[11-14]．本文不關(guān)注幻方在數(shù)學(xué)中的性質(zhì)研究或構(gòu)造，主要關(guān)注幻方在統(tǒng)計(jì)中，特別是在工程試驗(yàn)設(shè)計(jì)領(lǐng)域中的應(yīng)用．事實(shí)上，早在1964年，Phillips[1,2]發(fā)現(xiàn)了幻方可以用于帶有時(shí)間趨勢干擾的設(shè)計(jì)，然而這類設(shè)計(jì)僅包含兩個(gè)因子，其水平數(shù)正好等于幻方的階數(shù)．Hedayat[15]也研究了幻方的若干統(tǒng)計(jì)性質(zhì)，指出在對干擾因子的參數(shù)估計(jì)上，幻方可以使得相應(yīng)最小二乘估計(jì)的方差達(dá)到最小．本文將進(jìn)一步說明幻方(純幻方)可以應(yīng)用在兩因子(多因子)試驗(yàn)設(shè)計(jì)中，其對試驗(yàn)因子的參數(shù)估計(jì)與干擾因子是1階趨勢無關(guān)或線性無關(guān)(linear trend free)的，并給出相應(yīng)的構(gòu)造方法．

第2節(jié)主要介紹幻方的基本概念及其在帶有趨勢干擾效應(yīng)的設(shè)計(jì)中的應(yīng)用，借助幻方，可以構(gòu)造出一類與趨勢效應(yīng)無關(guān)的設(shè)計(jì)，且這類設(shè)計(jì)主要包括兩個(gè)多水平因子．第3節(jié)給出主要敘述純幻方在帶有趨勢干擾效應(yīng)的設(shè)計(jì)中的應(yīng)用．借助純幻方，可以構(gòu)造出一類含有多個(gè)與趨勢效應(yīng)無關(guān)因子的設(shè)計(jì)，同時(shí)這類設(shè)計(jì)也可以應(yīng)用在含有交互作用項(xiàng)的設(shè)計(jì)中．第4節(jié)是結(jié)束語．

2 幻方與兩因子趨勢無關(guān)設(shè)計(jì)

定義1[11-13]一個(gè)由1,2,…,n2數(shù)字組成的n行，n列方陣M=(mij)，稱為幻方(magic square)，如果它的每行、列和對角線上的數(shù)字和都相等，即有

表1給出了一個(gè)3階幻方，可以驗(yàn)證，不同行、列和對角線元素之和都相等．即

表1:3階幻方

下面給出趨勢無關(guān)設(shè)計(jì)的相關(guān)性質(zhì)．考慮一個(gè)包含q個(gè)因子，n2次試驗(yàn)的設(shè)計(jì)D，其中，各因子分別取p1,p2,…,pq水平，記 w0(·)=1,w1(·),w2(·),…,ws(·)為 s階正交多項(xiàng)式，特別地，它們定義在等間隔的時(shí)間點(diǎn)上，即有

比如，可以取

于是，在時(shí)間點(diǎn)u的趨勢效應(yīng)可以表示為

其中θ1,θ2,…,θs為相應(yīng)的趨勢項(xiàng)參數(shù)，從而，試驗(yàn)結(jié)果向量Y=(y1,y2,…,yn2)可以表示為一個(gè)可加、同方差線性模型[4,5,7-10,16]

其中β為感興趣參數(shù)，Xd為系數(shù)矩陣，θ =(θ1,θ2,…,θs)′為趨勢效應(yīng)向量，Xt為n2×s階矩陣，其(u,k)元素值為wk(u),1≤u≤n2,1≤k≤s．在模型(1)下，文獻(xiàn)[16]給出了趨勢無關(guān)設(shè)計(jì)所需要滿足的一個(gè)充要條件．

引理1[16]令D為一個(gè)n2行q列的設(shè)計(jì)，其水平數(shù)分別為p1,p2,…,pq．記χd為設(shè)計(jì)D與第i個(gè)因子的關(guān)聯(lián)矩陣，i=1,2,…,q．當(dāng)?shù)趗次試驗(yàn)對應(yīng)的第i個(gè)因子取水平j(luò)i時(shí)，其元素取值為1，否則取值為0，則設(shè)計(jì)D的第i列為s階趨勢無關(guān)，當(dāng)且僅當(dāng)

下面的定理建立了幻方與趨勢無關(guān)設(shè)計(jì)的關(guān)系．即該因子的各個(gè)水平所對應(yīng)的試驗(yàn)順序之和均等于

定理1設(shè)計(jì)D的第i列為1階趨勢無關(guān)的，當(dāng)且僅當(dāng)，這與幻方的性質(zhì)保持一致．

證明 由引理1的充要條件(2)可知

其值即為．

定理2一個(gè)n階幻方可以用來安排含有兩個(gè)n水平因子的設(shè)計(jì)，其兩個(gè)因子均是1階趨勢無關(guān)的．

證明 設(shè)M=(mij)是一個(gè)n階幻方，將幻方的行和列分別看成兩個(gè)n水平因子，分別記為A和B，取值為{0,1,…,n?1}，其中，行數(shù)i減1記做因子A的水平，列數(shù)j減1記做因子B的水平，對應(yīng)(i,j)的幻方元素值mij作為水平組合(i?1,j?1)的試驗(yàn)順序，則從定義1知

由定理1

于是，這兩個(gè)因子A和B均是s=1階趨勢無關(guān)的．

例1考慮一個(gè)化工廠生產(chǎn)的一種化工產(chǎn)品，影響采收率的因素有多個(gè)，包括催化劑種類、反應(yīng)溫度等．并且這類試驗(yàn)往往和反應(yīng)時(shí)間有關(guān)．下面用表1給出的三階幻方，來設(shè)計(jì)一個(gè)與1階時(shí)間趨勢效應(yīng)無關(guān)的設(shè)計(jì)順序，見表2，其包含兩個(gè)三水平因子，取值為{0,1,2}．且

即因子A,B的各個(gè)水平所對應(yīng)的試驗(yàn)順序值之和均為15．其是一個(gè)趨勢無關(guān)設(shè)計(jì)．由表2，第一次試驗(yàn)的水平組合為(2,1)，第二次試驗(yàn)的水平組合為(0,0)，依次類推．

表2: 與時(shí)間趨勢效應(yīng)無關(guān)的采收率試驗(yàn)設(shè)計(jì)(含兩個(gè)因子)

由于試驗(yàn)結(jié)果往往帶有一定的隨機(jī)性，所得試驗(yàn)結(jié)論需要經(jīng)過一定的統(tǒng)計(jì)顯著性檢驗(yàn)．交叉驗(yàn)證可以進(jìn)一步對試驗(yàn)結(jié)果進(jìn)行比較研究，如果一個(gè)設(shè)計(jì)能夠進(jìn)行交叉驗(yàn)證，此時(shí)的結(jié)論具有更高的可靠性．

定義2[11,13]同心幻方也稱嵌套幻方，其可以從中心開始向外輻射，依次剝?nèi)プ钔庖粚雍蟮姆疥嚲腔梅剑?/p>

文獻(xiàn)[13]給出了一個(gè)奇數(shù)階同心幻方的遞推構(gòu)造法，表3為一個(gè)5階同心幻方．其剝?nèi)プ钔庖粚樱梢缘玫揭粋€(gè)3階幻方．

定理3當(dāng)n為奇數(shù)時(shí)，一個(gè)n階同心幻方可以安排含有二個(gè)1階趨勢無關(guān)因子的設(shè)計(jì)，且可以進(jìn)行次交叉驗(yàn)證．

證明 設(shè)M=(mij)是一個(gè)n同心幻方，由定理2，幻方的行和列分別可以看作兩個(gè)n水平因子，記為A,B，它們均是1階趨勢無關(guān)的．所得試驗(yàn)結(jié)果可以進(jìn)行相應(yīng)的顯著性檢驗(yàn)．

然后，將該同心幻方剝?nèi)プ钔庖粚?，此時(shí)可以用來安排一個(gè)含有2個(gè)n?2水平因子的設(shè)計(jì)，且仍是趨勢無關(guān)的．于是，可以采用該部分?jǐn)?shù)據(jù)進(jìn)行相應(yīng)的顯著性檢驗(yàn)，并與第一次結(jié)論進(jìn)行交叉驗(yàn)證．以此類推，將該同心幻方剝?nèi)サ诙?，并進(jìn)行第2次交叉驗(yàn)證．綜上，可以進(jìn)行次交叉驗(yàn)證．

表3:5階同心幻方

3 純幻方與多因子趨勢無關(guān)設(shè)計(jì)

定義3[12,17,18]一個(gè)由1,2,…,n2數(shù)字組成的n行，n列方陣M=(mij)稱為純幻方，或者完美幻方(perfect magic square)、全對稱幻方．如果它的每行、列和泛對角線上的數(shù)字和都相等，即滿足

其中mod(x,y)表示對x取模y運(yùn)算．

表4給出了一個(gè)5階純幻方，容易驗(yàn)證

定理4當(dāng)n為素?cái)?shù)或素?cái)?shù)冪時(shí)，一個(gè)n階純幻方至少可以安排含有四個(gè)n水平因子的因析設(shè)計(jì)，且四個(gè)因子均是1階趨勢無關(guān)的．

證明 設(shè)M=(mij)是一個(gè)n純幻方，由定理2，幻方的行和列分別可以看作兩個(gè)因子，記為A,B，它們均是1階趨勢無關(guān)的．另外，對于純幻方，存在另外兩個(gè)1階趨勢無關(guān)的因子C,D．設(shè)A,B的水平取值分別為a,b，令因子C取值為mod(a+b,n)，因子D取值為mod(a?b,n)，即因子C由幻方中從右上往左下的對角線得到，D由幻方的從左上往右下的對角線得到，則由純幻方定義3可知，所有泛對角線上元素的和均相等．于是，有

由引理1知，因子C,D也是1階趨勢無關(guān)的．即對于任意純幻方，至少可以構(gòu)造含有四個(gè)1階趨勢無關(guān)因子的設(shè)計(jì)順序．

表4:5階純幻方

例2對于采收率試驗(yàn)，例1給出了一個(gè)與時(shí)間干擾效應(yīng)無關(guān)的設(shè)計(jì)，但所包含的試驗(yàn)因子只能為兩個(gè)，如果影響因子多于兩個(gè)，可以考慮下面的設(shè)計(jì)，見表5，其包含四個(gè)因子，水平取值{0,1,2,3,4}．比如，第一次試驗(yàn)的水平組合為(1,2,3,4)，即因子A取1水平，因子B取2水平，因子C取值水平為mod(1+2,5)=3，因子D取值水平為mod(1?2,5)=4．

表5: 與時(shí)間趨勢效應(yīng)無關(guān)的采收率試驗(yàn)設(shè)計(jì)(含四個(gè)因子)

定理5一個(gè)n階純幻方可以用來同時(shí)考察2個(gè)n水平因子及2個(gè)交互效應(yīng)．

證明 設(shè)2個(gè)n水平因子為A,B，相應(yīng)的水平取值為a,b，則取交互效應(yīng)[19]AB的水平取值為mod(a+b,n)和交互效應(yīng)ABn?1的水平取值為

從定義3可知，兩個(gè)交互效應(yīng)均是1階趨勢無關(guān)的，即純幻方也可以用來構(gòu)造交互效應(yīng)是趨勢無關(guān)的設(shè)計(jì)．

當(dāng)n為合數(shù)時(shí)，雖不能像素?cái)?shù)情形直接推出四個(gè)與趨勢無關(guān)的因子，但純幻方仍然可能給出包含多個(gè)因子的設(shè)計(jì)．下面給出一個(gè)包含三個(gè)因子的采收率試驗(yàn)設(shè)計(jì)方案．見表6，因子水平取值為{0,1,2,3}．比如，第一次試驗(yàn)的水平組合為(3,3,0)，即因子A取3水平，因子B取3水平，因子C取0水平．

表6: 與時(shí)間趨勢效應(yīng)無關(guān)的采收率試驗(yàn)設(shè)計(jì)順序(含三個(gè)因子)

4 結(jié)束語

構(gòu)造含有高水平因子的趨勢無關(guān)設(shè)計(jì)是比較困難的一個(gè)問題，本文給出了一種用幻方和純幻方構(gòu)造的簡單方法，所有這些設(shè)計(jì)均可以應(yīng)用在帶有時(shí)間或者空間趨勢干擾效應(yīng)的工程、生物醫(yī)藥等試驗(yàn)中．其可以避免在數(shù)據(jù)分析時(shí)，試驗(yàn)結(jié)果受時(shí)間干擾效應(yīng)的影響，從而造成在參數(shù)估計(jì)[20]和變量選擇上的誤判．由于幻方(純幻方)的特殊性，本文僅討論了幻方(純幻方)在等水平設(shè)計(jì)中的應(yīng)用，幻方(純幻方)是否可以進(jìn)一步應(yīng)用在混合水平趨勢無關(guān)的設(shè)計(jì)中，有待進(jìn)一步的探討．

參考文獻(xiàn)：

[1]Phillips J P N.The use of magic squares for balancing and assessing order effects in some analysis of variance designs[J].Journal of the Royal Statistical Society:Series C,1964,13(2):67-73

[2]Phillips J P N.A simple method of constructing certain magic rectangles of even order[J].The Mathematical Gazette,1968,52:9-12

[3]Phillips J P N.Methods of constructing one-way and factorial designs balanced for trend[J].Journal of the Royal Statistical Society:Series C,1968,17(2):162-170

[4]Daniel C,Wilcoxon F.Fractional 2p?qplans robust against linear and quadratic trends[J].Technometrics,1966,8(2):259-278

[5]Draper N R,Stoneman D M.Factor changes and linear trends in eight-run two-level factorial designs[J].Technometrics,1968,10(2):301-311

[6]Cheng C S,Steinberg D.Trend robust two-level factorial designs[J].Biometrika,1991,78(2):325-336

[7]Cheng C S,Jacroux M.The construction of trend-free run orders of two-level factorial designs[J].Journal of the American Statistical Association,1988,83(404):1152-1158

[8]Jacroux M.On the construction of trend resistant mixed level factorial run orders[J].The Annals of Statistics,1994,22(2):904-916

[9]Wang P C.Level changes and trend resistance in LN(2p4q)orthogonal arrays[J].Statistica Sinica,1996,6(2):471-479

[10]Wang P C,Wu S S.Trend-free designs in blocked factorial experiments[J].Communications in Statistics-Theory and Methods,2013,42(10):1870-1878

[11]Andrews W S.Magic Squares and Cubes(2nd Edition)[M].New York:Dover,1960

[12]徐桂芳，曹敏謙.純幻方的構(gòu)造原理和方法[M].西安：西安交通大學(xué)出版社，1993 Xu G F,Cao M Q.Construction Theory and Method of Pure Magic Squares[M].Xi’an:Xi’an Jiaotong University Press,1993

[13]劉興祥，楊文兵.奇階同心幻方的遞推構(gòu)造法[J].工程數(shù)學(xué)學(xué)報(bào)，2000,17(S1):146-148 Liu X X,Yang W B.Recurrence method in the structure of odd order’s concentric magic square[J].Chinese Journal of Engineering Mathematics,2000,17(S1):146-148

[14]Li W,Wu D H,Pan F C.A construction for doubly pandiagonal magic squares[J].Discrete Mathematics,2012,312(2):479-485

[15]Hedayat A.On a statistical optimality of magic squares[J].Statistics and Probability Letters,1987,5(3):191-192

[16]Dey A,Mukerjee R.Fractional Factorial Plans[M].New York:John Wiley&Sons,1999

[17]徐桂芳.4階純幻方和純幻立方[J].西安交通大學(xué)學(xué)報(bào)，1999,33(2):1-6 Xu G F.Four order of pure magic square and pure magic cube[J].Journal of Xi’an Jiaotong University,1999,33(2):1-6

[18]鄭格于，徐桂芳.5階及6階全對稱幻方[J].上海交通大學(xué)學(xué)報(bào)，2000,34(8):1134-1138 Zheng G Y,Xu G F.Perfect magic squares of 5 and 6 orders[J].Journal of Shanghai Jiao Tong University,2000,34(8):1134-1138

[19]馬海南，陳雪平.直積表中的純凈兩因子交互作用[J].高校應(yīng)用數(shù)學(xué)學(xué)報(bào)，2013,28(3):287-291 Ma H N,Chen X P.Clear two-factor interactions in cross arrays[J].Applied Mathematics:A Journal of Chinese Universities,2013,28(3):287-291

[20]戴維紀(jì)，王慶祥.用波動(dòng)因素試驗(yàn)法進(jìn)行參數(shù)設(shè)計(jì)試驗(yàn)[J].工程數(shù)學(xué)學(xué)報(bào)，1998,15(2):139-142 Dai W J,Wang Q X.Making experiments of designing parameters with fluctuation factor method[J].Chinese Journal of Engineering Mathematics,1998,15(2):139-142