具有下界約束混料試驗的Monte-Carlo漸近最優(yōu)設(shè)計

楊曉珍

（凱里學院理學院 貴州凱里 556011）

自20世紀二十年代初，由Fisher.R.A等人提出試驗設(shè)計的統(tǒng)計方法，許多統(tǒng)計學家根據(jù)不同的情況提出了各種試驗設(shè)計方案，例如，正交試驗設(shè)計，均勻試驗設(shè)計，回歸混料試驗設(shè)計等[1-2]。這些重要的設(shè)計方法已在社會生產(chǎn)和實踐中得到了廣泛的應用和推廣。

近年來，最優(yōu)試驗設(shè)計逐漸發(fā)展，稱為一門新學科，基于最優(yōu)試驗設(shè)計所產(chǎn)生的數(shù)據(jù)不僅能夠便于統(tǒng)計推斷，而且還能簡化計算，提高效率。最優(yōu)試驗設(shè)計的應用非常廣泛，遍及工農(nóng)業(yè)的生產(chǎn)和科學實驗的各項領(lǐng)域[3-4]。其中D-準則和MV-準則是評價最優(yōu)設(shè)計的常用準則，關(guān)于D-準則下的最優(yōu)設(shè)計也有很多的推廣。然而，對于含有較多未知變量的試驗設(shè)計，由于難以求得其信息矩陣的行列式或逆矩陣的顯式表達，所以難以得到各類最優(yōu)準則下的精確最優(yōu)設(shè)計解。為此，許多學者提出了不同的算法，F(xiàn)edorov最先提出了D-最優(yōu)設(shè)計的算法[5]，Evans提出了Dn-最優(yōu)設(shè)計的直接擴張法[6]，1974年，Mitchell給出了D-最優(yōu)設(shè)計的DETMAX算法[7]，通過增加或刪除設(shè)計點的方式使得信息矩陣行列式達到最大，直到它滿足試驗的精度要求.后來，Galil和Kiefer對該算法進行了優(yōu)化[8]，進一步節(jié)省了算法空間和時間。這一算法能有效獲得設(shè)計域內(nèi)的D-最優(yōu)準則下的設(shè)計點，并且誤差可控。對于其他一些準則，例如MV-準則，G-準則，Ds-準則，E-準則[9]等等，關(guān)于MV-準則下的有效算法的文獻不多，而不同的準則下所涉及到的算法也不同，在使用的過程中需要區(qū)別使用。

許多文獻都嘗試構(gòu)造隨機搜索算法來得到最優(yōu)設(shè)計[10-14]，本文在此基礎(chǔ)上，嘗試在D-準則和MV-準則下，使用相同的方法——Monte-Carlo方法，通過大量隨機數(shù)模擬的方式，求出不同準則下的近似最優(yōu)解，并且通過驗證，說明模擬效果良好。

一、具有下界約束的混料模型

混料是指若干種不同成分的物質(zhì)混合或合成一種穩(wěn)定的物質(zhì)或產(chǎn)品。這些產(chǎn)品的每種成分的多少是用相對量表示的，這種相對量就是所用成分在總量中所占比例。很多情形下，我們所需要考慮的是完全混料模型，即各成分所占比例都要大于零。有甚者，通過實踐經(jīng)驗可以首先規(guī)定各中成分所占比例的下界，這一類模型稱之為具有下界約束混料模型。

對于q分量有下界約束的混料系統(tǒng)，響應是各分量比例x1,x2,…,xq的函數(shù)，各個分量的約束下界分別為a1,a2,…,aq，由各分量比例所確定的q-1維單純形[15]可表示為:

η(xj,Θ) 是含有未知參數(shù)且形式已知的函數(shù)，Θ=(θ1,θ2,…,θm)T是未知待估的參數(shù)向量，ε=(ε1,ε2，…,εN)T為隨機誤差向量，并且假設(shè)E(εj)=0，Var(εj)=σ2，j=1,2,…,N，且Cov(εi,εj)=0(i≠j)。一般還假設(shè)隨機誤差向量服從正態(tài)分布，即有ε～N(0 ,σ2IN)，其中IN為N階單位陣。

