王丙參，魏艷華，賈金平

（天水師范學院 數(shù)學與統(tǒng)計學院，甘肅 天水 741001）

0 引言

蒙特卡羅方法（M-C，即隨機模擬）是利用計算機生成隨機數(shù)并通過模擬隨機現(xiàn)象來進行近似計算的方法。在工程、科學等領域中，實驗數(shù)據(jù)很難獲得或成本太高，這時蒙特卡羅方法就是最經(jīng)濟、實用的方法[1-4]。另外，對于一些復雜問題，比如方程組求解、最優(yōu)化、偏微分方程的解等，蒙特卡羅方法也非常有效，但它最常用的還是計算積分，是利用數(shù)值方法計算多重積分的首選方法，且積分維數(shù)越高，效率也越高。一般情況下，在蒙特卡羅算法中，利用均勻隨機數(shù)計算積分簡單，但精度不高，誤差與被積函數(shù)方差和模擬次數(shù)有關。通過增加模擬次數(shù)提高精度的速度太慢，且耗時過多，尤其對于多重積分問題更加嚴重。如果改進蒙特卡羅方法，可以縮減積分估計的方差，那就能提高估計的精度和可靠性[5-8]。

本文比較研究了計算重積分的隨機數(shù)法，并利用重點抽樣技術改進重積分的蒙特卡羅計算，提高了抽樣效率，得出有利隨機數(shù)法是重點抽樣法的特例，誤差較小。

1 利用隨機數(shù)法計算重積分

很多物理問題都可轉化為多重積分的計算問題，在實際計算中，經(jīng)常發(fā)現(xiàn)被積函數(shù)的原函數(shù)很難求出，或原函數(shù)不是初等函數(shù)，對于這樣問題，最好設計一種蒙特卡羅方法進行近似計算[7，8]。

方法1：均勻隨機數(shù)法

（1）取包含D（積分區(qū)域）的矩形區(qū)域Ω：a≤x≤b，c≤y≤d，其面積SΩ=(b-a)(d-c)。（2）生成 Ω 上均勻隨機數(shù)列 (xi，yi)，i=1，2，…，n，不妨設前k個隨機數(shù)落入積分區(qū)域D，則當n充分大時，

方法2：一般隨機數(shù)法

（1）取包含D（積分區(qū)域）的矩形區(qū)域Ω：a≤x≤b，c≤y≤d，其面積SΩ=(b-a)(d-c)。

（2）取概率密度函數(shù)g(x，y)，使得

（3）生成隨機數(shù)列 (xi，yi)～g(x，y)，i=1，2，…，n，不妨設前k個隨機數(shù)落入積分區(qū)域D，則當n充分大時，

方法3：有利隨機數(shù)法

由于f(x，y)=ln(1+2x+2y)的原函數(shù)為所以積分I的精確解為F(1，1)-F(1，0)-F(0，1)+F(0，0)=1.057615826853317。

利用 MATLAB 內(nèi)置函數(shù) dblquad(′log(1+2*x+2*y)′,0,1,0,1)可得數(shù)值積分為1.057615735740697。雖然數(shù)值積分的精度已經(jīng)很高，但隨著維數(shù)增加，它的計算效率顯著降低，而蒙特卡羅方法的精度與積分維數(shù)無關。下面利用蒙特卡羅方法計算積分I。

對f(x，y)在進行泰勒展開，可得f(x，y)=ln(1容易估算出I≈ln3，故取有利概率密度函數(shù)

程序中取100個隨機數(shù)，做了1000次模擬，其中，平均誤差等于模擬值減去真實值的絕對值的平均值，模擬結果如下：

圖1蒙特卡羅積分值的直方圖（左：均勻隨機數(shù)法，右：有利隨機數(shù)法）

表1 1000次模擬的均值、方差和平均誤差

由圖1和表1可知，不論用均勻隨機數(shù)序列還是有利隨機數(shù)列，蒙特卡羅積分值的分布都近似服從正態(tài)分布，分別為：

均勻隨機數(shù)法的方差與平均誤差明顯大于有利隨機數(shù)法的方差與平均誤差，這說明有利隨機數(shù)法的模擬結果更靠近真值（1.057615826853317），即有利隨機數(shù)法的誤差更小，更可靠。

2 利用重點抽樣技術改進重積分的蒙特卡羅計算

假設包絡函數(shù)g(x)在分布族{g(x;λ)}中取得，其中參數(shù)λ為限制條件，則目的就是在給定限制λ的條件下，求解：

其中，supp(g)表示g的支撐。

即Jenson不等式等號成立的充要條件?？梢?，當x(i)～g∝|h|π(x)時，有：

由于Eπ[h(X)]為常數(shù)，故最優(yōu)包絡函數(shù)g*不同于密度函數(shù)π。如果π(x)?g(x)，其中x為來自g(x)的抽樣，則權重會變得很大，從而估計方差也會變得很大，因此選取合適的包絡函數(shù)g(x)才能降低估計的方差。要使得估計方差足夠小，試驗分布g(x)的形狀要盡可能接π(x)h(x)，特別地，g應有一個比π(x)更長的尾部，并可以很方便從g(x)中抽取樣本。盡管有時候無法或很難在高維空間找到一個恰當?shù)脑囼灧植?，但這卻是重點抽樣的主要任務。

