常用正交表的構(gòu)造原理及SAS實(shí)現(xiàn)

華中科技大學(xué)同濟(jì)醫(yī)學(xué)院公共衛(wèi)生學(xué)院流行病與衛(wèi)生統(tǒng)計(jì)學(xué)系(430030)

陳遠(yuǎn)方 林曦晨 徐利華 湯洪秀 汪宏晶 尹 平△

正交試驗(yàn)是通過一套規(guī)格化的正交表將各試驗(yàn)因素、各水平進(jìn)行均勻搭配，只需較少的樣本量即可找出最優(yōu)組合，具有試驗(yàn)次數(shù)少、效率高、設(shè)計(jì)簡便等特點(diǎn)，已被廣泛應(yīng)用于醫(yī)學(xué)等研究領(lǐng)域〔1－3〕。但正交表種類繁多，各種類型的構(gòu)造原理不同且復(fù)雜，特別是混合正交表，操作過程很是繁瑣〔4－5〕。同時(shí)，目前缺乏采用通用軟件生成正交表的文獻(xiàn)介紹。為此，本文詳細(xì)分析了正交表的構(gòu)造原理并進(jìn)行了分類，采用SAS軟件構(gòu)建了常用正交表的程序模塊，極大地降低了正交試驗(yàn)中建立正交表的難度，方便實(shí)用。

方法和原理

1．2水平正交表

(1)哈達(dá)瑪矩陣直積法構(gòu)造N=2s型

N=2s型的2水平正交表任意兩列都是正交的，而且每一列都與全1列正交。假設(shè)A，B分別為n×m與s×t矩陣，將A的每一個(gè)元素都乘以矩陣B而得到的ns×mt矩陣，用A?B表示，叫做矩陣A與B的直積。當(dāng)A與B都是元素為±1的哈達(dá)瑪矩陣時(shí)，由直積的意義和運(yùn)算規(guī)則可知:

則A與B的直積也是哈達(dá)瑪矩陣。應(yīng)用上述構(gòu)造標(biāo)準(zhǔn)哈達(dá)瑪矩陣中，除全1列外，其任意兩列與全1列都正交。因此，去掉全1列，就得到N=2s型的2水平正交表〔6〕。

(2)哈達(dá)瑪矩陣特征函數(shù)法構(gòu)造N≠2s且N=2(p+1)(p為素?cái)?shù)且p≡3(mod4))型

特征值函數(shù)法利用的是構(gòu)造2(p+1)階哈達(dá)瑪矩陣的思想，建立有限域 GF(p)={0，P-1，P-2，…，1}，利用有限域GF(p)的元根，構(gòu)造R矩陣，然后構(gòu)造K矩陣，再構(gòu)造出2(p+1)階方陣為哈達(dá)瑪矩陣，最終構(gòu)造出2(p+1)階正交表。

當(dāng)N≠2s且 N=2(p+1)型(p為素?cái)?shù)且 p≡1(mod4))時(shí)，令N=2(p+1)，構(gòu)造pn階有限域GF(p)={0，P-1，P-2，…，1}，并給出 GF(p)的一個(gè)元根。把有限域GF(p)中每一個(gè)非零元素都表示為元根的方冪，由特征函數(shù)求出所有的

=2(p+1)E2(P+1)

即2(p+1)階方陣為哈達(dá)瑪矩陣。

在上面構(gòu)造的2(p+1)階哈陣中，把P+2行和P+2列，全部乘以-1。然后去掉全1列，再把矩陣中-1改成2，便構(gòu)造出了N≠2s且N=2(p+1)(p為素?cái)?shù)且 p≡1(mod4))類型的2水平正交表〔7〕。

2．Lt2(tm)型正交表

Lt2(tm)型正交表一般可采用正交拉丁方法進(jìn)行構(gòu)造。先建立正交拉丁方完全組，然后帶入基本列，就可以得到相應(yīng)的正交表了。建立正交拉丁方完全組的方法有多種，常用的有如下兩種。