記Y=(y1,y2,…,yN)T為試驗結(jié)果，(x1,x2,…,xN)∈X為給定的試驗條件，X為拓撲空間上的緊集，稱之為設(shè)計域。我們將N次試驗的試驗點記作，其中共有k個不同的設(shè)計點，記為x1,x2,…,xk.假設(shè)在試驗點xi處進行了ri次試驗，記wi=ri/N為試驗點xi的權(quán)系數(shù)。

稱ξ∈為一個具有測度的k點設(shè)計，稱為設(shè)計空間.如果函數(shù)η(x,Θ )是關(guān)于參數(shù)與各分量函數(shù)的一個線性組合，即可以表示為：

當det(M(ξ))≠0時，稱ξ為一個非奇異的設(shè)計。一般考慮的設(shè)計都假定為非奇異的。在試驗設(shè)計中，很多準則都是由關(guān)于信息矩陣的判定函數(shù)來衡量的。記在準則φ下的判定函數(shù)為Φφ(M(ξ))。

D-準則[1]的原理是將信息矩陣的行列式最大化，從而使得參數(shù)向量Θ的置信橢球體積最小。D-準則可表述為：

在D-準則下，關(guān)于信息矩陣的判定函數(shù)為:

若模型(2)可以表示為參數(shù)向量線性組合的形式，即:

以Mij表示M-1(ξ)中第i行j列的元素，則有

在MV-準則下，關(guān)于信息矩陣的判定函數(shù)為:

信息矩陣在最優(yōu)設(shè)計中扮演這重要的角色。當變量較多時，信息矩陣的行列式或是逆矩陣難以求得顯式結(jié)果，所以要解得精確的最優(yōu)設(shè)計是比較困難的。以下先使用Monte-Carlo方法來構(gòu)造隨機試驗點樣本，然后在隨機樣本中求取信息矩陣的行列式的近似極大值，并找到相對應的試驗點。

二、Monte-Carlo方法構(gòu)造試驗點樣本

對于(1)所示的q-1維單純形Sq-1，若Sq-1不是一點或是空集，則稱混料問題是非退化的有下界約束問題。構(gòu)成非退化的混料問題的充要條件是：

且各分量還應具有隱上界約束：

設(shè)A是非退化的有下界約束的混料單純形中的任意一點，我們可以通過線性變換的方式將Sq-1中的自然分量坐標轉(zhuǎn)換為無約束的混料單純形中的坐標。記點A的自然分量坐標及擬分量坐標分別為：

令Ιq為q階單位陣，a=(a1,a2,…,aq)T為約束下界值構(gòu)成的向量，再記1表示元素全為1的q維向量。

則兩種坐標間的變換關(guān)系為：

若D-1存在，(6)式的逆變換為：

這樣，經(jīng)過線性變換后擬分量的全集就構(gòu)成一個無約束的混料單純形，記作：

MC1步:生成n組無約束的隨機混料試驗點。

由U[0,1]生成n個隨機數(shù)，記作Z=(z,z,…,z)T。

111121n由于總體Z為連續(xù)型隨機變量，為了盡可能的使抽取的隨機數(shù)互不相同，我們可以將抽取隨機數(shù)的精度提高，例如要求抽取的隨機數(shù)保留小數(shù)點后6位，這樣以來能使得抽到相同的隨機數(shù)的概率變得很小。再令

Z1,Z2,…,Zq-1都是隨機向量，且滿足j=2,3,…,q-1；Zq是由前q-1個隨機向量所確定的常值向量。

矩陣Z=(Z1,Z2,…,Zq)是一個n×q階矩陣。容易驗證Z各行元素和為1，記矩陣Z的第k行為：