當直接從目標分布π(x)抽樣很難，而從包絡函數(shù)g(x)中抽樣和計算權重w(x)容易時，重點抽樣法有效，但是有效的高低很難估計，一個不太精確但有用的方法是：利用有效樣本的大小進行衡量?？山忉尀椋簭脑囼灧植贾谐槿個權重樣本相當于從目標函數(shù)中近似抽取EES(n)個獨立同分布的樣本。

在高維重點抽樣中尋找一個好的試驗分布是很困難的。一個解決辦法就是循序地構造試驗分布函數(shù)。假定x=(x1，…，xp)，其中每個xi仍可為多維隨機向量，則試驗分布可構造如下：

其中x≤p=(x1，x2，…，xp)。同樣，根據(jù)x的分塊，給出目標分布的分解：

這樣，未標準化權重可寫為：

該式建議從g1(x1)，g2(x2|x≤1)，…，gp(xp|x≤p-1)中序貫抽取X的各分量。令，則有遞推關系：

顯然，wp(x≤p)=w(x)。這個方法具有如下優(yōu)點：

（1）如果序貫產(chǎn)生的部分樣本計算得到的權重值太小，則可停止進一步產(chǎn)生x的分量。

（2）可利用 π(xt|x≤t-1)來設計g(xt|x≤t-1)，即有目標分布π(x)的邊際分布指導x的產(chǎn)生。

這個方法聽起來很令人興奮，但是不切實際的，因為條件分布π(xt|x≤t-1)往往是得不到的，甚至得到條件分布比原問題還困難。為了實現(xiàn)序貫抽樣，引入更復雜的一些步驟。假定可找到一系列輔助分布π?(x≤t-1)來近似邊際分布π(x≤t-1)。需要強調(diào)的是，僅需知道除了歸一化常數(shù)外的 π?(x≤t-1)，并且在構造整個樣本x=(x1，…，xp)時僅起到指導作用，而不一定需要從中抽樣。這樣，可用如下的遞推過程來進行序貫重點抽樣（SIS）：

（1）從g(xt|x≤t-1)中出去樣本xt，并令x≤t=(x≤t-1，xt)。

在序貫重點抽樣中，ut稱為增量權。使用序貫重點抽樣的原因是：利用輔助分布π?(x≤t-1)可構造更有效的試驗分布，并將一個極其困難的任務分解為多個易處理的小任務，從而增加了可行性。

重點抽樣利用不是來自π(x)中的隨機樣本X(1)，…，X(n)估計目標μ=Eπ(h(X))，假如存在一個正確的權重集與隨機樣本相關，同時這些權重不是太不對稱，則樣本幾乎可從任何分布中抽樣。如果對于任意二次可積函數(shù)h(x)成立：

其中，a為全部n個樣本的共同歸一化常數(shù)，則稱加權隨機樣本集關于 π 是正確的。

容易證明，如果wt-1(x≤t-1)是x≤t-1關于 π≤t-1的正確權，則wt(x≤t)是x≤t關于 π≤t的正確權。因此，用最后得到的正確權wp對序貫方式得到的目標分布π(x)的整個樣本x賦予了正確權。

加權樣本的實質(zhì)是為了強調(diào)對任意給定的樣本X存在多種加權w的方法，在重點抽樣中，權重w對應樣本X的一個確定函數(shù)，比如，這種情況下，(w，X)s得加權變量w為非退化隨機變量。

（1）原始蒙特卡羅算法：生成n對隨機樣本(x(1)，y(1))，…，(x(n)，y(n))，(x(i)，y(i))～π ，其中 π 服從[-1，1]×[-1，1]上的均勻分布，因為：

所以該積分的估計為：

在每次模擬中，令n=5000，共模擬m=100次,模擬結果見圖2左。

圖2原始蒙特卡羅（左）與重點抽樣（右）的100次模擬結果直方圖

（2）對函數(shù)f(x，y)做直觀分析后，決定選用試驗分布g(x，y)進行重點抽樣，其分布形式為：

其中 (x，y)∈χ=[-1，1]×[-1，1]。這是一截尾混合高斯分布：利用如下步驟從這一混合高斯分布中抽樣：

生成均勻隨機數(shù)u～U(0，1)，如果u≤0.464，則生成隨機數(shù)否則，生成隨機數(shù)

落入?yún)^(qū)域χ外時，舍去。

在每次模擬中，令n=5000，共模擬m=100次,模擬結果見圖2右。

θ?=0.1259，std(θ?)=3.7462e-004

顯然，利用重點抽樣技術，蒙特卡羅的標準誤差減少了很多，大大提高了抽樣效率，這是因為始蒙特卡羅存在與確定性方法同樣的問題：將大量的時間浪費在計算那些函數(shù)值幾乎為0的隨機樣本上。