(1)素?cái)?shù)或素?cái)?shù)冪法

設(shè)n為任意一個(gè)素?cái)?shù)或素?cái)?shù)冪，GF(n)為n階有限域，而 a1，a2，…，an與 b1，b2，…，bn為域中全體元素任意兩列排列;c1，c2，…，ci為域中的全體非零元素，則可建立n階正交拉丁完全組為

為構(gòu)成的n階正交拉丁方完全組〔8〕。

對建立的正交拉丁方完全組帶入Li2(tm)型前兩列的基礎(chǔ)列便構(gòu)造出了Li2(tm)型的正交表，即

3．水平數(shù)為素?cái)?shù)冪的混合正交表

水平數(shù)為素?cái)?shù)冪的混合正交表的構(gòu)造通常采用并列法，是一種由標(biāo)準(zhǔn)表構(gòu)造水平不同正交表，可以安排水平數(shù)不等的正交試驗(yàn)的常用方法。數(shù)學(xué)原理是并列性定理。

設(shè)Ln(S1×S2×…×Sn)型正交表A=(λij)的第j1和j2列(j1≠j2)的交互組{第 L1列，第 L2列，……，第Lι列}是完備的，將第j1和j2列按下列新規(guī)則φ合并成“新列”:

這里，φ 是集合 Γ ={(u，v):u=1，2，…，Sj1，v=1，2，…，Sj2}到結(jié)合 P={1，2，…，Sj1× Sj2}上的一一映射。則新列~λ和A中去除第j1，j2，L1，…，Lι后一起組成的矩陣就是利用并列法所構(gòu)造的混合型正交表。

SAS實(shí)現(xiàn)

1．N=2s型的2水平正交表的SAS實(shí)現(xiàn)

為了更加通俗介紹SAS軟件實(shí)現(xiàn)哈達(dá)瑪矩陣直積法構(gòu)造2水平正交表過程，本文以L16(215)為例進(jìn)行詳細(xì)闡述，運(yùn)用SAS中IML模塊創(chuàng)建標(biāo)準(zhǔn)哈陣，然后將標(biāo)準(zhǔn)哈陣去除第一列即全1列，且將-1寫成2，便得到了L16(215)的正交表，結(jié)果見表1。SAS程序與解釋如下，數(shù)據(jù)集zj2即為所求的正交表L16(215):

%let N=16;/*宏變量N表示試驗(yàn)次數(shù)*/

%let J=a@a@a@a;/*符號(hào)@表示矩陣直積*/

proc iml;/*運(yùn)用IML模塊創(chuàng)建標(biāo)準(zhǔn)哈達(dá)瑪矩陣*/

a={1 1，1－1};/*a為最簡單的二階哈達(dá)瑪矩陣*/

d=＆J;

create zj from d;/*用create語句的from選項(xiàng)創(chuàng)建SAS數(shù)據(jù)集zj，即所求的標(biāo)準(zhǔn)哈陣*/

append from d;

close zj;

quit;

data zj2(drop=i col1);/*將標(biāo)準(zhǔn)哈陣去除第一列，且將-1寫成2，便得到了L16(215)的正交表數(shù)據(jù)集zj2即為所求*/

set zj;

array col(＆N);

do i=2 to＆N;

if col{i}=－1 then col{i}=2;

end;run;

proc sort;

by col2 col3;

run;

在SAS程序中只需修改兩個(gè)宏變量N(試驗(yàn)次數(shù))和J(哈達(dá)瑪矩陣直積)，就可以構(gòu)造 L4(23)、L8(27)、L32(231)等等N=2s型的2水平正交表。例如正交表L8(27)，此時(shí)N=8時(shí)，J=a@a@a，即 a的個(gè)數(shù)等于s，a表示最簡單的哈陣，a與a之間用@連接。

表1 L16(215)正交設(shè)計(jì)表

