逆變換Monte Carlo 算法及其應(yīng)用

孟憲勇 王 晶 劉 彭

(山東農(nóng)業(yè)大學(xué)信息科學(xué)與工程學(xué)院,山東 泰安 271038)

在應(yīng)用數(shù)學(xué)和許多工程學(xué)科中,計(jì)算函數(shù)f(x)在某一范圍s 內(nèi)的積分,即

是普遍的事情,但這個(gè)積分一般不能用分析的方法解決或計(jì)算是復(fù)雜的,此時(shí)只能考慮近似算法,而Monte Carlo 就是大家經(jīng)常利用的方法之一。

Monte Carlo 算法,也稱為隨機(jī)模擬、隨機(jī)仿真或采樣算法,是一種只使用計(jì)算機(jī)對模型進(jìn)行采樣得到隨機(jī)數(shù),然后基于隨機(jī)數(shù)進(jìn)行計(jì)算的數(shù)學(xué)方法,也就是把現(xiàn)實(shí)的數(shù)據(jù)采集搬到計(jì)算機(jī)上進(jìn)行。算法由威勒蒙和馮紐曼在20 世紀(jì)40 年代為適應(yīng)原子彈事業(yè)的發(fā)展而首先提出,被稱為20 世紀(jì)十大最偉大算法之一。時(shí)至今日,算法已在各類文獻(xiàn)中司空見慣[1-2],數(shù)學(xué)中許多著名的算法與技巧其實(shí)就是采樣算法,這充分說明:今天的我們應(yīng)該學(xué)習(xí)和掌握這種算法。

于是求積分問題轉(zhuǎn)化為如何選擇隨機(jī)模型p(x)及由此產(chǎn)生隨機(jī)數(shù)。

隨機(jī)數(shù)的產(chǎn)生是Monte Carlo 算法的核心步驟,大多數(shù)情況下不是容易的事,其生成要遵循一定的方法,一般可分為:MCMC 法、篩選法,逆變換法等?？偨Y(jié)已有文獻(xiàn)不難看出,研究采樣的文獻(xiàn)大多著墨于前兩種方法,而對逆變換采樣往往一語帶過,缺乏系統(tǒng)而易理解的介紹,很不利于算法的學(xué)習(xí)和應(yīng)用,具體可參見文獻(xiàn)[3-7]。逆變換采樣法簡單直觀高效,思想深刻,是MCMC 方法的必要基礎(chǔ),對正確理解和掌握Monte Carlo 算法,提高實(shí)際應(yīng)用能力具有不可替代的作用。因此,借助R 軟件[8],全面而系統(tǒng)的探討逆變換Monte Carlo 算法及其應(yīng)用,以便于大家對此算法的學(xué)習(xí)是本文唯一的出發(fā)點(diǎn)。

1 逆變換Monte Carlo 算法

設(shè)X 是一隨機(jī)變量,其值域?yàn)锳?R,則X 的隨機(jī)性質(zhì)可用分布函數(shù)F(x)或密度函數(shù)p(x)精確刻畫,其中分布函數(shù)定義為F(x) =P{X≤x}。顯然,分布函數(shù)是隨機(jī)變量值的累積概率,它的自變量取值是實(shí)數(shù),值域是區(qū)間[0,1]。

我們不難想象如果在區(qū)間(0,1)內(nèi)均勻采樣,然后將得到的隨機(jī)數(shù)根據(jù)分布函數(shù)F(x)反映射到對應(yīng)的自變量點(diǎn)就可以實(shí)現(xiàn)對隨機(jī)變量的采樣,且這樣的采樣能充分反映隨機(jī)信息,即在概率增長快的地方多采樣,而在增長慢的地方少采樣。

然而,令人遺憾的是F(x)不存在一般意義上的反函數(shù),為此給出如下定義:

F-1(y) = inf{x:F(x) ≥y，x∈A} (0 ≤y≤1)

定理1 設(shè)隨機(jī)變量X 的分布函數(shù)為F(x),U～U(0,1),則隨機(jī)變量Y=F-1(U)的分布函數(shù)為F(y)。

