具有下界約束混料試驗的Monte-Carlo漸近最優設計

楊曉珍

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

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

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

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

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

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

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

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

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

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

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

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

在D-準則下，關于信息矩陣的判定函數為:

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

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

在MV-準則下，關于信息矩陣的判定函數為:

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

二、Monte-Carlo方法構造試驗點樣本

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

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

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

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

則兩種坐標間的變換關系為：

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

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

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

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

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

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中各元素的聯合密度函數可以表示為：

其中I{?}為示性函數。

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

三、模擬近似最優設計

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

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

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

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

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

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

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

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

表1 飽和近似D最優設計

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

圖2 飽和近似D最優設計的響應曲面與等高圖

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

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

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

表2 含測度的7點近似D最優設計

表3 飽和近似MV最優設計及試驗結果

下面考慮在MV-準則下的近似最優解，此時關于信息矩陣的判定函數為:

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

并按升序排列為:

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

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

通過計算得到：

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

圖3 飽和近似MV最優設計的響應曲面與等高圖

四、討論

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

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