這里，τk中各個元素間并不獨立，對于任意的zjk(j=2,3,…,q-1)都是由前j-1個元素所共同確定的。τk中各元素的聯(lián)合密度函數(shù)可以表示為：

其中I{?}為示性函數(shù)。

又因為Z1中各元素相互獨立，由此可知Z2,Z3,…,Zq-1中各個元素也相互獨立。則各τk，k=1,2,…,n是相互獨立的，即有。

當n→∞時，τ1,τ2,…,τn能均勻地充斥在整個單純形內(nèi)部。記

是τk為中心δ為半徑的鄰域，對于任意的δ＞0，應有

MC2步：由(7)將無約束的隨機混料試驗點轉(zhuǎn)換為具有下界約束的擬分量坐標。

例如，三分量和四分量混料系統(tǒng)取下界約束值分別取a=(0.1,0.2,0.1)T，a=(0.3,0.2,0.1,0.1)T，由MC1步生成500個無約束隨機混料試驗點，再轉(zhuǎn)換為具有下界約束的擬分量坐標，使用文獻[16]中的映射方法，得到如下圖所示。

圖1 無約束的隨機混料試驗點(a)(c)；轉(zhuǎn)換后具有下界約束的擬分量坐標(b)(d)

MCS步:對于混料試驗k點設(shè)計，在生成的隨機點τ1,τ2,…,τn中隨機抽取k點組成一個樣本，重復抽取m次，記第j次抽取的樣本為在τ(j)下的設(shè)計為：

將生成的各組隨機數(shù)按組分別代入設(shè)計陣中，再計算信息陣為

表示第j個隨機樣本下的信息矩陣。

設(shè)若存在一個在測度(w1,w2,…,wk)最優(yōu)設(shè)計ξ*，對于關(guān)于信息矩陣的某種最優(yōu)準則Φ(M(ξ))，應有：

MCM步:在各M(ξj)(j=1,2,…,m)分別計算Φφ(M(ξj))值，按準則φ在ξ1,ξ2,…,ξm中選擇出最優(yōu)解，記作，將其作為全局最優(yōu)解ξ*的個種近似解。

這里需要說明的是，如果在可行域上存在一個精確的最優(yōu)解

由MCS步生成的隨機點集為Ω={τ1,τ2,…,τn}，n＞k.F(Ω)表示由Ω生成的離散Borel域。下面分兩種情況討論。將簡記為ξ*∈F(Ω)。

（1）如果ξ*∈F(Ω)，那么一次抽樣恰好取到ξ*的概率為，m次抽樣中能取到ξ*的概率為

理論上，可以由固定的α值確定需要抽樣的次數(shù)m(α)，例如當m＞m(0.95)時，有α＞0.95，也就是我們總能找到一個確定的值m(α)使得取到最優(yōu)設(shè)計點的概率充分的大。

（2）如果ξ*?F(Ω)，進一步假設(shè)在準則φ下的判定函數(shù)Φφ(M(ξ))是關(guān)于變量(τ1,τ2,…,τk)的 連續(xù)函數(shù)，令是關(guān)于τ*的球形鄰域，若對于任意的ε＞0，都存在一個常數(shù)δ(ε)＞0，使得當(τ1,τ2,…,τk)∈U(τ*,ε)時，有恒成立。則一次取樣下的試驗點組合τ=(τ1,τ2,…,τk)能落入U(τ*,ε)的概率為:

m次取樣有試驗點組合落入U(τ*,ε)的概率為α'=1-(1-p)m。所以當ε固定時，我們可以確定一個正整數(shù)m(ε)，使得當m＞m(ε)時，有成立。