2．N≠2s且 N=2(p+1)(p為素?cái)?shù)且 p≡1(mod4))型的2水平正交表的SAS實(shí)現(xiàn)

以L12(211)為例，根據(jù) N=2(p+1)得 p=5，2為有限域GF(5)的一個(gè)元根，SAS程序構(gòu)造思路為:構(gòu)造R陣→K陣→H12陣→正交表L12(211)，生成數(shù)據(jù)集zj即為所求的正交表L12(211)，結(jié)果見表2，程序如下。

%let n=12;/*宏變量n為試驗(yàn)次數(shù)*/

%let yg=2/*宏變量yg代表有限域的最小元根*

%let k=%eval(＆n-1);/*宏變量 k為試驗(yàn)因素，由于構(gòu)造的是飽和正交表，其值也等于n-1*/

%let p=%eval(＆n/2-1);/* 宏變量 p，等同于GF(p)中的p，其值等于n/2-1*/

%let p1=%eval(＆p+1);

proc iml;/*根據(jù)R陣的數(shù)學(xué)排列規(guī)律運(yùn)用proc iml模塊構(gòu)造數(shù)據(jù)集a*/

a=j(＆p，＆p，0);do i=1 to ＆p;do j=1 to ＆p;

if i＞ j then a［i，j］=mod((i-j)，＆p);

else;if j＞ =i then a［i，j］=mod((＆p+i-j)，＆p);

end;end;create a from a;append from a;close a;quit;

data R(keep=b1-b＆p);/*利用數(shù)據(jù)集a生成數(shù)據(jù)集R，即為R陣*/

set a;array a{＆p};array col{＆p};array b{＆p};do i=1 to＆p;if col{i}^=0 then do;n=1;do until(y=col{i});y=mod(＆yg．＊＊n，＆p．);n+l;end;a{i}=mod(n －1，2);end;else if col{i}=0 then a{i}=．;if a{i}=0 then b{i}=1;else if a{i}=1 then b{i}=-1;else if a{i}=．then b{i}=0;end;run;

data c(drop=i);/*創(chuàng)建K陣中除R以外的部分*/

b0=0;array b{＆p};do i=1 to ＆p;b{i}=1;end;

run;

data K;/*合并數(shù)據(jù)集c、R生成K陣*/

set c R;if b0=．then b0=1;

run;

proc iml;/*生成單位矩陣E陣*/

I=I(＆p1);create E from I;append from I;quit;

data KE1(keep=k11-k1＆p1 k21-k2＆p1);/* 數(shù)據(jù)集KE1為H12的上部分，即K+E和K－E*/

merge K E;array b{＆p1}b0-b＆p;array col{＆p1};array k1{＆p1};array k2{＆p1};do i=1 to＆p1;k1{i}=b{i}+col{i};k2{i}=b{i}-col{i};end;k21=k21*-1;

run;

data KE2(drop=b0 － b＆p col1 － col＆pli);/* 數(shù)據(jù)集KE1為H12的下部分，即K-E和-K-E*/

merge K E;array b{＆p1}b0-b＆p;array col{＆p1};array k1{＆p1};array k2{＆p1};do i=1 to ＆p1;k1{i}=b{i}-col{i};k2{i}=-b{i}-col{i};if_n_=1 then k1{i}=k1{i}*-1;if_n_=1 then k2{i}=k2{i}*-1;end;k21=k21*-1;run;

data KE;/*合并KE1、KE2得KE，即為方陣H12*/

set KE1 KE2;

run;

proc transpose data=ke out=ke;/*對方陣H12進(jìn)行適當(dāng)?shù)霓D(zhuǎn)換，去掉全1列*/

run;

data ke(drop=c_NAME_);set ke;c=sum(of col1-col＆n);if c= ＆n then delete;

run;

proc transpose data=ke out=ke(drop=_NAME_);

run;

data zj(drop=i);/*數(shù)據(jù)集zj即為所求*/