證明 由F-1(U)的定義和均勻分布可得

P{Y≤y} =P{F-1(U) ≤y} =P{U≤F(y)} =F(y)

以上定理告訴我們?nèi)绾螐姆植己瘮?shù)為F(x)的隨機(jī)模型進(jìn)行計(jì)算機(jī)采樣。首先產(chǎn)生模型U(0,1)的隨機(jī)數(shù)u,然后通過函數(shù)F-1(u)計(jì)算即可,這樣的采樣方法就是逆變換采樣法。

逆變換采樣法是一種通用算法,但其缺點(diǎn)是非常明顯的,算法必須容易求得逆函數(shù),這對許多函數(shù)是不容易做到的。因此,在隨機(jī)數(shù)生成中是很少應(yīng)用的,但對某些分布還是非常實(shí)用,比如一般的均勻分布、離散分布、密率分布等。

實(shí)際抽樣過程中,可細(xì)分為連續(xù)型分布和離散型分布抽樣法。對離散型隨機(jī)變量,定理1 可以敘述如下:

定理2 清楚表明了離散隨機(jī)變量的逆變換抽樣過程,先產(chǎn)生(0,1)上的均勻隨機(jī)數(shù)u,再根據(jù)u 的數(shù)值確定隨機(jī)數(shù)。

對連續(xù)型隨機(jī)變量,其分布函數(shù)的逆變換就是一般意義下的反函數(shù)。定理1 可敘述如下:

定理3 設(shè)X 是連續(xù)型隨機(jī)變量,取值于區(qū)間(a,b),分布函數(shù)為F(x),U～U(0,1),則Y=F-1(U)～F(y)。

連續(xù)型隨機(jī)變量的分布函數(shù)的反函數(shù)不一定容易求出,但對一些特殊形式的分布,還是容易求的,如密度函數(shù)是分段常數(shù)的分布。設(shè)X 的密度函數(shù)為

F(x) =p+c(x-x)

令F(X)=U,則U～U(0,1),解得X=xj+(U-pj)/cj,此處j滿足pj≤U＜pj+1。因此,產(chǎn)生密度函數(shù)為分段水平的抽樣數(shù)的具體算法為:

(1)產(chǎn)生U(0,1)隨機(jī)數(shù)u；

(3)計(jì)算x=xj+(u-pj)/cj,則x 是隨機(jī)變量X 的隨機(jī)數(shù)。

U(0,1)的隨機(jī)數(shù)生成是隨機(jī)數(shù)生成的基礎(chǔ),非均勻分布隨機(jī)數(shù)往往是由U(0,1)上隨機(jī)數(shù)經(jīng)過運(yùn)算而產(chǎn)生。今天大多數(shù)計(jì)算機(jī)語言(軟件)中都有自帶的程序可方便用于從U(0,1)生成隨機(jī)數(shù),大家只要掌握這些函數(shù)即可,不必細(xì)究其產(chǎn)生過程。當(dāng)然,感興趣者可參閱文獻(xiàn)[5-7]。

2 逆變換抽樣技術(shù)程序?qū)崿F(xiàn)

下面舉例說明逆變換技術(shù)程序?qū)崿F(xiàn)及其在分段模型中的抽樣效果。我們使用的軟件是R 軟件。R 是用于計(jì)算和制圖的優(yōu)秀免費(fèi)統(tǒng)計(jì)軟件[9],非常有利于統(tǒng)計(jì)思想和統(tǒng)計(jì)方法的學(xué)習(xí)。

例1 設(shè)隨機(jī)變量X 的分布律為

試產(chǎn)生X 的隨機(jī)數(shù)。

解 設(shè)u1,u2,…,un是(0,1)上均勻分布的隨機(jī)數(shù),令

則數(shù)據(jù)x1,x2,…,xn是隨機(jī)變量X 的隨機(jī)數(shù)。

根據(jù)以上過程,產(chǎn)生2000 個(gè)隨機(jī)數(shù)的R 程序如下:

u＜-runif(2000)