雖然，不論是哪種情形，我們總可以得到這樣一個結(jié)論：就是當m充分大的時候，能保證取到最優(yōu)設(shè)計點或是隨機試驗點落入最優(yōu)試驗點的某個鄰域的概率趨近于1，但是存在一個嚴重的問題，當隨機點數(shù)n較大時，需要取的樣本數(shù)量會大得驚人，而n如果較小，就不能保證隨機試驗點充分接近最優(yōu)解。關(guān)于如何優(yōu)化隨機點數(shù)與抽樣次數(shù)的問題仍有待進一步研究。以下我們討論使用Monte-Carlo方法求出具有下界約束的三分量混料試驗的近似最優(yōu)解。

三、模擬近似最優(yōu)設(shè)計

對于形如(3)中所示的回歸模型的設(shè)計，按照Monte-Carlo方法構(gòu)造隨機試驗點樣本。經(jīng)過轉(zhuǎn)換后得到n個具有下界約束的隨機混料試驗點，令X為n×q階矩陣，其中每一行對應了一個試驗點的個分量坐標。然后將每一行元素代入設(shè)計陣，再由MCS步隨機取出m個樣本，計算出各組隨機樣本下的信息矩陣的行列式，并按升序排列為:

由于有隨機數(shù)精度的要求，由各detM(ξ(j))(j=1,2,…,m)組成的樣本中，一般不會出現(xiàn)打結(jié)的情形，如果該樣本中存在結(jié)點，只需滿足detM(ξ(m-1))＜detM(ξ(m))就可保證結(jié)點對結(jié)果不會造成影響，此時detM(ξ(n))的樣本秩就等于n。再由detM(ξ(n))在原樣本中的位置，例如detM(ξ(m))=detM(ξj)，那么就是對應的第j組隨機試驗點，使得樣本信息矩陣的行列式在n次模擬中達到最大，從而可以將τ(j)視作最優(yōu)設(shè)計點的一種近似.設(shè)若存在最優(yōu)精確設(shè)計為ξ*，應有detM(ξ(m))→ detM(ξ*)，τ(j)→τ*(n→ ∞)成立。

下面我們通過一個實例來說明Monte-Carlo方法在求解漸近最優(yōu)設(shè)計的有效性。

考慮三分量二階Scheffé多項式回歸模型

給定下界約束值向量為a=(0.1,0.2,0.1)T，在D-準則下，求解飽和設(shè)計的試驗點。

按照2中所給出的步驟，首先生成1000組無約束的混料隨機數(shù)，再將其轉(zhuǎn)換為具有下界約束的隨機數(shù)，由于設(shè)計的行列式值很小，為了便于比較，在計算行列式時將信息陣乘以100，即det(1 0 0M(ξj)).在MCS步隨機抽取10000個隨機試驗點的樣本，并重復MCM步若干次(由于隨機取樣不能完全取遍所有隨機數(shù)的組合)，以確定在各種組合下的最優(yōu)設(shè)計，經(jīng)過計算可得到：

假設(shè)模型(11)參數(shù)真值為:

在模型中加入標準差為0.1的高斯噪聲作為誤差項，所得y值作為試驗的結(jié)果，并假設(shè)在每個試驗點重復三次試驗，在所對應的試驗點下列出對應的響應值，結(jié)果如表1所示:

表1 飽和近似D最優(yōu)設(shè)計

此時響應曲面及等高圖如下圖所示。

圖2 飽和近似D最優(yōu)設(shè)計的響應曲面與等高圖

如果考慮帶有測度的7點設(shè)計，假設(shè)含測度的設(shè)計為:

通過計算可得到參數(shù)的最小二乘估計為:

記表2和表3的設(shè)計分別為ξ6和ξ7，以上結(jié)果中，兩種設(shè)計對一次項系數(shù)的估計值與真值偏離都不大，但交互項的系數(shù)估計值與真值之間都存在一定的偏離。導致這一結(jié)果的因素主要有兩個:一是源于隨機誤差項的干擾：二是由于試驗受約束條件限制，使得可行區(qū)域較小，所以得到的響應曲面的估計也是局部的。