set ke;array col{＆k};do i=1 to ＆k;if col{i}=-1 then col{i}=2;end;run;

proc sort;by col1 col2 col3;

run;

在SAS程序中只需修改宏變量n(試驗(yàn)次數(shù))，就可以構(gòu)造其他試驗(yàn)次數(shù)但滿足GF(p)的一個(gè)元根為2的N≠2s且N=2(p+1)(p為素?cái)?shù)且p≡1(mod4))型2水平正交表，如 L26(225)、L36(235)、L60(259)等。當(dāng)有限域的元根不包含2時(shí)，則需在修改宏變量n的基礎(chǔ)上，再修改宏變量yg(元根)即可。

3．Lt2(tm)型正交表的SAS實(shí)現(xiàn)

為了SAS操作的簡便，采用有限域的一個(gè)元根法建立通用的n階正交拉丁方完全組，并以L25(56)為例介紹SAS過程。其中，SAS生成的數(shù)據(jù)集LT即為所求，結(jié)果見表3，程序與注解如下。

表2 L12(211)正交設(shè)計(jì)表

%let t=5;/*宏變量t即L_(t2)(tm)中的t值*/

%let t1=%eval(＆t+1);

%let fz=%eval(＆t－1);/*宏變量fz表示正交拉丁方完全組方陣的個(gè)數(shù)，其值等于t-1*/

%macro iml;/*結(jié)合原理，運(yùn)用宏程序和proc iml過程構(gòu)造四個(gè)方陣*/

%do g=1%to＆fz;

proc iml;a=j(＆t，＆t);do i=1 to ＆t;

a［1，i］=i;/* 創(chuàng)建每個(gè)方陣的第一行*/

end;x=＆t*2-1;do i=2 to＆t;/*依據(jù)原理創(chuàng)建方陣的其他行*/

do j=1 to＆t;do k=3 to x;

m=i-1;l=j+1;

if i+j=k ＆ j^= ＆t then a［i，j］=a［m，l］+ ＆g-1;

if a［i，j］=0 then a［i，j］= ＆t;

if a［i，j］＞ ＆t then a［i，j］=mod(a［i，j］，＆t);

a［i，＆t］=A［1，+］-A［i，+］+a［i，＆t］;end;end;end;create c＆g from a;append from a;close c＆g;quit;

%end;%mend iml;%iml;

%macro cl;/*調(diào)用宏程序?qū)roc iml創(chuàng)建的每個(gè)方陣均轉(zhuǎn)換成列*/

%do k=1%to ＆fz;proc transpose data=C＆k．out=CC＆k．(drop=_NAME_);var_all_;run;%do i=1%to ＆t;data CC＆k．＆i．(keep=col＆i rename=(col＆i=col1));set CC＆k．;run;%end;data CC＆k．(rename=(col1=C＆k．));set%do i=1%to ＆t;CC＆k．＆i．

%end;;run;%end;%mend cl;%cl;

data CC＆t．;/*創(chuàng)建圖1中的基本列——第一列和第二列*/

do BASE1=1 to＆t;do BASE2=1 to＆t;output;end;end;

run;

%macro hb;/*將基本列的數(shù)據(jù)集和方陣轉(zhuǎn)換成列的各個(gè)數(shù)據(jù)集合并成一個(gè)數(shù)據(jù)集LT，即為所求的正交表*/

data LT;length base1 base2 8．;merge%do i=1%to ＆t．;CC＆i．%end;;

run;

%mend hb;

%hb;

表3 L25(56)正交設(shè)計(jì)表

在程序中只需修改宏變量t的值即可得到t為素?cái)?shù)的Lt2(tm)類型正交表。但當(dāng)t不為素?cái)?shù)時(shí)，基于生成Lt2(tm)類型正交表的正交拉丁方結(jié)構(gòu)不同，其SAS程序更為復(fù)雜，這里暫不介紹。

4.混合正交表的SAS實(shí)現(xiàn)

