逆變換Monte Carlo 算法及其應用

2022-02-21 01:55:42孟憲勇

科學技術創新 2022年3期

孟憲勇王晶劉彭

(山東農業大學信息科學與工程學院,山東泰安 271038)

在應用數學和許多工程學科中,計算函數f(x)在某一范圍s 內的積分,即

是普遍的事情,但這個積分一般不能用分析的方法解決或計算是復雜的,此時只能考慮近似算法,而Monte Carlo 就是大家經常利用的方法之一。

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

于是求積分問題轉化為如何選擇隨機模型p(x)及由此產生隨機數。

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

1 逆變換Monte Carlo 算法

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

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

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

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

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

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

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

以上定理告訴我們如何從分布函數為F(x)的隨機模型進行計算機采樣。首先產生模型U(0,1)的隨機數u,然后通過函數F-1(u)計算即可,這樣的采樣方法就是逆變換采樣法。

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

實際抽樣過程中,可細分為連續型分布和離散型分布抽樣法。對離散型隨機變量,定理1 可以敘述如下:

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

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

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

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

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

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

(1)產生U(0,1)隨機數u；

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

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

2 逆變換抽樣技術程序實現

下面舉例說明逆變換技術程序實現及其在分段模型中的抽樣效果。我們使用的軟件是R 軟件。R 是用于計算和制圖的優秀免費統計軟件[9],非常有利于統計思想和統計方法的學習。

例1 設隨機變量X 的分布律為

試產生X 的隨機數。

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

則數據x1,x2,…,xn是隨機變量X 的隨機數。

根據以上過程,產生2000 個隨機數的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)

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