張艷，楊金語

（1.天津財經大學珠江學院數據工程學院，天津 301811；2.南開大學統(tǒng)計與數據科學學院，天津 300071）

0 引言

在經濟、工程等領域的實際應用中，常需要估計某個指標的均值；在統(tǒng)計、數學等領域中，常需要估計某個復雜的定積分。然而，這些問題常常難以直接計算，此時可以借助蒙特卡羅方法模擬實現[1—3]，其主要思想是通過計算機軟件來模擬實踐中的隨機抽樣[2]。常用的蒙特卡羅方法有隨機投點法與平均值法，在相同的抽樣方法以及樣本容量（或試驗次數）下，平均值法得到的估計量的方差小于等于隨機投點法得到的估計量的方差[1]，即平均值法得到的估計量精度更高。以一維情形為例，簡述平均值法的主要過程如下：

設X～U(0,1)，且：

由式（1）知，可以用數學期望的估計值作為式（1）的積分的估計值。從(0,1)上的均勻分布中抽取獨立樣本x1,x2,…,xn，利用樣本均值估計μ的值，μ的估計值為：

其中，μ?是μ的無偏估計量，?的方差為：

當重復獲取m個樣本進行如上估計時，可得到估計量，則式（3）的方差可用式（4）進行估計：

為了獲得μ的精度較高的估計量，方差縮減技術得到了廣泛研究。常見的方差縮減技術主要包括對偶變量法、條件期望法、分層抽樣法、控制變量法、重要性抽樣法等[2]。這些方法的技巧主要分為以下四類：（1）只改變概率分布，比如重要性抽樣法是把概率分布改為重要概率分布，而分層抽樣是把整層概率分布改為分層概率分布；（2）只改變統(tǒng)計量，比如控制變量法與對偶變量法；（3）同時改變概率分布和統(tǒng)計量，比如條件期望法；（4）只改變統(tǒng)計特性，比如系統(tǒng)抽樣和擬蒙特卡羅方法是把隨機數改變?yōu)榉请S機數。這些方法均得到了廣泛的研究與應用[3—6]。事實上，單個技巧降低方差的效果往往有限，因此可以使用組合技巧[3]。

綜上所述，目前關于提高蒙特卡羅模擬估計精度的研究主要集中在方差縮減技術方面，而其中基于抽樣方法角度研究蒙特卡羅模擬估計精度的文獻較少。從式（2）至式（4）不難發(fā)現，樣本點x1,x2,…,xn的質量直接決定著估計量方差的大小，即決定著估計量的精度，而這n個樣本點質量的優(yōu)劣又受到抽樣方法的影響。此外，現有研究大多采用經典的抽樣方法，比如簡單隨機抽樣（Simple random sampling，SRS）、分層隨機抽樣（Stratified random sampling，STRS）、重要性抽樣等[1—7]，且將不同抽樣方法進一步結合起來的更少。本文將基于抽樣方法視角，借助試驗設計思想，將分層抽樣方法與拉丁超立方體抽樣方法[8]相結合，提出一種新的抽樣方法，以實現方差縮減，從而提高蒙特卡羅模擬的估計精度。

1 分層拉丁超立方體抽樣

1.1 拉丁超立方體抽樣的定義

在試驗設計中，一個設計可用一個矩陣表示，稱為設計陣。設計陣每一行對應一次試驗，每一列對應一個試驗因子，故設計陣的行數即為試驗次數，列數即為試驗包含因子的個數。拉丁超立方體抽樣（Latin Hypercube Sampling，LHS）是McKay等（1979）[8]提出的一種估計響應函數總體均值的抽樣方法，在計算機試驗中應用廣泛。與LHS對應的概念是拉丁超立方體設計（Latin Hypercube Design，LHD），LHS可由LHD變換而來。為此先介紹LHD的定義。

定義1：對于一個有n次試驗、q個因子的設計陣，若其每一列均是1,2,…,n的隨機排列，且各列是獨立生成的，則稱此設計為拉丁超立方體設計，記為LHD(n,q)。