表2 含測度的7點近似D最優(yōu)設(shè)計

表3 飽和近似MV最優(yōu)設(shè)計及試驗結(jié)果

下面考慮在MV-準則下的近似最優(yōu)解，此時關(guān)于信息矩陣的判定函數(shù)為:

與之前的討論類似，由MC1步和MC2步生成具有下界約束的隨機混料試驗點，再有MCS步隨機抽取其中的m個樣本，計算出各組隨機樣本的信息矩陣及對應的逆矩陣，再計算出各組隨機樣本下的信息矩陣的逆矩陣的跡，記為：

并按升序排列為:

與D-漸近最優(yōu)中的討論類似，由各trj(j=1,2,…,m)組成的樣本一般不會出現(xiàn)打結(jié)的情形，如果該樣本中存在結(jié)點，只需滿足tr(1)＜tr(2)，結(jié)點對結(jié)果就不會產(chǎn)生影響。此時tr(1)的樣本秩就等于1，再由tr(1)在原樣本中的位置，由tr(1)=trj=trM-1(ξj)確定對應的第j組隨機樣本，即使得。

此時，ξj下的設(shè)計點τ(j)可以作為MV-準則下最優(yōu)設(shè)計點的一種近似，設(shè)若存在最優(yōu)精確設(shè)計為ξ*，應有成立。我們?nèi)钥紤]模型(11)，假設(shè)回歸系數(shù)真值仍為(13)，生成1000個具有下界約束的混料隨機數(shù)，然后重復取樣10000次，重復MCM步若干次，計算出飽和設(shè)計下的的近似MV-最優(yōu)解，并在模型中加入標準差為0.1的高斯噪聲，由于信息矩陣中的元素較小，在循環(huán)過程中常會出現(xiàn)奇異設(shè)計，為避免這一情形，我們將設(shè)計矩陣乘以一個常數(shù)使之不易出現(xiàn)奇異情形，在此我們把信息矩陣乘以100再進行后續(xù)運算。

通過計算得到：

所得y值作為試驗的結(jié)果.若在每個試驗點重復三次試驗，在所對應的試驗點下列出對應的響應值，結(jié)果如表3所示:

圖3 飽和近似MV最優(yōu)設(shè)計的響應曲面與等高圖

四、討論

本文討論了Monte-Carlo方法來求解D-準則和MV-準則下的近似最優(yōu)解，并給出了實例模擬。雖然關(guān)于Scheffé多項式在各種最優(yōu)性下的設(shè)計的文獻已有很多，但對于含有約束條件的混料試驗設(shè)計，在不同的準則下求解最優(yōu)解的算法都是不同的，并且有的算法需要占用大量的內(nèi)存。而Monte-Carlo方法原理相對簡單，各種不同的最優(yōu)準則可視作“選取條件”，在大量隨機樣本中抽取出符合這一條件的最優(yōu)解，并且所得結(jié)果并不差，雖然這一結(jié)果并非精確解，但可以作為在該準則下最優(yōu)解的一種參考。

使用Monte-Carlo方法求解最優(yōu)準則下的精確解尚處于探索階段，有很多問題還不能有效的解決。比如，隨著分量的增加，需要生成的隨機樣本數(shù)就會呈指數(shù)次增長，只有這樣才能保證近似解的有效性；此外，隨機抽取混料試驗點樣本再進行比較的方法，效果不及全局搜索算法的結(jié)果，對于n個隨機混料試驗點，從中選取k點的設(shè)計，使用全局搜索需要運行的循環(huán)次數(shù)是，并且根據(jù)不同的信息陣需要編寫不同的字典序程序，而Monte-Carlo方法則是在n個隨機混料試驗點中隨機搜索，就像是“在茫茫數(shù)據(jù)中漫無目的地搜尋著最優(yōu)點?！彼赃@種方法仍存在很大的局限性.關(guān)于以上問題都有待進一步研究。