摘 要: MCMC 算法是一種基于隨機(jī)采樣的統(tǒng)計(jì)計(jì)算方法, 通過構(gòu)建馬爾可夫鏈, 從而實(shí)現(xiàn)對(duì)復(fù)雜系統(tǒng)的抽樣和模擬。文章基于蒙特卡洛馬爾可夫鏈的相關(guān)原理, 重點(diǎn)分析蒙特卡洛馬爾可夫鏈在數(shù)值模擬中的應(yīng)用。首先, 介紹了MCMC 算法及相關(guān)算法的原理, 包括蒙特卡洛方法的原理和三種常用的采樣方法—直接采樣、接受拒絕采樣和重要性采樣、馬爾可夫鏈的原理和平穩(wěn)狀態(tài)、MCMC 算法的原理及兩種經(jīng)典的采樣算法—Metropolis Hastings 算法和Gibbs 采樣; 其次, 著重研究MCMC算法的應(yīng)用, 包括蒙特卡洛方法在擬合參數(shù)、求解不規(guī)則面積以及接受拒絕采樣的應(yīng)用, 馬爾可夫鏈在租車還車問題中的應(yīng)用以及使用MCMC 算法估計(jì)回歸系數(shù)的應(yīng)用; 最后, 文章指出MCMC 算法的改進(jìn)方向, 促進(jìn)MCMC 算法更廣泛的應(yīng)用。
關(guān)鍵詞: 數(shù)值模擬; 蒙特卡洛方法; 馬爾可夫鏈; MCMC 算法
基金項(xiàng)目: 山西省高等學(xué)校教學(xué)改革創(chuàng)新項(xiàng)目“轉(zhuǎn)型發(fā)展背景下地方高校理工科專業(yè)《概率論與數(shù)理統(tǒng)計(jì)》教學(xué)改革與實(shí)踐” (J20221082); 長(zhǎng)治學(xué)院教學(xué)改革創(chuàng)新項(xiàng)目“新工科背景下《概率論與數(shù)理統(tǒng)計(jì)》課程的實(shí)踐教學(xué)研究” (JC202310)
中圖分類號(hào): O212 文獻(xiàn)標(biāo)識(shí)碼: A
文章編號(hào): 1674-537X (2024) 10. 0004-10
一、引言
隨信息時(shí)代的發(fā)展, 數(shù)值模擬技術(shù)在越來越多的領(lǐng)域被廣泛應(yīng)用, 然而, 隨著問題的復(fù)雜度和數(shù)據(jù)量的增加, 傳統(tǒng)的數(shù)值模擬方法往往面臨著計(jì)算效率低下、模型參數(shù)估計(jì)困難等挑戰(zhàn)。
作為基于隨機(jī)抽樣的統(tǒng)計(jì)計(jì)算方法, MCMC 算法可以有效的解決上述問題, 被廣泛應(yīng)用于從復(fù)雜概率分布中抽樣。在數(shù)值模擬的過程中, MCMC 不僅可以為數(shù)值模擬提供抽樣方法, 還能幫助深入理解復(fù)雜模型, 并進(jìn)行模型參數(shù)的優(yōu)化、評(píng)估不確定性、進(jìn)行貝葉斯推斷等。能夠有效克服傳統(tǒng)數(shù)值模擬方法的局限性, 幫助解決復(fù)雜問題。
隨著國內(nèi)外的學(xué)者們的探索研究, MCMC 算法在數(shù)值模擬中的應(yīng)用方面出現(xiàn)了更多研究成果, 涉及金融、地震預(yù)測(cè)、貝葉斯統(tǒng)計(jì)、機(jī)器學(xué)習(xí)等領(lǐng)域。在我國, 郭翠將MCMC 算法應(yīng)用于利率期限結(jié)構(gòu)模型, 確定連續(xù)時(shí)間資產(chǎn)定價(jià)模型中的后驗(yàn)分布, 為金融方面提供了重要參考[1] ; 趙晨等人利用MCMC 算法對(duì)地震彈性抗阻反演方法進(jìn)行優(yōu)化,結(jié)合地質(zhì)特征, 對(duì)地下儲(chǔ)存層參數(shù)進(jìn)行預(yù)測(cè), 為油氣勘探提供了參數(shù)支撐[2] 。在國外, Andrew Gelman等人在貝葉斯統(tǒng)計(jì)建模中廣泛應(yīng)用MCMC 算法, 提出了一系列基于MCMC 的貝葉斯建模方法, 為貝葉斯統(tǒng)計(jì)的發(fā)展做出了重要貢獻(xiàn)[3] ; Radford Neal 在機(jī)器學(xué)習(xí)領(lǐng)域利用MCMC 算法進(jìn)行模型訓(xùn)練和推斷, 提出了一種基于MCMC 的變分推斷方法, 能夠高效地處理復(fù)雜的概率模型[4] ; Christian Robert 等人提出了一種基于MCMC 的貝葉斯模型比較方法,能夠比較復(fù)雜模型之間的相對(duì)優(yōu)劣, 為統(tǒng)計(jì)學(xué)中的模型選擇問題提供了一種有效的解決方案[5] 。這些研究為MCMC 算法在數(shù)值模擬中的應(yīng)用提供了豐富的理論基礎(chǔ)和實(shí)踐經(jīng)驗(yàn), 為其在實(shí)際問題中的應(yīng)用奠定了堅(jiān)實(shí)的基礎(chǔ)。
通過深入研究MCMC 算法在數(shù)值模擬中的應(yīng)用, 不僅可以為理論研究提供新的思路和方法, 還可以為很多領(lǐng)域提供重要的參考和指導(dǎo)。因此, 深入探討MCMC 算法的原理和應(yīng)用現(xiàn)狀具有重要的理論和實(shí)踐意義, 對(duì)推動(dòng)相關(guān)領(lǐng)域的發(fā)展和應(yīng)用具有重要的促進(jìn)作用。
二、MCMC 算法及相關(guān)原理
MCMC 算法由蒙特卡羅方法和馬爾可夫鏈兩個(gè)MC 組成。因此, 要介紹MCMC 的原理, 需先引入蒙特卡羅方法和馬爾可夫鏈的相關(guān)原理, 為深入理解MCMC 做理論基礎(chǔ), 后介紹MCMC 算法的原理。
( 一) 蒙特卡洛方法
“曼哈頓計(jì)劃” 時(shí)期, 為了探究中子在原子彈內(nèi)擴(kuò)散和增值問題??茖W(xué)家們構(gòu)建了高維玻爾茲曼方程, 并將解決高維玻爾茲曼方程的問題轉(zhuǎn)化為相等形式的隨機(jī)問題, 并通過計(jì)算機(jī)模擬成功求得解, 這便是最初的蒙特卡洛方法[6] 。這種方法很像賭場(chǎng)里扔骰子的過程, 因此梅特羅波利斯提議用賭城的Monte Carlo 作為代號(hào), 這就誕生了蒙特卡羅方法[7] 。
蒙特卡洛方法的中心思想是通過從合適的分布上進(jìn)行抽樣, 生成樣本來模擬問題的不確定性或者隨機(jī)性。再基于生成的樣本, 使用統(tǒng)計(jì)方法對(duì)問題進(jìn)行推斷, 估計(jì)目標(biāo)函數(shù)的期望值。從密度函數(shù)p(·) 中構(gòu)造樣本{xt , t = 1, 2, , N} 估計(jì)E[f(xi )] :
即用樣本估計(jì) f(x) 的平均值。 當(dāng)樣本 {xi} 之間相互獨(dú)立, 由格里文科定理, 隨樣本量增加, 經(jīng)驗(yàn)分布函數(shù)收斂于真實(shí)的總體分布函數(shù)。樣本均值在樣本量趨于無窮時(shí)依分布收斂于總體期望是大數(shù)定理的核心。兩個(gè)理論相結(jié)合可得, 采集樣本的數(shù)量增加, 估計(jì)值就越趨于真實(shí)值, 對(duì)分布的模擬就更加準(zhǔn)確[8] 。
在深厚的理論基礎(chǔ)支撐下, 面臨著一個(gè)核心挑戰(zhàn): 如何對(duì)目標(biāo)分布進(jìn)行采樣, 以顯著提升擬合效果。為此, 就要引入蒙特卡洛方法。蒙特卡洛采樣技術(shù)通常包括三種主要策略: 直接采樣、接受-拒絕采樣以及重要性采樣。以下給出了三個(gè)抽樣方法的原理:
1. 接采樣方法
假如p(z) 是要采樣的分布, 若可以得到p(z)的概率密度(probability density function, PDF), 對(duì)概率密度函數(shù)求積分可以得到分布函數(shù)的累積分布函數(shù)(cumulative distribution function, CDF), 假設(shè)u(i) ~ U(0, 1) , 即u 是(0, 1) 的均勻分布, 則取得樣本x(i)
x(i) = CDF-1(u(i) ) (2)
對(duì)PDF 求積分后, 得到分布函數(shù)CDF , 則可求出u(i) 對(duì)應(yīng)的樣本x(i) , 同理可以抽樣出N 個(gè)樣本點(diǎn)x(1) x(2) , …, x(N) , 但是, 現(xiàn)實(shí)中很難得到問題的概率密度函數(shù), 也就無法求出累計(jì)分布函數(shù)。因此, 這種方法只能適用于一些非常有限的分布, 很多實(shí)際問題很難通過這種方法去抽樣得到樣本點(diǎn)。
2. 接受拒絕采樣方法
接受拒絕采樣的核心思想, 是對(duì)較復(fù)雜的概率分布p(z) , 引入簡(jiǎn)單的提議分布q(z) 和一個(gè)常數(shù)k , 使得對(duì)于任意的Zi , 有kq(zi ) ? p(zi ) , 然后對(duì)q(z) 進(jìn)行采樣獲得樣本, 具體的采樣方法步驟如下:
第一步, 從q(z) 取出一個(gè)樣本, 記為zi ; 第二步, 按照提議分布q(z) 隨機(jī)抽樣得到樣本zi ,計(jì)算接收率α =p(zi )/kq(zi ) , 再按照均勻分布在(0, 1)范圍內(nèi)抽樣得到樣本ui ; 第三步, 如果ui ≤(p(zi )/kq(zi ) ),則zi 作為抽樣結(jié)果; 否則, 拒絕zi 并返回第二步,直到獲得N 個(gè)樣本, 則采樣結(jié)束。
每次采樣的接受概率計(jì)算如下:
則為了防止舍棄過多樣本降低采樣效率, k 的選取應(yīng)該在滿足ui ≤(p(zi )/kq(zi ) )的基礎(chǔ)上盡可能小。
拒絕采樣的優(yōu)點(diǎn)是容易實(shí)現(xiàn), 缺點(diǎn)是采樣效率可能不高。如果p(zi ) 的涵蓋體積與kq(zi ) 的涵蓋體積比值低, 就會(huì)增大拒絕的比例, 降低抽樣效率。尤其是在高維空間抽樣, 即使p(zi ) 與kq(zi )很接近, 兩者涵蓋體積的差異也可能很大, 出現(xiàn)維數(shù)災(zāi)難, 因此, 科學(xué)家們想出了其他采樣方法。
3. 重要性采樣方法
累積分布函數(shù)不可的情況可由接受拒絕采樣解決, 但對(duì)于無法找到適當(dāng)提議分布或者是在高維中采樣的情況, 也無法使用接受拒絕采樣法。由此,出現(xiàn)了重要性采樣法, 步驟如下:
假設(shè)隨機(jī)變量的概率密度函數(shù)為p(x) , 要估計(jì)目標(biāo)分布下某個(gè)函數(shù)f(x) 的期望, 但是從目標(biāo)分布中直接采樣可能很困難, 即可以選擇一個(gè)容易采樣的提議分布q(x) , 然后利用重要性權(quán)重來修正從提議分布中采樣得到的樣本, 以更好地逼近目標(biāo)分布。
對(duì)于一個(gè)函數(shù)f(x) , 其在目標(biāo)分布下的期望值E[f(x)] 可以表示為:
Ep(x) [f(x)] = ∫f(x)p(x)dx (4)
再引入較為簡(jiǎn)單的提議分布q(x) , 對(duì)上式進(jìn)行變換
記重要性權(quán)重wi =p(xi )/q(xi ) , 其中, xi 是從提議分布 q(x) 中采樣得到的樣本。 估計(jì) E [f(x) ] 的公式為:
重要性采樣解決了概率密度復(fù)雜的情況, 還能經(jīng)常被應(yīng)用于一些需要減小方差的情況, 加權(quán)不僅減小了拒絕接受采樣中提議分布和概率分布相差過大而帶來的誤差, 還能有效減小方差[9] 。
蒙特卡洛方法要求樣本{xi} 相互獨(dú)立, 而由于無法判斷p(·) 是否標(biāo)準(zhǔn)化, 因而從p(·) 中構(gòu)造獨(dú)立樣本不太可行。若要強(qiáng)制要求相互獨(dú)立, 可以從任意分布過程中產(chǎn)生樣本{xt} , 其中一種方法是將含有p(·) 的馬爾可夫鏈看作一個(gè)固定分布, 這就是上文所闡述的馬爾可夫鏈蒙特卡洛算法。
(二) 馬爾可夫鏈
馬爾可夫鏈(Markov Chain) 的發(fā)現(xiàn)可以追溯到19 世紀(jì)末和20 世紀(jì)初, 以俄國數(shù)學(xué)家安德烈·馬爾可夫(Andrey Markov) 的名字命名, 他是這一概念的創(chuàng)始人和最早的研究者之一。
1880 年, “馬爾可夫過程” 的概念首次被馬爾可夫提出。1906 年, 馬爾可夫發(fā)表了一篇名為《關(guān)于離散值連續(xù)變化的正定連續(xù)函數(shù)的問題》(Problems on the Theory of Probability) 的論文, 其中他詳細(xì)討論了一種離散時(shí)間的隨機(jī)過程, 這就是現(xiàn)在所稱的馬爾可夫鏈。隨著后續(xù)數(shù)學(xué)家的不斷探索, 馬爾可夫鏈理論被不斷完善發(fā)展[8] ?,F(xiàn)如今,馬爾可夫鏈在許多領(lǐng)域, 如概率論、統(tǒng)計(jì)學(xué)、物理學(xué)、計(jì)算機(jī)科學(xué)等方面都有廣泛應(yīng)用。為深入理解馬爾可夫鏈, 需要如下前提:
1. 馬爾可夫性質(zhì)
給定當(dāng)前狀態(tài), 未來狀態(tài)的概率分布僅依賴于當(dāng)前狀態(tài), 與其他狀態(tài)無關(guān)。是馬爾可夫鏈的關(guān)鍵特征, 條件概率表示為,
P(Xn+1 = j ∣X0, X1, …, Xn)= P(Xn+1 = j ∣Xn)(7)
其中, Xn 表示第n 個(gè)狀態(tài)。
2. 狀態(tài)空間
x(n) = (x1 (n), x2 (n), … , xk (n)) (8)
其中, 概率向量的每個(gè)元素都是概率, 元素和為1, 系統(tǒng)的可能狀態(tài)數(shù)有k 個(gè), 向量中各個(gè)元素分別表示第n 次觀測(cè)時(shí)第i 個(gè)狀態(tài)的概率, x0 (n) 被稱為初始狀態(tài)。
3. 隨機(jī)過程
用T 表示一個(gè)無限實(shí)數(shù)集, 把依賴于參數(shù)t ∈ T的一族隨機(jī)變量{Xt , t ∈ T} 稱為隨機(jī)過程。
4. 轉(zhuǎn)移概率矩陣
對(duì)于馬爾可夫鏈中任意兩個(gè)狀態(tài)i 和j , pij(i,j = 1, 2, … , k) 表示從狀態(tài)i 到狀態(tài)j 的概率, P矩陣元素非負(fù), 且每一行元素之和都為1。
5. 馬爾可夫鏈
假設(shè)序列狀態(tài)是… Xt -2, Xt -1, Xt , Xt +1… ,那么依據(jù)馬爾可夫鏈性質(zhì), 在Xt +1 時(shí)刻的狀態(tài)的條件概率只依賴于前一刻的狀態(tài)Xt , 即:
P(Xt +1 | … Xt -2, Xt -1, Xt ) = P(Xt +1 | Xt )(10)
求出系統(tǒng)中任兩個(gè)狀態(tài)間的轉(zhuǎn)換概率, 即可確定此馬爾可夫鏈模型。
6. 遍歷性
設(shè)有一個(gè)狀態(tài)空間為E 的馬爾可夫鏈{Xn , n =0, 1, 2, … } , 若對(duì)于任意的i, j ∈ E , 有極限
存在, 則代表馬爾可夫鏈具有遍歷性, 遍歷性代表, 在有限步驟內(nèi), 可以從任何一個(gè)狀態(tài)轉(zhuǎn)移到另一個(gè)狀態(tài), 則認(rèn)為這個(gè)馬爾可夫鏈?zhǔn)潜闅v的。
7. 平穩(wěn)分布
設(shè)齊次馬爾可夫鏈{Xn , n = 0, 1, 2, … } 的狀態(tài)空間為{E = 0, 1, 2,… } , P = (pij ) 代表其轉(zhuǎn)移概率矩陣, 若存在一個(gè)概率分布{πj , j ? 0}滿足
則稱{πj , j ?0} 為該馬爾可夫鏈的平穩(wěn)分布。馬爾可夫鏈的平穩(wěn)分布{πj , j ? 0} 又可以表示為π ={π0, π1, π2,…} , 則有π = πP , 又由于平穩(wěn)分布{πj , j ?0} 是一個(gè)概率分布, 則有Σ∞i = 0πj = 1。
馬爾可夫鏈作為一種隨機(jī)過程, 顯著特征是具有馬爾可夫性質(zhì), 并存在于離散的指數(shù)集和狀態(tài)空間之中。當(dāng)馬爾可夫鏈的指數(shù)集變?yōu)檫B續(xù)時(shí), 它通常被稱作馬爾可夫過程, 有時(shí)也被視作馬爾可夫鏈的一個(gè)特定子集。當(dāng)一個(gè)隨機(jī)過程符合馬爾可夫性時(shí), 將其歸類為馬爾可夫過程[10] 。
此外, 不可約性、常返性、周期性和遍歷性也可能是馬爾可夫鏈具有的性質(zhì)。不可約和正常返是嚴(yán)格平穩(wěn)的馬爾可夫鏈性質(zhì), 它只有一個(gè)平穩(wěn)分布[8] 。遍歷馬爾可夫鏈的極限分布收斂于它的平穩(wěn)分布。馬爾可夫鏈達(dá)到收斂后, 它的狀態(tài)分布保持在某一特定的分布上不隨時(shí)間變化, 即滿足馬爾可夫鏈的穩(wěn)定性。這種穩(wěn)定性表明鏈的長(zhǎng)期行為是可預(yù)測(cè)的, 可以通過鏈的平穩(wěn)分布得到。
(三) MCMC 算法
馬爾可夫鏈蒙特卡羅法(Markov Chain MonteCarlo, MCMC), 是蒙特卡羅方法和馬爾可夫鏈的結(jié)合。在蒙特卡洛模擬的框架下, 從后驗(yàn)分布中抽取樣本時(shí), 若這些樣本滿足獨(dú)立性, 依據(jù)大數(shù)定律, 樣本均值會(huì)逐漸逼近其期望值。若樣本之間不獨(dú)立性, 為了進(jìn)行有效的抽樣, 需要采用馬爾可夫鏈的方法, 即MCMC 方法。
在已知的分布過于復(fù)雜導(dǎo)致無法直接生成樣本的情況下, 可根據(jù)它的條件分布生成樣本。這樣生成的隨機(jī)樣本不相互獨(dú)立, 生成的樣本與前一個(gè)樣本相關(guān), 而由于遍歷性條件下, 這些樣本可以看作獨(dú)立樣本, 這種基于馬爾可夫鏈基本性質(zhì)的方法就是馬爾可夫鏈蒙特卡洛算法[7] 。
MCMC 算法的核心原理是利用馬爾可夫鏈的收斂性質(zhì), 通過定義合適的轉(zhuǎn)移核函數(shù)和接受—拒絕策略, 使得從該鏈中抽取的樣本可以近似地表示目標(biāo)概率分布, 基本步驟如下:
首先, 需要選擇初始狀態(tài); 第二步, 確定狀態(tài)轉(zhuǎn)移; 第三步, 判斷接受或者拒絕新狀態(tài), 通常是通過比較新狀態(tài)的概率與當(dāng)前狀態(tài)的概率來完成,如果新狀態(tài)的概率更高, 則接受該狀態(tài), 否則, 根據(jù)一定的概率規(guī)則, 可能會(huì)拒絕該狀態(tài)。第四步,重復(fù)進(jìn)行狀態(tài)轉(zhuǎn)移和接受或拒絕的過程, 直到滿足停止條件為止。第五步, 在MCMC 算法收斂之后,根據(jù)得到的馬爾可夫鏈, 可以進(jìn)行采樣。采樣過程中, 通常會(huì)丟棄一定數(shù)量的初始狀態(tài), 以保證采樣結(jié)果的有效性。
在第五步中, 當(dāng)模擬的馬氏鏈?zhǔn)諗坑谄椒€(wěn)分布時(shí), 才可以進(jìn)行后續(xù)步驟。判斷MCMC 算法收斂性的常用策略主要包括兩種。首先, 可以通過觀察基于未知參數(shù)模擬得到的馬爾可夫鏈的迭代圖來進(jìn)行初步判斷。若該鏈能迅速擺脫初始值的影響, 并在某個(gè)特定值附近保持隨機(jī)波動(dòng), 無明顯的趨勢(shì)或周期性變化, 則可視為算法收斂。其次, 另一種有效的方法是應(yīng)用Gelman-Rubin 方差比檢驗(yàn)。此方法涉及構(gòu)建多條基于不同初始值的平行鏈, 并計(jì)算它們的后驗(yàn)方差和鏈內(nèi)方差。若兩者值接近或相等,那么可以認(rèn)為MCMC 算法具有良好的收斂性[11] 。
MCMC 方法將蒙特卡羅方法與馬爾可夫鏈相結(jié)合, 有效的解決了復(fù)雜分布難以獲取樣本, 在高維中維數(shù)災(zāi)難的問題。下面介紹兩種MCMC 算法的經(jīng)典算法, Metropolis-Hasting 算法(MH 算法) 以及Gibbs Sampling (吉布斯采樣)。
1. Metropolis Hastings 算法
為了解Metropolis Hastings 算法, 需要先詳細(xì)闡述Metropolis 算法。對(duì)于一個(gè)十分復(fù)雜的過程, 其中變量分布為π , 并潛在地對(duì)應(yīng)著一個(gè)復(fù)雜的轉(zhuǎn)移概率矩陣P , 希望通過一個(gè)對(duì)稱的轉(zhuǎn)移矩陣Q , 從簡(jiǎn)單分布θ 中采樣, 以一定概率α 接受采得樣本作為對(duì)(π, P) 過程的細(xì)致平穩(wěn)分布近似, 用于統(tǒng)計(jì)分析。首先, 給定一個(gè)對(duì)稱的提議轉(zhuǎn)移概率矩陣Q , 但該矩陣一般不滿足細(xì)致平穩(wěn)條件, 即
πi × Pi, j ≠ πj × Pj, i (13)
為了構(gòu)造細(xì)致平穩(wěn)條件, 對(duì)上式兩端同時(shí)乘以接受概率α , 得到:
πi × Qi, j × αi, j = πj × Qj, i × αj, i (14)
使得接受概率α 與提議矩陣Q 的組合滿足如下概率轉(zhuǎn)移矩陣:
P ^i, j = Qi, j × αi, j
P ^j, i = Qj, i × αj, i(15)
這樣一來, 所構(gòu)造的轉(zhuǎn)移概率P ^ 滿足細(xì)致平穩(wěn)條件:
πi ×P ^i, j = πj ×P ^j, i (16)
在Metropolis 算法中, 假設(shè)在第t - 1 步得到狀態(tài)樣本x(t -1) = si , 接下來便從θ(s | si ) 中新采集一個(gè)樣本s ^ 。
以概率α = min{1, π(s ^ ) / π(si )} 接受x(t) =s ^ ,否則x(t) = x(t -1) = si 。該過程對(duì)應(yīng)的轉(zhuǎn)移概率滿足:
P ^i, j = Qi, j × αi, j (17)
那么
πi ×P ^i, j = πi × Qi , j × αi , j
= πi × min{1, πj / πi } × Qi , j
= min{πi × Qi , j , πj × Qj, i }
= πj × min{1, πi / πj } × Qj, i
= πj ×P ^j, i (18)
因此, 如上方式構(gòu)造所得狀態(tài)轉(zhuǎn)移矩陣P ^ 滿足細(xì)致平穩(wěn)條件, 這樣獲得的采樣分布就是π 。
而Metropolis 算法中的接受概率α 的值可能會(huì)很小, 延長(zhǎng)算法的迭代次數(shù), 影響效率。為了解決這一問題, Metropolis Hastings 采樣算法誕生。具體來說, MH 算法將式(20) 等式兩邊的接受概率放大。將其中的一個(gè)接受概率設(shè)置為1, 這樣就能保證每次迭代過程中接收新狀態(tài)的概率越大, 加速算法收斂。
αi, j = min{1, πj Qi, j/πi Qj, i } (19)
從上述的抽樣方法可以得到, Metropolis -Hasting 抽樣基于現(xiàn)有的樣本, 即抽樣需要基于的條件概率。這樣取得的樣本不獨(dú)立, 這就是MCMC 方法與蒙特卡洛方法抽樣的一個(gè)很大不同[12] 。
2. Gibbs 采樣
Metropolis-Hasting 算法可以很好的解決蒙特卡羅方法需要的任意概率分布的樣本集的問題, 但是仍有兩個(gè)缺點(diǎn): 一個(gè)是需要計(jì)算接受率, 計(jì)算量大, 效率低。第二個(gè)是對(duì)于高維數(shù)據(jù), 取得特征的聯(lián)合分布比較困難, 吉布斯抽樣(Gibbs Sampling)的誕生就是M-H 采樣的改進(jìn)。
假設(shè)已知隨機(jī)變量X 服從具有參數(shù)Θ = (θ1,θ2) 的分布X ~ p(X | Θ) , Θ 的先驗(yàn)分布為p(θ1,θ2) 。設(shè){x1, x2, … , xn } 是從該分布抽取出的一組樣本, 可以通過貝葉斯定理計(jì)算得到后驗(yàn)分布:
此后驗(yàn)分布就是需要進(jìn)行抽樣的目標(biāo)分布。有兩種抽樣方法, 一種是在多維角度直接對(duì)Θ 進(jìn)行抽樣, 一種是對(duì)Θ 的每一個(gè)維度θi 逐一地進(jìn)行一維的抽樣。第一種是Metropolis-Hasting 算法的多維推廣, 第二種方案就是吉布斯抽樣。
吉布斯采樣通過全概率分布的目標(biāo)分布對(duì)每一個(gè)維度的θi 進(jìn)行抽樣。θ1, θ2 的滿概率分布形式如下:
p(θ1 | θ2. x1, … , xn ) =p(θ1, θ2 | x1, … , xn )/pθ2(θ2)(21)
其中, p(θ1, θ2 | x1, … , xn ) 為Θ 的聯(lián)合后驗(yàn)概率分布, pθ2(θ2) 為后驗(yàn)分布下關(guān)于θ2 的邊緣分布。由于在這個(gè)條件分布中假設(shè)θ2 = a 已知, 因此p(θ2) 實(shí)際為一個(gè)常數(shù), 可以把上式轉(zhuǎn)變?yōu)椋?/p>
p(θ1 | θ2 = a, x1, … , xn )∝ p(x1, …, xn | θ1, θ2)p(θ1, θ2)∝ p(x1, … , xn | θ1, θ2 = a)p(θ1, θ2 = a)(22)
在θ2 = a 已知的前提下, 使用M-H 算法, 甚至更簡(jiǎn)單的算法, 就可以完成對(duì)θ1, θ2 的抽樣。至此, 還需要解決對(duì)不已知θ2 的抽樣。首先, 得到第一個(gè)樣本Θ1 = (θ11, θ12) , 再利用θi2, 從p(θ1 | θ2 =θi2, x1, … , xn ) 中抽取θi +11 。然后, 利用θi +11 , 從p(θ2 | θ1 = θi +11 , x1, … , xn ) 中抽取θi +12 。構(gòu)成Θi +1 , 最后, 重復(fù)上述步驟, 直到抽取的樣本數(shù)量滿足需求。
可以看到, 抽樣的過程中, 每一個(gè)樣本Θi 總會(huì)受到前一個(gè)樣本的影響, 因此通過吉布斯抽樣得到的序列{Θ1, … , Θn } 是一條馬爾可夫鏈。在此基礎(chǔ)上, 可以類似M-H 算法, 利用馬爾卡夫鏈相關(guān)的性質(zhì)能夠證明, 該馬爾可夫鏈的平穩(wěn)分布是目標(biāo)聯(lián)合概率分布p(θ1, θ2) 。
Gibbs 采樣是基于M-H 采樣的發(fā)展, 在高維更具有優(yōu)勢(shì), 目前應(yīng)用更為廣泛。此外, Gibbs 采樣要求數(shù)據(jù)至少是二維的, 一維概率分布采樣仍然需要使用M-H 采樣。Gibbs 采樣可以獲取概率分布的樣本集, 蒙特卡羅方法可以用樣本集模擬, 二者結(jié)合, 突出了MCMC 算法在大數(shù)據(jù)時(shí)代, 尤其在處理高維數(shù)據(jù)模擬求和任務(wù)時(shí), 獨(dú)特的優(yōu)勢(shì)和作用。
三、MCMC 算法的應(yīng)用
(一) 蒙特卡洛方法的應(yīng)用
蒙特卡洛方法被廣泛應(yīng)用于估計(jì)、優(yōu)化、積分等問題, 在統(tǒng)計(jì)學(xué)、物理學(xué)、金融學(xué)等領(lǐng)域廣泛應(yīng)用, 尤其在高維、復(fù)雜問題上表現(xiàn)出色。以下介紹幾種蒙特卡羅方法在數(shù)值模擬中的經(jīng)典應(yīng)用。
1. 蒙特卡洛方法擬合正態(tài)分布的參數(shù)
利用Python 軟件分別生成包含100 個(gè), 1000,10000 個(gè)樣本的正態(tài)分布數(shù)據(jù), 均值為10, 標(biāo)準(zhǔn)差為2。后定義了一個(gè)函數(shù), 用于執(zhí)行蒙特卡洛擬合。這個(gè)函數(shù)接受樣本數(shù)據(jù)和模擬次數(shù)作為參數(shù), 然后在樣本數(shù)據(jù)中進(jìn)行隨機(jī)的有放回地抽樣, 計(jì)算均值和標(biāo)準(zhǔn)差, 并返回估計(jì)的均值和標(biāo)準(zhǔn)差, 結(jié)果如表1:
可以看出, 隨著樣本數(shù)量的增加, 蒙特卡洛模擬對(duì)參數(shù)的擬合更準(zhǔn)確。在實(shí)際應(yīng)用中, 可以通過增加抽樣次數(shù)和改變初始參數(shù)值來提高擬合的準(zhǔn)確性。
2. 蒙特卡洛方法求解不規(guī)則圖形的面積
如圖2, 若想對(duì)紅色的面積進(jìn)行計(jì)算, 需要在x∈[0, 2], y ∈[0, 8] 的區(qū)域均勻隨機(jī)撒點(diǎn), 如果點(diǎn)進(jìn)如紅色區(qū)域, 則計(jì)數(shù)為1, 否則為0, 最后只要將撒在紅色區(qū)域的數(shù)量除以撒點(diǎn)的總數(shù)N 除以滿足撒點(diǎn)總數(shù)再乘上x ∈ [0, 2], y ∈ [0, 8] 的面積,則可得到紅色區(qū)域的面積。
以下分別進(jìn)行100 次, 1000 次, 10000 次,100000 次的隨機(jī)撒點(diǎn), 得到模擬的結(jié)果以及根據(jù)撒點(diǎn)情況計(jì)算的紅色面積如圖3;
圖3 顯示: 隨著撒點(diǎn)數(shù)量的增加, 區(qū)域內(nèi)被擊中部分的面積也不斷增加, 可以得到更加準(zhǔn)確的數(shù)字, 為了更精確的求得紅色區(qū)域面積, 將N 次撒點(diǎn)計(jì)數(shù)為1 的累計(jì)和除以撒點(diǎn)的總數(shù)N 再乘上x ∈[0, 2], y ∈ [0, 8] 的面積, 則可得到紅色區(qū)域的面積, 結(jié)果如表2:
需要注意, 隨機(jī)擬合得出的面積是一個(gè)近似值, 且隨著樣本數(shù)量的增加, 近似值更接近真實(shí)值, 上述應(yīng)用的采樣要求均勻隨機(jī), 但并不是所有MC 采樣都需要應(yīng)用均勻采樣, 如重要性采樣,因此, 通過應(yīng)用特定的非均勻采樣策略, 可以顯著提高算法的收斂速度。
3. 接受拒絕采樣的應(yīng)用
定義函數(shù)f1(x) = 0. 3exp( - (x - 0. 3)2) +0. 7exp( - (x - 2)2/0. 3 )/k , 生成從-4 到6 的橫坐標(biāo)范圍, 步長(zhǎng)為0. 01, 將橫坐標(biāo)范圍x 和對(duì)應(yīng)的概率密度函數(shù)f1(x) 繪制成圖像設(shè)置為紅色。
設(shè)置抽樣數(shù)量10000000, 設(shè)置正態(tài)分布的標(biāo)準(zhǔn)差為1. 2, 均值為1. 4, 使用函數(shù)生成符合指定均值和標(biāo)準(zhǔn)差的正態(tài)分布隨機(jī)樣本z。根據(jù)生成的正態(tài)分布隨機(jī)樣本z, 使用概率密度函數(shù)計(jì)算對(duì)應(yīng)的概率密度函數(shù)值qz 。設(shè)置k = 2. 5, 根據(jù)k 和qz 生成均勻分布的隨機(jī)樣本u , 根據(jù)生成的正態(tài)分布隨機(jī)樣本z , 使用定義的概率密度函數(shù)計(jì)算對(duì)應(yīng)的真實(shí)概率密度函數(shù)值pz , 根據(jù)生成的隨機(jī)樣本u 和真實(shí)概率密度函數(shù)值pz , 選取滿足條件pz ? u 的樣本, 得到抽樣結(jié)果sample, 繪制直方圖, 對(duì)直方圖進(jìn)行歸一化處理, 使得面積為1, 結(jié)果如圖4。
則采樣擬合的分布與原分布相比基本一致, 且擬合效果隨著樣本數(shù)量增加越來越好, 以上是蒙特卡洛方法中接受拒絕采樣在抽樣中的應(yīng)用。
(二) 馬爾可夫鏈的應(yīng)用
馬爾可夫鏈在實(shí)踐中有許多應(yīng)用, 如: 馬爾可夫鏈可用于文本生成, 通過學(xué)習(xí)文本中單詞或字符之間的轉(zhuǎn)移概率來生成新的文本, 還可以用來建模語音信號(hào)中的音素序列, 達(dá)到語音識(shí)別。在金融領(lǐng)域, 馬爾可夫鏈可以用來模擬資產(chǎn)價(jià)格的隨機(jī)演變過程, 被用于股票價(jià)格預(yù)測(cè)、風(fēng)險(xiǎn)管理和期權(quán)定價(jià)等。在生物信息學(xué)中, 馬爾可夫鏈可以用來建模DNA、RNA 和蛋白質(zhì)序列的序列變化和序列結(jié)構(gòu)[14] 。在網(wǎng)絡(luò)分析中, 馬爾可夫鏈可以用來建模網(wǎng)絡(luò)中節(jié)點(diǎn)之間的轉(zhuǎn)移過程。此外, 馬爾可夫鏈還可以用于預(yù)測(cè)天氣和傳染病病情。下面以一模型為例, 有A, B, C 有三家租車的門店, 租賃和歸還可以選擇任何一個(gè)門店, 從不同門店借出和歸還的概率如表3:
由表3 的數(shù)據(jù)可以寫出轉(zhuǎn)移矩陣為
則可確定還車模型, 接下來用此模型進(jìn)行預(yù)測(cè), 已知一輛車位于B 門店, 三次借還以后, 計(jì)算這輛車最可能在哪個(gè)門店。
若初始狀態(tài)向量為:
X(0) = (0, 1, 0) (24)
則
X(1) = X(0) P = (0. 300, 0. 100, 0. 600)
X(2) = X(1) P = (0. 360, 0. 430, 0. 210)
X(3) = X(2) P = (0. 372, 0. 241, 0. 387)
(25)
故三次借還后這輛車最可能在C 門店, 用Python 軟件模擬, 借還1 次至20 次后, 此車停在三個(gè)門店的概率結(jié)果如表4 所示:
表4 顯示: 當(dāng)n 增大, 狀態(tài)向量趨近于
X(n) = (0. 375, 0. 300, 0. 325) (26)
則多次歸還后, 從A, B, C 三家門店歸還的概率分別趨近于0."375、0. 3、0. 325, 達(dá)到平穩(wěn)狀態(tài)。
(三) MCMC 算法的應(yīng)用
使用MCMC 算法進(jìn)行數(shù)值模擬被廣泛應(yīng)用。在統(tǒng)計(jì)學(xué)習(xí)領(lǐng)域, MCMC 算法在廣泛應(yīng)用與貝葉斯統(tǒng)計(jì)推, 用于從后驗(yàn)分布中抽樣。在數(shù)學(xué)領(lǐng)域,MCMC 算法可以用于計(jì)算高維積分, 特別是對(duì)于復(fù)雜的概率分布。在計(jì)算機(jī)領(lǐng)域, MCMC 算法用于對(duì)馬爾可夫隨機(jī)場(chǎng)(MRF-Markov Random Field , 簡(jiǎn)稱MRF) 進(jìn)行建模和推斷[15] 。在生物信息學(xué)領(lǐng)域: MCMC 算法可以用來模擬蛋白質(zhì)的折疊過程,從而找到最穩(wěn)定的構(gòu)象。在金融學(xué)領(lǐng)域, MCMC 算法對(duì)金融模型進(jìn)行參數(shù)估計(jì)、風(fēng)險(xiǎn)評(píng)估等。在機(jī)器學(xué)習(xí)領(lǐng)域, MCMC 算法可以用于擬合深度學(xué)習(xí)模型的參數(shù)、進(jìn)行變分推斷等。
下面以線性回歸模型為例, 用MCMC 算法來估計(jì)回歸系數(shù)。先定義一個(gè)線性模型linear_ regression(x, beta), 其中x 是自變量, beta 是回歸系數(shù), 選擇了一組真實(shí)的回歸系數(shù)true_ beta = [2, 1], 為了生成觀測(cè)數(shù)據(jù), 在自變量x 的范圍內(nèi)生成了一組數(shù)據(jù)點(diǎn), 并通過線性回歸模型加上正態(tài)分布的噪聲來生成因變量y_ observed。具體來說, 在0 到10之間生成了50 個(gè)等間隔的自變量數(shù)據(jù)點(diǎn)。然后,使用linear_ regression 函數(shù)和真實(shí)的回歸系數(shù)true_ beta 來計(jì)算相應(yīng)的因變量y_ true。接著, 生成了服從正態(tài)分布的噪聲數(shù)據(jù)noise, 然后將其加到真實(shí)的因變量y_ true 上, 得到最終的觀測(cè)數(shù)據(jù)y_ ob?served。將合成后的數(shù)據(jù)進(jìn)行可視化, 并對(duì)散點(diǎn)進(jìn)行擬合, 結(jié)果如圖5 所示。
接下來, 使用MCMC 算法估計(jì)線性回歸模型中的參數(shù), 得到參數(shù)的后驗(yàn)分布樣本之后, 可計(jì)算參數(shù)的后驗(yàn)均值、標(biāo)準(zhǔn)差等統(tǒng)計(jì)量, 也可以使用這些樣本來進(jìn)行后續(xù)的推斷和預(yù)測(cè)。利用Python 軟件進(jìn)行數(shù)值模擬, 結(jié)果如圖6 所示。
可得到截距和斜率的后驗(yàn)分布。這就是MCMC算法在數(shù)字領(lǐng)域的應(yīng)用實(shí)例, 顯示了MCMC 算法在解決各種復(fù)雜問題中的強(qiáng)大能力, 學(xué)習(xí)MCMC 算法在各個(gè)領(lǐng)域都有著重要的意義。
四、總結(jié)
使用MCMC 算法進(jìn)行數(shù)值模擬中在各個(gè)領(lǐng)域都有強(qiáng)大的應(yīng)用。通過構(gòu)建馬爾可夫鏈, MCMC 算法能夠從復(fù)雜的概率分布中抽樣, 解決了許多傳統(tǒng)方法難以處理的問題。在過去幾十年里, MCMC 算法已經(jīng)被廣泛應(yīng)用于各個(gè)領(lǐng)域, 包括統(tǒng)計(jì)學(xué)、機(jī)器學(xué)習(xí)、貝葉斯推斷、物理學(xué)、生物學(xué)等。其應(yīng)用涵蓋了參數(shù)估計(jì)、模型選擇、數(shù)據(jù)合成、圖像處理等多個(gè)方面, 為科學(xué)研究和工程應(yīng)用提供了重要支持。
雖然MCMC 算法已經(jīng)取得了巨大成功, 并在許多領(lǐng)域得到了廣泛應(yīng)用, 但仍然存在一些挑戰(zhàn)和改進(jìn)的空間。首先, 隨著數(shù)據(jù)的增加, MCMC 算法和并行化方法需要被開發(fā), 應(yīng)對(duì)大規(guī)模數(shù)據(jù)帶來的挑戰(zhàn)。因此, 需要開發(fā)適用于大規(guī)模數(shù)據(jù)的, 以加速計(jì)算和提高效率。其次, 改進(jìn)MCMC 算法的采樣效率是當(dāng)前研究的一個(gè)重要方向。通過引入新的采樣策略、優(yōu)化算法參數(shù)或并行計(jì)算等手段, 提高算法的收斂速度和采樣效率, 尤其是在高維參數(shù)空間中的應(yīng)用。再者, MCMC 算法在貝葉斯推斷中扮演著重要角色, 但在復(fù)雜模型和大規(guī)模數(shù)據(jù)情況下, 后驗(yàn)推斷可能變得困難。因此, 需要進(jìn)一步研究高效的后驗(yàn)推斷方法, 包括變分推斷、近似推斷等, 以應(yīng)對(duì)復(fù)雜模型和大規(guī)模數(shù)據(jù)的挑戰(zhàn)。最后, MCMC算法已經(jīng)在多個(gè)學(xué)科領(lǐng)域得到了廣泛應(yīng)用, 但不同領(lǐng)域之間的交叉應(yīng)用仍然有待深入探索。未來的研究可以關(guān)注跨學(xué)科合作, 將MCMC 算法與其他領(lǐng)域的方法結(jié)合, 為更廣泛的應(yīng)用場(chǎng)景提供解決方案。
參考文獻(xiàn):
[1]郭翠 MCMC 方法在利率期限結(jié)構(gòu)模型中的應(yīng)用[D]. 暨南大學(xué),2012.
[2]趙晨,郭淑文,金鳳鳴等. 基于優(yōu)化MCMC 算法的地震彈性阻抗反演方法[J]."礦產(chǎn)與地質(zhì),2023,37(03):587-596.
[3]Gelman A , Carlin J B , Stern H S ,et al. Bayesian data analysis, third edition[ J] . Journal of the American Statistical Association, 2003, 45(2):21-100.
[4]Neal R M . Bayesian leaning for neural networks[J]. Mackay Neura/ Computation, 1996:46-152.
[5]Reviewer D P F . Monte Carlo Methods: Monte Carlo Methods in Statistical Physics[J]. Computing in Science & Engineering, 2000, 2(6):73-74.
[6]劉朝暉. 計(jì)算機(jī)誕生,推動(dòng)人類文明進(jìn)階[J]. 新民周刊,2023,(05):36-37.
[7]周將. 基于蒙特卡洛模擬的防洪設(shè)施可靠性分析[J]."價(jià)值工程,2024,43(10):150-153.
[8]周艷斌 基于馬爾可夫鏈蒙特卡洛法的螞蟻游走模型[D]. 電子科技大學(xué),2023
[9]邵豪,王倫文,朱然剛等 基于重要性采樣的超圖網(wǎng)絡(luò)高效表示方法[J]. 軟件學(xué)報(bào),1-18.
[10]肖北芳. 基于MCMC 算法的貝葉斯分位數(shù)回歸方法的實(shí)證應(yīng)用[D]. 湖南師范大學(xué),2015.
[11]楊新輝. 隨機(jī)抽樣的方差縮減以及MCMC 收斂診斷[D]. 北京交通大學(xué),2018.
[12]張廣智,趙晨,涂奇催,等. 基于量子退火Metropolis-Hastings 算法的疊前隨機(jī)反演[J]. 石油地球物理勘探,2018,53(01):2-9.
[13]趙啟鈐. 基于吉布斯采樣算法的虛擬樣本生成技術(shù)研究[D]. 北京化工大學(xué),2023.
[14]劉洪艷,劉良森. 生物化學(xué)教學(xué)案例的構(gòu)建—人工智能預(yù)測(cè)蛋白質(zhì)結(jié)構(gòu)[J]. 化學(xué)教育(中英文),2024,45(02):92-97.
[15]李兮嘉. 基于馬爾科夫鏈蒙特卡洛算法對(duì)條件生成對(duì)抗網(wǎng)絡(luò)生成圖片質(zhì)量的改進(jìn)[D]. 山東大學(xué),2023.