由定義1可知，LHD中每個因子的水平數（因子取值的個數）均與試驗次數相同，這使得LHD能更好地適應計算機試驗不存在試驗誤差的特點，較大的水平數使其更具探索因子與響應之間復雜曲面關系的潛力。定義2說明了如何由一個LHD得到LHS。

定義2：給定一個LHD(n,q)，不妨記為A=(aij)n×q，令：

其中，uij是[0,1)上的均勻分布隨機數，則稱B=(bij)n×q為有q個因子、樣本容量為n的拉丁超立方體抽樣，記為LHS(n,q)。

由定義2知，第j個因子的第i個水平經過式（5）的變換，由{1,2,…,n}標準化到區(qū)間(0,1]內，從而由一個LHD(n,q)得到空間(0,1]q中的拉丁超立方體抽樣，即LHS(n,q)。相當于將空間(0,1]q等分成nq個小的超立方體，然后根據A=(aij)n×q的n行選出n個小的超立方體，再由第i行的q個均勻隨機數ui1,ui2,…,uiq決定樣本點在n個小的超立方體中的具體位置。這n個樣本點滿足如下性質：若把每一維區(qū)間(0,1]等分成n個小區(qū)間(0,1/n],(1/n,2/n],…,((n-1)/n,1]，則將n個樣本點向任一維度上投影時，每個小區(qū)間中包含且僅包含一個樣本點。第j個因子的第i個水平決定了在第j個維度上投影點落在哪個區(qū)間中，而隨機數uij決定了投影點的具體位置。下面給出一個例子，演示如何由一個LHD變換為相應的LHS。

例1：隨機生成一個包含2個因子、8次試驗的拉丁超立方體設計LHD(8,2)，并通過式（5）將其變成一個LHS(8,2)。

將此LHS(8,2)的樣本點標記在圖1中，其中(0,1]2被等分成8×8=64個小正方形。從圖1可知，LHS(8,2)的8個樣本點所在的小正方形由LHD(8,2)的行決定，而在小正方形中的具體位置則受隨機數uij（i=1,2,…,8;j=1,2）的影響。當把8個樣本點向兩個維度上投影時，每個維度上的8個等分區(qū)間(0,1/8],(1/8,2/8],…,(7/8,1]中均包含且只包含一個點。為了對比，隨機生成一個包含8個樣本點、2個變量的簡單隨機抽樣，記為SRS(8,2)，其散點圖如圖2所示。易見，在兩個維度上的8個等分區(qū)間(0,1/8],(1/8,2/8],…,(7/8,1]中，有的區(qū)間中沒有點，而有的區(qū)間中卻有多個點。顯然，在一維投影均勻性方面，LHS(n,q)要優(yōu)于SRS(n,q)。

圖1 LHS(8,2)的散點圖

圖2 SRS(8,2)的散點圖

1.2 分層拉丁超立方體抽樣的定義

不妨記有n個樣本點、s個層，每層樣本點數分別為t1,t2,…,ts的STRS為STRS(n,s;t1,t2,…,ts)，其中n=t1+t2+…+ts。下面結合STRS與LHS的定義，給出分層拉丁超立方體抽樣的定義。

定義3：在分層抽樣的基礎上，如果各層抽樣是按拉丁超立方體抽樣獨立地進行的，那么這種抽樣稱為分層拉丁超立方體抽樣（Stratified Latin Hypercube Sampling，STLHS）。記有n個樣本點、s個層，每層樣本點數分別為t1,t2,…,ts的STLHS為STLHS(n,s;t1,t2,…,ts)，其中n=t1+t2+…+ts。

由定義3可知，STLHS與STRS的區(qū)別在于每層中的抽樣方式不同。對于STRS，每層中的抽樣是按簡單隨機抽樣進行的；而對于STLHS，每層中的抽樣是按LHS進行的。下面給出STLHS(n,s;t1,t2,…,ts)的抽樣算法。

1.3 分層拉丁超立方體抽樣的算法

本文將給出一維STLHS(n,s;t1,t2,…,ts)的兩種算法。

（1）算法1

步驟1：給定分層個數s，將區(qū)間(0,1]等分為s個小區(qū)間(0,1/s],(1/s,2/s],…,((s-1)/s,1]。

