王丙參，魏艷華，賈金平

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

0 引言

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

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

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

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

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

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

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

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

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

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

方法3：有利隨機(jī)數(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ù)增加，它的計(jì)算效率顯著降低，而蒙特卡羅方法的精度與積分維數(shù)無關(guān)。下面利用蒙特卡羅方法計(jì)算積分I。

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

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

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

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

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

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

2 利用重點(diǎn)抽樣技術(shù)改進(jìn)重積分的蒙特卡羅計(jì)算

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

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

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

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

當(dāng)直接從目標(biāo)分布π(x)抽樣很難，而從包絡(luò)函數(shù)g(x)中抽樣和計(jì)算權(quán)重w(x)容易時(shí)，重點(diǎn)抽樣法有效，但是有效的高低很難估計(jì)，一個(gè)不太精確但有用的方法是：利用有效樣本的大小進(jìn)行衡量?？山忉尀椋簭脑囼?yàn)分布中抽取n個(gè)權(quán)重樣本相當(dāng)于從目標(biāo)函數(shù)中近似抽取EES(n)個(gè)獨(dú)立同分布的樣本。

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

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

這樣，未標(biāo)準(zhǔn)化權(quán)重可寫為：

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

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

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

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

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

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

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

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

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

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

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

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

所以該積分的估計(jì)為：

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

圖2原始蒙特卡羅（左）與重點(diǎn)抽樣（右）的100次模擬結(jié)果直方圖

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

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

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

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

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

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

顯然，利用重點(diǎn)抽樣技術(shù)，蒙特卡羅的標(biāo)準(zhǔn)誤差減少了很多，大大提高了抽樣效率，這是因?yàn)槭济商乜_存在與確定性方法同樣的問題：將大量的時(shí)間浪費(fèi)在計(jì)算那些函數(shù)值幾乎為0的隨機(jī)樣本上。