河南師范大學(xué)數(shù)學(xué)與信息科學(xué)學(xué)院(453007)
龐善起 鹿姍姍
正交表的構(gòu)造方法及Matlab實(shí)現(xiàn)*
河南師范大學(xué)數(shù)學(xué)與信息科學(xué)學(xué)院(453007)
龐善起 鹿姍姍
正交試驗(yàn)設(shè)計(jì)是研究多因素水平的試驗(yàn)設(shè)計(jì)方法,它是一種高效率、快速、經(jīng)濟(jì)的試驗(yàn)設(shè)計(jì)方法。正交表在衛(wèi)生、統(tǒng)計(jì)、醫(yī)藥、農(nóng)業(yè)、制造業(yè)等方面都有非常重要的應(yīng)用,尤其是高強(qiáng)度正交表可以研究較多因素之間的相互作用,但是正交表的強(qiáng)度越高越難構(gòu)造;另外,混合正交表的魅力在于它有較大的靈活性,允許試驗(yàn)因素具有不同水平數(shù),得到了越來越多的關(guān)注。因此如何快速、高效地得到需要的正交表是有必要和有意義的。
至今為止,已經(jīng)有許多組合數(shù)學(xué)家和統(tǒng)計(jì)學(xué)家致力于正交表的研究,例如楊子胥、龐善起、張應(yīng)山、杜蛟等分別在文獻(xiàn)[1-8]中給出了不同類型的正交表的構(gòu)造方法,這些正交表不僅在試驗(yàn)設(shè)計(jì)中有廣泛的應(yīng)用,隨著研究的不斷深入,它們還可以被應(yīng)用于編碼學(xué)、密碼學(xué)、計(jì)算機(jī)科學(xué)中。
1.Ltu(tm)型正交表[1]
(1,0,0,…,0,0,0)
1個(gè)
(bu-1,1,1,0,…,0,0,0)
t個(gè)
(bu-2,1,bu-2,2,1,…,0,0,0)
t2個(gè)
?
?
(b2,1,b2,2,b2,3,…,b2,u-2,1,0)
tu-2個(gè)
(b1,1,b1,2,b1,3,…,b1,u-2,b1,u-1,1)
tu-1個(gè)
因?yàn)槊總€(gè)bij為域GF(t)的所有元素,即有t種取法,所以以上的向量,就是Vu的全部標(biāo)準(zhǔn)向量,共有
構(gòu)造Ltu(tm)型正交表的方法和步驟如下:
第1步,先給出t階有限域GF(t);
第3步,以全體向量為行號(hào),以全部的標(biāo)準(zhǔn)向量為列號(hào),做出所有內(nèi)積即可得Ltu(tm)型正交表。
2.混合水平正交表
混合水平正交表[1]可由Ltu(tm)型正交表通過并列而產(chǎn)生。構(gòu)造混合水平正交表的方法和步驟如下:
第1步,可借助上文Ltu(tm)型正交表的方法得到Ltu(tm)型正交表;
第2步,從Ltu(tm)中任意取定k個(gè)獨(dú)立的列(k
第4步,再將這k列同行元素所構(gòu)成的tk個(gè)互異的有序組,按字典排列法與tk個(gè)自然數(shù)1,2,…,tk建立一一對應(yīng),而把這k列所構(gòu)成的每個(gè)有序組均換成它對應(yīng)的自然數(shù),便得到一個(gè)新的tk水平列;
3.高強(qiáng)度正交表
一個(gè)第j 列的元素是0,1,…,qj-1的ω×n矩陣A,稱為一個(gè)強(qiáng)度為2的正交表,如果滿足以下條件[4]:
(1)每一列中每個(gè)元素出現(xiàn)的次數(shù)相同;
(2)在任兩列ai,aj(1≤i,j≤n)中的每一個(gè)數(shù)對(0,0),…,(0,qj-1),(1,0),…,(1,qj-1),…,(qj-1,0),…,(qj-1,qj-1)出現(xiàn)的次數(shù)相同。
只滿足上述條件(1)的矩陣叫做強(qiáng)度1的正交表。如果A的任何ω×d的子矩陣包含所有可能的1×d行向量,且這些行向量出現(xiàn)的次數(shù)相同,稱為強(qiáng)度為d的正交表。強(qiáng)度d≥3的正交表稱為高強(qiáng)度正交表[4-6]。
引理1 設(shè)M是一個(gè)選自Qn的任意t(1≤t≤qn-1)個(gè)行向量構(gòu)成的矩陣,則
是一個(gè)強(qiáng)度為1的正交表。
引理2 對于任意的整數(shù)1≤k≤q總存在Qn上的1-正交分劃A1,A2,…,Ak。
證明:由于Qn本身就是一個(gè)正交表,所以當(dāng)k=1時(shí)定理顯然成立; 當(dāng)k 由于Qn中有qn個(gè)q元n維向量,考慮將Qn分為行數(shù)相等的分劃容量為q的A1,A2,…,Ak,其中的每個(gè)Ai(1≤i≤k)中都有qn-1個(gè)向量,為了得到Ai: ? 第2q-1步,將剩下的qn-1個(gè)向量,組成矩陣Aq,它當(dāng)然是一個(gè)強(qiáng)度1的正交表。 這樣就構(gòu)造出了Qn的一個(gè)1-正交分劃A1,A2,…,Aq。 引理3 若有Qn上的正交分劃A1,A2,…,Aq,如果Ai(i=1,2,…,q)都是qn-1行n列的強(qiáng)度為t的正交表,那么下列矩陣都是強(qiáng)度為t+1的正交表,并且是Qn+1的(t+1)-正交分劃。 其中的0,1,2,…,q-1均指與其對應(yīng)的Ai有相同的行的元素全為0,1,2,…,q-1的向量,下文同。 高強(qiáng)度正交表的構(gòu)造可由引理1,2得到強(qiáng)度為1的正交表,和Qn的1-正交分劃,在這個(gè)前提下,再由引理3得到強(qiáng)度為2的正交表,和Qn+1的2-正交分劃,接著由引理3得到強(qiáng)度為3的正交表,和Qn+2的3-正交分劃,依次下去,逐步構(gòu)造出高強(qiáng)度的正交表。 為了更加直觀詳細(xì)地說明Matlab程序的可行性,下面我們將選取幾個(gè)較為簡單的正交表作為例子說明利用Matlab程序?qū)崿F(xiàn)正交表的構(gòu)造方法。 1.Ltu(tm)型正交表 clc;clear t=2;u=3;%初始條件 p=t;V=[ ];K=[ ];L=[ ];D=[ ];L1=[ ];L2=[ ]; G=transpose(0:(t-1));%t階有限域 i=1; while(i<=u) R=kron(kron(ones(t^(i-1),1),G),ones(t^(u-i),1)); V=[V,R]; i=i+1; end V%所有的u維向量 U=rref(V′);%標(biāo)準(zhǔn)化 [m,n]=find(U==1); [s,t]=find(U((m+1):u,n)<=1); B=U(:,unique(n(t)))′;%標(biāo)準(zhǔn)向量 L=V*B′; L=mod(L,p)%得到正交表L8(27) 在Matlab程序中只需修改初始條件t,u的值,我們就可以生成形如Ltu(tm)的正交表。 2.混合水平正交表 以L8(4×23)為例,則23=8,此時(shí)t=2,u=3,k=2,由上文的程序得到L8(27),取定2個(gè)獨(dú)立的列(2<3),去掉能用這2列線性表示的所有列,再將這2列所構(gòu)成的22=4個(gè)互異的有序組,把這2 列所構(gòu)成的每個(gè)有序組均換成它對應(yīng)的自然數(shù),從而生成正交表L8(4×23),結(jié)果如表2,Matlab程序與解釋如下: 表1 正交設(shè)計(jì)表L8(27) clc;clear t=2;u=3;k=3;%初始條件 p=t;[z,z1]=size(L);%接上節(jié)程序 K=L(:,1:k);%默認(rèn)1:k列,也可以換成其他列 i=1; while(i<=z1) d=rank([K,L(:,i)]′)%m個(gè)n維向量線性相關(guān)?秩 D=[D,d]; i=i+1; end D a=find(D(1,:)>=(k+1)) L1=L(:,a)%得到正交表Ltu(tr) %以下生成混合水平列 Y=unique(K,′rows′)%刪除K中重復(fù)的行,得到有序組 [y,y1]=size(Y); j=1; while(j<=y) b=find(K(:,1)==Y(j,1)&K(:,2)==Y(j,2)); [z,z1]=size(b); K(b,:)=(y-(j-1))*ones(z,k);%有序組與自然數(shù)建立一一對應(yīng) j=j+1; end L2=[flipud(K(:,1)),L1]%得到混合水平正交設(shè)計(jì)表 3.高強(qiáng)度正交表 以L8(24)為例,則Q={0,1},q=2,N=4,t=3,首先生成強(qiáng)度為1的正交表,和Q2的1-正交分劃,再生成強(qiáng)度為2的正交表,和Q3的2-正交分劃,接著生成強(qiáng)度為3的正交表,和Q4的3-正交分劃,這樣就可以生成2水平強(qiáng)度為3的正交表。結(jié)果如表3,和表4,Matlab程序與解釋如下: 表2 正交設(shè)計(jì)表L8(4×23) clc;clear; q=2;%水平(初始條件) t=3;%強(qiáng)度(初始條件) N=4;%列數(shù)(初始條件) n=N-(t-1);%q水平強(qiáng)度為1的正交表的列數(shù) Q=[];a0=[];a1=[];a=[];A1=[];A=[];L=[];L0=[];L1=[]; %第一部分:生成Q^n即所有的q元n維向量 i=1; while(i<=n) R=kron(kron(ones(q^(i-1),1),[0:(q-1)]′),ones(q^(n-i),1)); Q=[Q,R];%Q^n即所有的q元n維向量 i=i+1; end %第二部分:生成1-正交分劃 j=1; while(j<=(q-1)) for i=0:(q-1); a0=Q((1+q^(n-2)*(j-1)):q^(n-2)*j,:); a1=mod(i+a0,q);%正交分劃 a=[a;a1]; end Q=setdiff(Q,repmat(a,2,1),′rows′);%Q^n-a j=j+1; end A=[a;Q];%生成強(qiáng)度為1的正交表 %第三部分:正交分劃的遞歸構(gòu)造 while(n<=N-1)%Q正交分劃q^(n-2)行向量的個(gè)數(shù) L0=kron(ones(q,1),kron([0:q-1]′,ones(q^(n-1),1))); j=0; while(j<=q-1) A1=circshift(A,[q^(n-1)*j,0]);%每次向下移q^(n-1) L1=[L1;A1]; j=j+1; end A=[L0,L1];%所有正交分劃的并列 L1=[]; n=n+1; end %第四部分:生成強(qiáng)度為t的正交表 for i=1:q; L=A((1+q^(n-1)*(i-1)):q^(n-1)*i,:)%強(qiáng)度為t的正交表 end 在Matlab程序中,保證n=N-(t-1)≥2的前提下,只需修改初始條件q,N,t的值,我們就生成q水平,強(qiáng)度為t,列數(shù)為N,行數(shù)為q^(N-1)的正交表。比如修改初始條件,使得q=4,t=3,N=4便可得到水平為4,強(qiáng)度為3,列數(shù)為4,行數(shù)為64的正交表;修改初始條件,使得q=5,t=3,N=4,便可得到水平為5,強(qiáng)度為3,列數(shù)為6,行數(shù)為125的正交表。 表3 正交設(shè)計(jì)表L8(24) 表4 正交設(shè)計(jì)表L8(24) 正交表研究正在引起眾多領(lǐng)域研究者的廣泛關(guān)注。本文在原有的理論基礎(chǔ)上使用Matlab 生成三種類型的正交表,但是仍有一些問題值得改進(jìn)和研究,本文只給出了生成單一水平的高水平高強(qiáng)度正交表構(gòu)造方法,所以接下來還需要繼續(xù)尋找混合水平的高水平高強(qiáng)度正交表的構(gòu)造方法。 [1]楊子胥.正交表的構(gòu)造.濟(jì)南:山東人民出版社,1978. [2]龐善起.正交表的構(gòu)造及其應(yīng)用.成都:電子科技大學(xué),2004. [3]張應(yīng)山.正交表的數(shù)據(jù)分析及其構(gòu)造.上海:華東師范大學(xué),2006. [4]杜蛟,王守印,溫巧燕,等.強(qiáng)度m的對稱正交表的遞歸構(gòu)造.應(yīng)用數(shù)學(xué)學(xué)報(bào),2012,35(2):232-244. [5]Jiang L,Yin JX.An approach of constructing mixed-level orthogonal arrays of strength≧3.Science China(Mathematics),2013,6(56):1109-1115. [6]Du J,Wen QY,Zhang J.New Construction of Symmetric Orthogonal Arrays of Strength t.IEICE Transactions on Fundamentals of Electronics Communications and Computer Sciences,2013,9(96):1901-1904. [7]李學(xué)斌,余小領(lǐng),李斌.運(yùn)用Excel進(jìn)行正交表的構(gòu)建.河南科技學(xué)院學(xué)報(bào)(自然科學(xué)版),2008,36(4):115-118. [8]陳遠(yuǎn)方,林曦晨,徐利華,等.常用正交表的構(gòu)造原理及SAS實(shí)現(xiàn).中國衛(wèi)生統(tǒng)計(jì),2012,29(4):470-474. (責(zé)任編輯:劉 壯) 國家自然科學(xué)基金(11571094)正交表構(gòu)造的Matlab實(shí)現(xiàn)
結(jié) 論