步驟2：將第i層((i-1)/s,i/s]等分為ti個小區(qū)間，i=1,2,…,s，即：

步驟3：在式（6）中的每個小區(qū)間中，隨機抽取一個樣本點，得到第i層的ti個樣本點，并將ti個樣本點的順序隨機化。

步驟4：將各層樣本點組合得到最終的總樣本。

由算法1的步驟2以及定義2可知，第i層產生的樣本點組成一個LHS(ti,1)，i=1,2,…,s。再由定義3可知，所得抽樣結果為一個分層拉丁超立方體抽樣，具體為STLHS(n,s;t1,t2,…,ts)。

（2）算法2

步驟1：給定分層個數s，將區(qū)間(0,1]等分為s個小區(qū)間(0,1/s],(1/s,2/s],…,((s-1)/s,1]。

步驟2：對于第i層，隨機生成一個LHS(ti,1)，i=1,2,…,s。

步驟3：將第i層的LHS(ti,1)按式（7）進行變換：

步驟4：將各層樣本點組合得到最終的總樣本。

基于LHS，算法2更易實現，但算法2和算法1所得的抽樣結果等價。

下面用例2對以上算法的結果進行展示。

例2：考慮隨機生成一個STLHS(10,2;4,6)，記為A，易見s=2,t1=4,t2=6。

先利用算法1生成。根據步驟1，將區(qū)間(0,1]等分為2個區(qū)間(0,1/2],(1/2,1]；由步驟2，將(0,1/2]等分為4個區(qū)間(0,0.125]，(0.125,0.250]，(0.250,0.375]，(0.375,0.500]，將(1/2,1]等分為6個區(qū)間(0.500,0.583]，(0.583,0.667],(0.667,0.750]，(0.750,0.833],(0.833,0.917]，(0.917,1]；根據步驟3，在以上10個小區(qū)間中分別抽取一個樣本點；最后，由步驟4將所有樣本點合在一起即可得到STLHS(10,2;4,6)。

再根據算法2生成。根據步驟1，將區(qū)間(0,1]等分為2個區(qū)間(0,1/2],(1/2,1]；根據步驟2，在兩層區(qū)間(0,1/2]和(1/2,1]中分別生成LHS(4,1)和LHS(6,1)；根據步驟3，將LHS(4,1)變換為[LHS(4,1)+1-1)]/2，將LHS(6,1)變換為[LHS(6,1)+2-1)]/2；最后，由步驟4將所有樣本點合在一起即可得到STLHS(10,2;4,6)。

為了展示STLHS的均勻性，現隨機生成一個STRS(10,2;4,6)，記為B。兩種方法下的抽樣結果A與B如下所示：

為了更清晰地了解兩種抽樣的區(qū)別，現繪制A與B的散點圖，分別如圖3與圖4所示。

圖3 STLHS(10,2;4,6)樣本點的分布圖

圖4 STRS(10,2;4,6)樣本點的分布圖

從圖3可以看出，對于本文新提出的分層拉丁超立方體抽樣STLHS(10,2;4,6)來說，第一層(0,0.5]的4個樣本點分別分布在4個等距區(qū)間(0,0.125],(0.125,0.25],(0.25,0.375],(0.375,0.5]中，而第二層(0.5,1]的6個樣本點分別分布在(0.5,0.583],(0.583,0.667],(0.667,0.75],(0.75,0.833],(0.833,0.917],(0.917,1]中，即在每一層中都實現最大程度的分層均勻性。從圖4可以看出，對于分層隨機抽樣STRS(10,2;4,6)來說，在某些小區(qū)間中沒有樣本點（比如(0,0.125]），而在某些小區(qū)間中有多個樣本點（比如(0.375,0.500]中有2個，(0.750,0.833]中有3個），即每層的樣本點都無法保證達到最大程度的分層均勻性。

1.4 分層拉丁超立方體抽樣的性質

為討論抽樣方法的性質，給出如下抽樣方法優(yōu)劣的定義。

定義4：若在抽樣方法甲下比在抽樣方法乙下得到的μ的估計量方差更小，則稱抽樣方法甲優(yōu)于抽樣方法乙。

