常用正交表的構造原理及SAS實現

華中科技大學同濟醫學院公共衛生學院流行病與衛生統計學系(430030)

陳遠方 林曦晨 徐利華 湯洪秀 汪宏晶 尹 平△

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

方法和原理

1．2水平正交表

(1)哈達瑪矩陣直積法構造N=2s型

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

則A與B的直積也是哈達瑪矩陣。應用上述構造標準哈達瑪矩陣中，除全1列外，其任意兩列與全1列都正交。因此，去掉全1列，就得到N=2s型的2水平正交表〔6〕。

(2)哈達瑪矩陣特征函數法構造N≠2s且N=2(p+1)(p為素數且p≡3(mod4))型

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

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

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

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

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

2．Lt2(tm)型正交表

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

(1)素數或素數冪法

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

為構成的n階正交拉丁方完全組〔8〕。

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

3．水平數為素數冪的混合正交表

水平數為素數冪的混合正交表的構造通常采用并列法，是一種由標準表構造水平不同正交表，可以安排水平數不等的正交試驗的常用方法。數學原理是并列性定理。

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

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

SAS實現

1．N=2s型的2水平正交表的SAS實現

為了更加通俗介紹SAS軟件實現哈達瑪矩陣直積法構造2水平正交表過程，本文以L16(215)為例進行詳細闡述，運用SAS中IML模塊創建標準哈陣，然后將標準哈陣去除第一列即全1列，且將-1寫成2，便得到了L16(215)的正交表，結果見表1。SAS程序與解釋如下，數據集zj2即為所求的正交表L16(215):

%let N=16;/*宏變量N表示試驗次數*/

%let J=a@a@a@a;/*符號@表示矩陣直積*/

proc iml;/*運用IML模塊創建標準哈達瑪矩陣*/

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

d=＆J;

create zj from d;/*用create語句的from選項創建SAS數據集zj，即所求的標準哈陣*/

append from d;

close zj;

quit;

data zj2(drop=i col1);/*將標準哈陣去除第一列，且將-1寫成2，便得到了L16(215)的正交表數據集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程序中只需修改兩個宏變量N(試驗次數)和J(哈達瑪矩陣直積)，就可以構造 L4(23)、L8(27)、L32(231)等等N=2s型的2水平正交表。例如正交表L8(27)，此時N=8時，J=a@a@a，即 a的個數等于s，a表示最簡單的哈陣，a與a之間用@連接。

表1 L16(215)正交設計表

2．N≠2s且 N=2(p+1)(p為素數且 p≡1(mod4))型的2水平正交表的SAS實現

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

%let n=12;/*宏變量n為試驗次數*/

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

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

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

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

proc iml;/*根據R陣的數學排列規律運用proc iml模塊構造數據集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);/*利用數據集a生成數據集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);/*創建K陣中除R以外的部分*/

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

run;

data K;/*合并數據集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);/* 數據集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);/* 數據集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進行適當的轉換，去掉全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);/*數據集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(試驗次數)，就可以構造其他試驗次數但滿足GF(p)的一個元根為2的N≠2s且N=2(p+1)(p為素數且p≡1(mod4))型2水平正交表，如 L26(225)、L36(235)、L60(259)等。當有限域的元根不包含2時，則需在修改宏變量n的基礎上，再修改宏變量yg(元根)即可。

3．Lt2(tm)型正交表的SAS實現

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

表2 L12(211)正交設計表

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

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

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

%macro iml;/*結合原理，運用宏程序和proc iml過程構造四個方陣*/

%do g=1%to＆fz;

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

a［1，i］=i;/* 創建每個方陣的第一行*/

end;x=＆t*2-1;do i=2 to＆t;/*依據原理創建方陣的其他行*/

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;/*調用宏程序將proc iml創建的每個方陣均轉換成列*/

%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．;/*創建圖1中的基本列——第一列和第二列*/

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

run;

%macro hb;/*將基本列的數據集和方陣轉換成列的各個數據集合并成一個數據集LT，即為所求的正交表*/

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

run;

%mend hb;

%hb;

表3 L25(56)正交設計表

在程序中只需修改宏變量t的值即可得到t為素數的Lt2(tm)類型正交表。但當t不為素數時，基于生成Lt2(tm)類型正交表的正交拉丁方結構不同，其SAS程序更為復雜，這里暫不介紹。

4.混合正交表的SAS實現

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

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

length col23 col59 col611 col810 col712 6．;set zj2;/*調用哈達瑪矩陣法中生成的數據集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)正交設計表

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

結 論

構造正交表的方法繁多，但不同類型的正交表的構造原理與方法往往不同。本文根據不同類型正交表的構造原理，完成了四種類型正交表SAS通用構造程序。這樣，可以解決正交試驗設計中獲得相應正交表的困擾，方便實用，具有較強的應用價值。

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，Fang 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．楊子胥．正交表的構造．山東:山東人民出版社，1978．