利用蒙特卡羅計算重積分的改進算法

王丙參，魏艷華，賈金平

（天水師范學院 數學與統計學院，甘肅 天水 741001）

0 引言

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

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

1 利用隨機數法計算重積分

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

方法1：均勻隨機數法

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

方法2：一般隨機數法

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

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

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

方法3：有利隨機數法

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

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

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

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

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

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

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

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

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

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

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

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

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

當直接從目標分布π(x)抽樣很難，而從包絡函數g(x)中抽樣和計算權重w(x)容易時，重點抽樣法有效，但是有效的高低很難估計，一個不太精確但有用的方法是：利用有效樣本的大小進行衡量。可解釋為：從試驗分布中抽取n個權重樣本相當于從目標函數中近似抽取EES(n)個獨立同分布的樣本。

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

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

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

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

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

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

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

這個方法聽起來很令人興奮，但是不切實際的，因為條件分布π(xt|x≤t-1)往往是得不到的，甚至得到條件分布比原問題還困難。為了實現序貫抽樣，引入更復雜的一些步驟。假定可找到一系列輔助分布π?(x≤t-1)來近似邊際分布π(x≤t-1)。需要強調的是，僅需知道除了歸一化常數外的 π?(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))，假如存在一個正確的權重集與隨機樣本相關，同時這些權重不是太不對稱，則樣本幾乎可從任何分布中抽樣。如果對于任意二次可積函數h(x)成立：

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

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

加權樣本的實質是為了強調對任意給定的樣本X存在多種加權w的方法，在重點抽樣中，權重w對應樣本X的一個確定函數，比如，這種情況下，(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）對函數f(x，y)做直觀分析后，決定選用試驗分布g(x，y)進行重點抽樣，其分布形式為：

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

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

落入區域χ外時，舍去。

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

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

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