比如，STRS優(yōu)于SRS[9]，即與SRS相比，基于STRS得到的μ的估計量的方差更小。本文將從理論上闡述STLHS相較于STRS在估計μ時所帶來的改進。

假設分層隨機抽樣中的層數與分層拉丁超立方體抽樣中的層數是相同的，記為s，設第c層積分值為μc，且利用式（2）得到的μc的估計量為?c，則μ的估計量為：

其中，0<λc<1為權重，且=1。關于權重λc的選擇，可采用等比例分配、最優(yōu)分配、奈曼最優(yōu)分配等原則確定。特別地，如果將每層樣本看作同等重要，那么也可以令λ1=λ2=…=λs=1/s。若分別采用STRS和STLHS兩種方法獲取樣本點來估計μ的值，有以下結論：

定理1：記?STLHS為采用STLHS方法獲取樣本點對μ的估計量，則?STLHS是μ的無偏估計量。

由于利用分層抽樣得到的μ的估計量是無偏的[9]，因此定理1的證明是顯然的。為了給出進一步的結論，先介紹引理[10]。

引理1：若函數f()x在每一維變量xj（j=1,2,…,q）上均具有單調性，則LHS優(yōu)于SRS，即VLHS(?)≤VSRS(?)，其中，VLHS(?)與VSRS(?)分別表示在LHS與SRS下得到的μ?的方差。

在引理1的基礎上，有如下結論：

定理2：若每一層中對應的響應函數fc(x)(c=1,2,…,s）關于x均具有單調性，則利用式（8）估計μ時，STLHS優(yōu)于STRS。

證明：根據STRS與STLHS的定義可知，STRS中每一層的樣本是采用SRS得到的，而STLHS中每一層的抽樣方法屬于LHS。因此，基于引理1的結論，在定理2的條件下，利用第c層中的樣本估計時，LHS優(yōu)于SRS，即：

由式（8）可知：

即STLHS優(yōu)于STRS。

2 數值模擬

本文通過一個例子對定理2的結論進行數值模擬驗證。為了比較方便，不妨假設各層具有相同的樣本量n，且λ1=λ2=…=λs=1/s。

例3：設f(x)=lnx為定義在區(qū)間(0,1]上的響應函數，X為在(0,1]上均勻分布的隨機變量，計算f(x)在(0,]1的積分值為：

下面分別采用STLHS與STRS兩種抽樣方法獲取樣本對μ進行估計，以驗證定理2的結論。為了簡便，只比較層數s為2，每層樣本量n均為5、10、20與50的四種組合下的情況。在每種組合下，重復抽樣10000次，每次得到樣本后均利用式（8）計算μ的估計量，其中，,模擬結果在表1中給出。

表1 兩種抽樣方法的估計量及其估計方差

從表1不難看出，在四種參數組合下，STLHS下得到的估計方差總是小于STRS下得到的估計方差，即STLHS優(yōu)于STRS，這就驗證了定理2的結論。此外，在兩種抽樣方法下，隨著樣本量的增加，估計方差都有下降的趨勢，這也與常識相符。

3 結束語

本文基于拉丁超立方體抽樣在估計總體均值時優(yōu)于簡單隨機抽樣的事實，對分層隨機抽樣進行了改進，提出了分層拉丁超立方體抽樣的概念，給出了兩種抽樣算法。此外，不僅從理論上證明了在估計總體均值方面，若被積函數滿足單調性時，則分層拉丁超立方體抽樣優(yōu)于分層隨機抽樣，還利用數值模擬對理論結果進行了驗證。本文的結果為蒙特卡羅方法估計高難度積分提供了很好的抽樣方法。值得指出的是，本文僅針對一維情形給出了抽樣算法，即只有一個分層變量情形。對于包含兩個或兩個以上分層變量的情形，若每一維的分層變量是獨立的，則可以基于本文的算法1或算法2在每一維進行獨立抽樣，從而獲得所需樣本；而對于分層變量不獨立的情形，并不能基于一維結果進行簡單推廣，這也是接下來值得研究的問題。此外，式（8）中權重λc的選擇也是值得探索的問題。