x＜-rep(0,length(u))

for(i in 1:length(u)){

if(u[i] 0.2) x[i]＜-0

else if(u[i] 0.7) x[i]＜-1

else x[i]＜-2}

hist(x)

產(chǎn)生2000 個(gè)隨機(jī)數(shù)的直方圖如圖1 左圖所示,該圖顯示隨機(jī)數(shù)分布與理論分布已十分接近,說明算法合理。事實(shí)上,R 軟件中的抽樣函數(shù)sample()也可以輕松完成取值有限的離散模型的抽樣,程序d＜-sample(0:2,2000,prob=c(0.2,0.5,0.3))即可完成抽樣任務(wù),其直方圖如圖1 右圖所示。

圖1 兩種方式產(chǎn)生離散型隨機(jī)變量的隨機(jī)數(shù)

例2 設(shè)連續(xù)型隨機(jī)變量X 的密度函數(shù)為

試生成X 的2000 個(gè)隨機(jī)數(shù)

解 由前面的分析,可給出R 程序代碼如下:

u＜-runif(2000)

x＜-rep(0,length(u))

for(i in 1:length(u)){

if(u[i]＜0.2) x[i]＜-u[i]/0.2

else if(u[i]＜0.7) x[i]＜-1+(u[i]-0.2)/0.5

else x[i]＜-2+(u[i]-0.7)/0.3}

hist(x,freq=F)

產(chǎn)生2000 隨機(jī)數(shù)的直方圖如圖2 所示,有圖不難看出隨機(jī)數(shù)分布規(guī)律與理論分布十分接近,說明隨機(jī)數(shù)產(chǎn)生方法合理。

圖2 產(chǎn)生2000 個(gè)連續(xù)分布隨機(jī)數(shù)的直方圖

3 逆變換抽樣算法的實(shí)際應(yīng)用

積分運(yùn)算是實(shí)際研究中非常基本的運(yùn)算,有些積分運(yùn)算是不能通過解析求解,只能得到近似解。Monte Carlo 算法可輕松用于積分近似求解。Monte Carlo 算法許多應(yīng)用的根本原因在于求解積分運(yùn)算。

由前面的分析知:

因此,積分的采樣近似算法可轉(zhuǎn)化為如下3 步工作:

(1)選取密度函數(shù)P(x)；

(2)抽取隨機(jī)數(shù)x1,…,xn；

解 我們分別選取密度函數(shù)為(0,3)上的均勻分布p1(x)和分段常數(shù)密度函數(shù)p2(x)為抽樣分布,以說明分段常數(shù)密度函數(shù)的優(yōu)點(diǎn),p2(x)定義如下。

表1 分別列出了抽樣數(shù)分別是10,100,1000,10000,100000 時(shí)的計(jì)算結(jié)果,其精確結(jié)果為1444.53,兩個(gè)密度函數(shù)均可用于積分運(yùn)算,p2(x)的抽樣優(yōu)于p1(x)。

表1

4 結(jié)論

本文主要對Monte Carlo 算法中的一類通用算法,即逆變換Monte Carlo 算法進(jìn)行了分析討論。算法理論簡單明了,易于掌握和理解。從實(shí)際計(jì)算來說,盡管有時(shí)不能應(yīng)用,但對一些特殊模型采樣還是非常高效的。

積分計(jì)算幾乎是所有學(xué)科的基礎(chǔ)工作,許多計(jì)算不得不依賴于Monte Carlo 算法,這是Monte Carlo 算法應(yīng)用廣泛的根源。本文提出了利用分段常數(shù)模型的逆變換采樣進(jìn)行具體積分的計(jì)算,算法高效且易于用R 軟件實(shí)現(xiàn)。

逆變換Monte Carlo 算法的缺點(diǎn)也一目了然,即必須容易確定分布函數(shù)的反函數(shù),這必然影響算法的實(shí)際應(yīng)用。更好的計(jì)算機(jī)抽樣技術(shù)MCMC 和篩選法將應(yīng)運(yùn)而生,我們會在今后繼續(xù)研究探討這兩種方法。