本文以混合正交表L16(43×26)(即L16(4k×2m)型混合正交表)為例，介紹如何通過水平數(shù)相等的正交表L16(215)通過“并列法”進(jìn)行構(gòu)造，結(jié)果見表4，SAS程序與注解如下。%let K4=3;/*宏變量K4表示水平數(shù)為4的因素個(gè)數(shù)*/

data hh1(drop=i col1 col2 col3 col4);/*運(yùn)用“并列法”將L16(215)進(jìn)行改造*/

length col23 col59 col611 col810 col712 6．;set zj2;/*調(diào)用哈達(dá)瑪矩陣法中生成的數(shù)據(jù)集zj2，即L16(215)*/

col23=col3*(col2=1)+(col2+col3)*(col2=2);

col59=col9*(＆K4＞ =2 and col5=1)+(col5+col9)*(＆K4＞ =2 and col5=2);

if col59^=0 then do;col5=0;col9=0;col13=0;end;col611=col11*(＆K4＞ =3 and col6=1)+(col6+col11)*(＆K4＞ =3 and col6=2);

if col611^=0 then do;col6=0;col11=0;col16=0;end;col810=col10*(＆K4＞ =4 and col8=1)+(col8+col10)*(＆K4＞ =4 and col8=2);

if col810^=0 then do;col8=0;col10=0;col15=0;end;col712=col12*(＆K4=5 and col7=1)+(col7+col12)*(＆K4=5 and col7=2);

if col712^=0 then do;col7=0;col12=0;col14=0;end;

run;

proc transpose data=hh1 out=hh2;var_all_;

run;

data hh1;set hh2;if col1=0 then delete;

run;

proc transpose data=hh1 out=hh2;

var_all_;

run;

data hhzj(drop=_NAME__LABEL_);

set hh2;

if_NAME_=“_NAME_”then delete;

run;

表4 L 16(43×26)正交設(shè)計(jì)表

上述程序中混合正交表L16(4k×2m)，當(dāng)K不同時(shí)，程序中只需修改宏變量K4即可，K4的取值為整數(shù)，范圍是〔1，5〕。不同的試驗(yàn)次數(shù)，水平數(shù)為素?cái)?shù)冪的混合正交表也可參照同樣的原理及方法生成。

結(jié) 論

構(gòu)造正交表的方法繁多，但不同類型的正交表的構(gòu)造原理與方法往往不同。本文根據(jù)不同類型正交表的構(gòu)造原理，完成了四種類型正交表SAS通用構(gòu)造程序。這樣，可以解決正交試驗(yàn)設(shè)計(jì)中獲得相應(yīng)正交表的困擾，方便實(shí)用，具有較強(qiáng)的應(yīng)用價(jià)值。

1．Zhang YS，Li WG，Mao SS，et al．Orthogonal arrays obtained by generalized difference matrices with g levels．SCIENCE CHINA，2011，54:133-143．

2．Ma CX，F(xiàn)ang KT，Erkki Liski．A new approach in constructing orthogonal and nearly orthogonal arrays．Metrika，2000，50:255-268．

3．Aloke Dey，Midha CK．Construction of some asymmetrical orthogonal arrays．Statistics ＆ Probability Letters，1996，28:211-217．

4．Liu ZW．A survey of orthogonal arrays of strength two．Acta Mathematicae Applicatae Sinica，1995:308-317．

5．Man VM．Nguyen．Some new constructions of strength 3 mixed orthogonal arrays．Journal of Statistical Planning and Inference，2008，138:220-233．

6．Zhang YS，Pang SQ，Wang YP．Orthogonal arrays obtained by generalized Hadamard product．Discrete Mathematics，2001，238:151-170．

7．Sloane NJA，Hedayat AS，John Stufken．Orthogonal arrays:theory and applications．Springer-verlag New York Inc，1999．

8．楊子胥．正交表的構(gòu)造．山東:山東人民出版社，1978．