实例介绍
【实例简介】
对马尔科夫蒙特卡洛模拟算法的描述，比较好的参考资料
法(一般通稱MCMC演算法),配合練習題,將使讀者對MCMC有淸楚與概括的暸 解。 Metropolis method: Metropolis演算法的穊念很簡單,可以下列步驟來說明 合X0爲任何變數的樣本 2.第k+1個變數的樣本局 y k=0.1 其中s服從均等分配,即s~U(-a,a 3.產生一{個亂數~U(0,1) 4.如果 a< r(y) (X) Xk+1=y,否則Xk+1=Xk(即下一個變數仍稚持Xk),其中丌()焉目標 機率函數。 Metropolis溟算法有幾個地方值得注意 步驟2的均等分配鮑圍(-a,a)的選擇決定了收斂的速度、若a過小,馬可夫 鍊的移動太慢,收斂比較慢。過大的α容易導致步驟4的拒絕率升高,Xk+1傾 向不變動。 ·透過上的演算法將產生一系列的樣本,最初的樣本與目標機窣匦數闌倸不大, 在不斷的透過步驟4的限制,樣本將慢慢具備服從目標禨羍函數的特性。所以當 我們希望達穠定狀態後的馬可夫鍊樣本,Xm+1,Xm+2,…,服從目標機夲函數π()分配。 以樣本邳均估計期望值時,通常會猞棄前m個樣本(俗稱burn- cin samples), 即 n- 7r ∑X k=m+1 這當然是因爲前面的樣本還在調整修正中,未達穩定狀態,不適合當作目標機率 函數的樣本。 Metropolis-Hastings method Hastings修正了 Metropolis對樣本的產生(步驟2)與認定(也就是步驟4的「按 受或拒絕」的準則)。 Metropolis-Hastings演算法如下 1.合X0鳥任何變數樣本 2.第k+1個變數的樣本爲 (YXk) 3.產生一數~U(0,1) 4.如果 us Tk 則Xk+1=y,否則X+1=Xk(即下一個變數仍雜持Xk),其中 (Y)qXkY 7k=m(·x(x)(YXA) g(·-)爲提議機羍函數( proposal distribution)或轉換函數( transition func tion). Markov chain的轉換機羍寫成?(-)=q(|-)r 新樣本來自提議機羍函數q(Y|X),這個函數的選擇當然扮演重大的角色,神奇的是, q(YX)可以是任何型態的函數,其所產生的 Markov chain達穩定狀態後(sta tionary)的分配將是丌()。前面 Metropolis演算法所探用的方式一般稱鳥 random walk chain.也是q(YX)的選擇之一。不過適當的選擇卻關係到收斂( (stationary) 的速度,一般而言,選擇q(YX)的原則不外乎 q(YX)愈按近目標機函數丌(Y)愈好(如範例1所示)。 容易自q(YX)產生樣本 Metropolis演算法是 Metropolis- Hastings演算法的特例,主要的差别在其具對稱 性的提議機羍函數,即q(YX)=q(X1Y),這個假設影響了樣本的接受/拒絕率。2以 上沈 Metropolis演算法步驟2鴦例,新樣本的產生可以表小鳥 y Nq(Y XR=q(Y-XkD=s s~U(0,a) 另外,q(Y|X)=N(X,a2)的常態假設也是常見的特例,因盒q(Y1X)=q(XY) 以下利用簡單的實作問題來瞭解MCMC的運作方式。 範例1:假設目熛函數π()鹪標凖常態N(0,1),寫一支程式利用 Metropolis演算 法產生 Markov chain樣本並測武當提議分配q(YX)分别焉 1.如 Metropolis演算法步驟2的 Random walk q(Y|X)=q(Y-X1) 2.q(YX)=N(X,0.5),N(X,0.1),N(X,10)三種 畫圖觀察 Markov(hain達穩定狀態的情況。 圖1展示了三種常態分配作爲提議分配的表現,不難看岀當提議分配愈接近目標分配 時,其收斂的速度比較理想。睛注意這些常態分配都是以 Markov chain變數ⅹ做 鳥中心點,並非固定。即新旳樣本從「以前一個樣本做篤均數」的常態分配中抽樣。 從圖1也可以觀察到經過足夠艮(多)的時問(樣本)以後,譬如100個以後. Markov hain產生的相依性樣本愈接近從目標函數所產生的。這時便可以利用樣本不均數來 計算母體(日標囪數)均數的期望值F(X),或是進一步的流計量估計F(f(X)。而 前一段被捨棄的100個標本稱 Burn-in samples MATLAB也提供指合「 sample產生樣本,試著使用看看,觀察與自己寫的程式 所產生的結果有何不同。譬如 2 r(Yq(Xky 丌(xk)q(Y|Xk)-m(Xk) 100270040500001 圖1: Metropolis演算法:不同提議分配的收斂表現 sigma=sqrt(.5) tarpdf =@(x)normpdf(x); target pdf proppdf=@(x, y) normpdf(y, x, sigma ) %proposal pdf qlyl r proprnd=@(x)normrnd (x, sigma ) generate samples from proposal pdf nsamples= 1500 z= mhsample(1, nsamples,pdf;. tarpdf, proprnd', proprnd, 'symmetric, 1) histfit(z, 50) 上逃「 sample」指命選擇了對稱性的提議分配,所以第三行對提議分配( propp) 的設定實際上是用不的,其秸果與一般的用法相同。 z= mhsamplc(1, nsamplcs, ' pdf, tarpdf, proprnd, proprnd, proppdf', proppdf ); sample」指爷除輸出樣本外,也可以同時轍出「 Acceptance rate,以便觀察提 議分配的表現,讀者不妨參考線上手册關於 sample的其他選項。 1. 4 Data Augmentation 當後驗機窣密度函數(1)沒有分析解,需要訴諸數值解時, Data. Augmentation的 概念十分合適。如同FM演算法所引入的隱藏變數造成較箐單的概似函數型態,Data Δ augmentation也是試圖在後驗機率的計算上引入隱藏蘷數,讓這個計算變得可能與 簡罩,其觀念如下式 P(ey) 0z, yp(zly ) 從p(θy)到p(θz,y)引進了隱藏變數z使得p(θz.y)變得更谷易分析、計算或是 抽樣,當然最好是一個已知型態的分配。引進了z,卻也同時多出個積分式,當式(4 的積分不容易做到時 Monte carlo積分是不錯的選擇,於是式(4)可以看成 p(By P(az, yp(zy)dx ((6z,y)) Ave(p(ezi, y)) 其中藏變數的資料z;是從所謂的 preditive distribution p(zy)抽樣而來。p(zy) 可以寫成 (zy)=/(z0,y)(y)d 上式意謂著z;的產生通常從p(zθ,y)更為方便,當利用 Monte carlo積分時,積 分變數θ的樣本必須從γ(θy)抽樣而來。這於是形成了「雞生蛋,蛋生雞」的循環遞 迴問題:p(硎y)的計算依靠隱藏變數樣本z,而z;從p(θy)抽樣而交。這個Data augmentation演算法可以寫成下列的兩個步驟 I step( Imputation Step):執行下列兩個步驟m次,產生z1,Z2,…,Zm (1)從p(y)抽出一個樣本θ (2)依步驟(1)產生的6;但,從p(z,y)抽出一個z P step( Posterior Step):更新後驗機率鳥 p(y)= ∑(6z,y) 在 I step中有一個m需要決定, Tanner&wong的論文[2中說明m的選擇可以 是變動的,甚至低到m=1都可以保證迥圈移動的方向是正確的。另外在 P step中 更新的後驗機羍可以被看成是一個混合機羍密度囡數( mixture distribution),而在 L-step中的p(z,y)必須是容易計算或容易抽樣的密度函數,這個演算法扌有實質 的意義。這個演算法被證明凹2在大部分的假設情況下,都能收斂到理論的後驗機羍密 度函數。 範例2:舉著名的 Genetic Linkage Model鳥例,197隻動物以多項分配的分佈屬於 四種類別,分別是 y=(y,y,y3,y4)=(125,18,20,34) 其各項的楼李分別是 6× 61b6 欲估計未知參數θ 最瘟接的作法是訴諸最大概似函數的佔計, L(y 0)= max L(y 6) 0<6<1 其中概似數 Ly0)∝(2+6)12(1-0)0p 參數θ的最大概似估計可以由解概似函數的一次導數爲0,或是輾轉的以EM涣算 法計算可得。另一種作法是更直接的以其件式期望值作估計,即 E(0ly)=/p(ely)de 其中p(0ly)焉後驗機羍。上式可以迤一步寫成 ∫ep(y|O)xr(6)1 6=E(O|y)= π(θ)爲先驗機羍。上式又稱貝氐估計。貝氏估計經常面臨到積分的問題,特別當變 數維度很高時,積分的難度便高岀許多,往往得不到分析解( analytical solution)。以 本題鶯例,假設對先驗機羍一無所知,即π(0)-L.,0≤0≤I,式(10寫成 he2+0)1201-0)3sed (2+0)25(1-6)3834a0 僅是單一變數,這樣的積分也夠嚇人的了 這類的問題譏MCMC找到發揮的舞台, Monte carlo積分透過抽取大量樣本並計算 其吓均來近似期望值的作法,可以避開繁瑱的積分。但是當抽取樣本所憑藉的後驗機 牽不是典型的機羍密度巫數時,以本範例鳥例 p(|y)∝(2+b)12(1-0)034 樣本的產生便是個問題。還好馬可夫鍊( Markov chain)的概念提供一個解套的方 式,以下先利用1.3介紹的兩種MCMC演算法,看看效果如何? Metropolis algorithm 圄2根據式(12)的機羍密度函數以 Metropolis演算法所產生1000個樣本畫出的 直方圖。直方圖的分佈與式(I2)非常的接近。從樣本計算式(9)的積分式,得到的估 計值如圖3中的紅線所示。當懞本數愈大時,其估計值愈穩定,表示從均等分配產生的 Histogram of the randor semple 10 圄2:利用 Metropolis演算法產生的樣本畫出的直方圖 樣本趨近日標機牽酾數(12)。請注意圖3的紅線是所謂的「累加夲均偵,可以使用 MATLAB的指合「 cumsum」進行累加後冉做4均,譬如 y=cmsm(x)./(1:m)%向量x代表1xm的樣本 Metropolis-Hastings Algorithm 當應用 Metropolis- Hastin⑧s演算法時轉換機羍函數的選擇很關鍵。以本題的目標 機羍函數鳥例(如圖4實線),可以選擇常態分配或選擰Beta(56,38)(如圖4虛線)當 作轉換機率囪數,其實Beta(56,34)作鴦q(Y|X)最駑完美,幾乎與目標機率鹵數 致.以此轉換函數產生的樣本與均數展示在圖5。程式如下 125182034 A=56B=38: dr=((x)(2+x)^y(1)*(1-x)^(y(2)+y(3)*^y(4); ppdf =@(a, y) betapdf(r, A, B) nd=o() betarnd(A, B) samples= mhsample(0.5, 1500,, tpdf, proprnd', prnd,'proppdf', ppdf) histfit(samples, 50) 【实例截图】
【核心代码】
对马尔科夫蒙特卡洛模拟算法的描述，比较好的参考资料
法(一般通稱MCMC演算法),配合練習題,將使讀者對MCMC有淸楚與概括的暸 解。 Metropolis method: Metropolis演算法的穊念很簡單,可以下列步驟來說明 合X0爲任何變數的樣本 2.第k+1個變數的樣本局 y k=0.1 其中s服從均等分配,即s~U(-a,a 3.產生一{個亂數~U(0,1) 4.如果 a< r(y) (X) Xk+1=y,否則Xk+1=Xk(即下一個變數仍稚持Xk),其中丌()焉目標 機率函數。 Metropolis溟算法有幾個地方值得注意 步驟2的均等分配鮑圍(-a,a)的選擇決定了收斂的速度、若a過小,馬可夫 鍊的移動太慢,收斂比較慢。過大的α容易導致步驟4的拒絕率升高,Xk+1傾 向不變動。 ·透過上的演算法將產生一系列的樣本,最初的樣本與目標機窣匦數闌倸不大, 在不斷的透過步驟4的限制,樣本將慢慢具備服從目標禨羍函數的特性。所以當 我們希望達穠定狀態後的馬可夫鍊樣本,Xm+1,Xm+2,…,服從目標機夲函數π()分配。 以樣本邳均估計期望值時,通常會猞棄前m個樣本(俗稱burn- cin samples), 即 n- 7r ∑X k=m+1 這當然是因爲前面的樣本還在調整修正中,未達穩定狀態,不適合當作目標機率 函數的樣本。 Metropolis-Hastings method Hastings修正了 Metropolis對樣本的產生(步驟2)與認定(也就是步驟4的「按 受或拒絕」的準則)。 Metropolis-Hastings演算法如下 1.合X0鳥任何變數樣本 2.第k+1個變數的樣本爲 (YXk) 3.產生一數~U(0,1) 4.如果 us Tk 則Xk+1=y,否則X+1=Xk(即下一個變數仍雜持Xk),其中 (Y)qXkY 7k=m(·x(x)(YXA) g(·-)爲提議機羍函數( proposal distribution)或轉換函數( transition func tion). Markov chain的轉換機羍寫成?(-)=q(|-)r 新樣本來自提議機羍函數q(Y|X),這個函數的選擇當然扮演重大的角色,神奇的是, q(YX)可以是任何型態的函數,其所產生的 Markov chain達穩定狀態後(sta tionary)的分配將是丌()。前面 Metropolis演算法所探用的方式一般稱鳥 random walk chain.也是q(YX)的選擇之一。不過適當的選擇卻關係到收斂( (stationary) 的速度,一般而言,選擇q(YX)的原則不外乎 q(YX)愈按近目標機函數丌(Y)愈好(如範例1所示)。 容易自q(YX)產生樣本 Metropolis演算法是 Metropolis- Hastings演算法的特例,主要的差别在其具對稱 性的提議機羍函數,即q(YX)=q(X1Y),這個假設影響了樣本的接受/拒絕率。2以 上沈 Metropolis演算法步驟2鴦例,新樣本的產生可以表小鳥 y Nq(Y XR=q(Y-XkD=s s~U(0,a) 另外,q(Y|X)=N(X,a2)的常態假設也是常見的特例,因盒q(Y1X)=q(XY) 以下利用簡單的實作問題來瞭解MCMC的運作方式。 範例1:假設目熛函數π()鹪標凖常態N(0,1),寫一支程式利用 Metropolis演算 法產生 Markov chain樣本並測武當提議分配q(YX)分别焉 1.如 Metropolis演算法步驟2的 Random walk q(Y|X)=q(Y-X1) 2.q(YX)=N(X,0.5),N(X,0.1),N(X,10)三種 畫圖觀察 Markov(hain達穩定狀態的情況。 圖1展示了三種常態分配作爲提議分配的表現,不難看岀當提議分配愈接近目標分配 時,其收斂的速度比較理想。睛注意這些常態分配都是以 Markov chain變數ⅹ做 鳥中心點,並非固定。即新旳樣本從「以前一個樣本做篤均數」的常態分配中抽樣。 從圖1也可以觀察到經過足夠艮(多)的時問(樣本)以後,譬如100個以後. Markov hain產生的相依性樣本愈接近從目標函數所產生的。這時便可以利用樣本不均數來 計算母體(日標囪數)均數的期望值F(X),或是進一步的流計量估計F(f(X)。而 前一段被捨棄的100個標本稱 Burn-in samples MATLAB也提供指合「 sample產生樣本,試著使用看看,觀察與自己寫的程式 所產生的結果有何不同。譬如 2 r(Yq(Xky 丌(xk)q(Y|Xk)-m(Xk) 100270040500001 圖1: Metropolis演算法:不同提議分配的收斂表現 sigma=sqrt(.5) tarpdf =@(x)normpdf(x); target pdf proppdf=@(x, y) normpdf(y, x, sigma ) %proposal pdf qlyl r proprnd=@(x)normrnd (x, sigma ) generate samples from proposal pdf nsamples= 1500 z= mhsample(1, nsamples,pdf;. tarpdf, proprnd', proprnd, 'symmetric, 1) histfit(z, 50) 上逃「 sample」指命選擇了對稱性的提議分配,所以第三行對提議分配( propp) 的設定實際上是用不的,其秸果與一般的用法相同。 z= mhsamplc(1, nsamplcs, ' pdf, tarpdf, proprnd, proprnd, proppdf', proppdf ); sample」指爷除輸出樣本外,也可以同時轍出「 Acceptance rate,以便觀察提 議分配的表現,讀者不妨參考線上手册關於 sample的其他選項。 1. 4 Data Augmentation 當後驗機窣密度函數(1)沒有分析解,需要訴諸數值解時, Data. Augmentation的 概念十分合適。如同FM演算法所引入的隱藏變數造成較箐單的概似函數型態,Data Δ augmentation也是試圖在後驗機率的計算上引入隱藏蘷數,讓這個計算變得可能與 簡罩,其觀念如下式 P(ey) 0z, yp(zly ) 從p(θy)到p(θz,y)引進了隱藏變數z使得p(θz.y)變得更谷易分析、計算或是 抽樣,當然最好是一個已知型態的分配。引進了z,卻也同時多出個積分式,當式(4 的積分不容易做到時 Monte carlo積分是不錯的選擇,於是式(4)可以看成 p(By P(az, yp(zy)dx ((6z,y)) Ave(p(ezi, y)) 其中藏變數的資料z;是從所謂的 preditive distribution p(zy)抽樣而來。p(zy) 可以寫成 (zy)=/(z0,y)(y)d 上式意謂著z;的產生通常從p(zθ,y)更為方便,當利用 Monte carlo積分時,積 分變數θ的樣本必須從γ(θy)抽樣而來。這於是形成了「雞生蛋,蛋生雞」的循環遞 迴問題:p(硎y)的計算依靠隱藏變數樣本z,而z;從p(θy)抽樣而交。這個Data augmentation演算法可以寫成下列的兩個步驟 I step( Imputation Step):執行下列兩個步驟m次,產生z1,Z2,…,Zm (1)從p(y)抽出一個樣本θ (2)依步驟(1)產生的6;但,從p(z,y)抽出一個z P step( Posterior Step):更新後驗機率鳥 p(y)= ∑(6z,y) 在 I step中有一個m需要決定, Tanner&wong的論文[2中說明m的選擇可以 是變動的,甚至低到m=1都可以保證迥圈移動的方向是正確的。另外在 P step中 更新的後驗機羍可以被看成是一個混合機羍密度囡數( mixture distribution),而在 L-step中的p(z,y)必須是容易計算或容易抽樣的密度函數,這個演算法扌有實質 的意義。這個演算法被證明凹2在大部分的假設情況下,都能收斂到理論的後驗機羍密 度函數。 範例2:舉著名的 Genetic Linkage Model鳥例,197隻動物以多項分配的分佈屬於 四種類別,分別是 y=(y,y,y3,y4)=(125,18,20,34) 其各項的楼李分別是 6× 61b6 欲估計未知參數θ 最瘟接的作法是訴諸最大概似函數的佔計, L(y 0)= max L(y 6) 0<6<1 其中概似數 Ly0)∝(2+6)12(1-0)0p 參數θ的最大概似估計可以由解概似函數的一次導數爲0,或是輾轉的以EM涣算 法計算可得。另一種作法是更直接的以其件式期望值作估計,即 E(0ly)=/p(ely)de 其中p(0ly)焉後驗機羍。上式可以迤一步寫成 ∫ep(y|O)xr(6)1 6=E(O|y)= π(θ)爲先驗機羍。上式又稱貝氐估計。貝氏估計經常面臨到積分的問題,特別當變 數維度很高時,積分的難度便高岀許多,往往得不到分析解( analytical solution)。以 本題鶯例,假設對先驗機羍一無所知,即π(0)-L.,0≤0≤I,式(10寫成 he2+0)1201-0)3sed (2+0)25(1-6)3834a0 僅是單一變數,這樣的積分也夠嚇人的了 這類的問題譏MCMC找到發揮的舞台, Monte carlo積分透過抽取大量樣本並計算 其吓均來近似期望值的作法,可以避開繁瑱的積分。但是當抽取樣本所憑藉的後驗機 牽不是典型的機羍密度巫數時,以本範例鳥例 p(|y)∝(2+b)12(1-0)034 樣本的產生便是個問題。還好馬可夫鍊( Markov chain)的概念提供一個解套的方 式,以下先利用1.3介紹的兩種MCMC演算法,看看效果如何? Metropolis algorithm 圄2根據式(12)的機羍密度函數以 Metropolis演算法所產生1000個樣本畫出的 直方圖。直方圖的分佈與式(I2)非常的接近。從樣本計算式(9)的積分式,得到的估 計值如圖3中的紅線所示。當懞本數愈大時,其估計值愈穩定,表示從均等分配產生的 Histogram of the randor semple 10 圄2:利用 Metropolis演算法產生的樣本畫出的直方圖 樣本趨近日標機牽酾數(12)。請注意圖3的紅線是所謂的「累加夲均偵,可以使用 MATLAB的指合「 cumsum」進行累加後冉做4均,譬如 y=cmsm(x)./(1:m)%向量x代表1xm的樣本 Metropolis-Hastings Algorithm 當應用 Metropolis- Hastin⑧s演算法時轉換機羍函數的選擇很關鍵。以本題的目標 機羍函數鳥例(如圖4實線),可以選擇常態分配或選擰Beta(56,38)(如圖4虛線)當 作轉換機率囪數,其實Beta(56,34)作鴦q(Y|X)最駑完美,幾乎與目標機率鹵數 致.以此轉換函數產生的樣本與均數展示在圖5。程式如下 125182034 A=56B=38: dr=((x)(2+x)^y(1)*(1-x)^(y(2)+y(3)*^y(4); ppdf =@(a, y) betapdf(r, A, B) nd=o() betarnd(A, B) samples= mhsample(0.5, 1500,, tpdf, proprnd', prnd,'proppdf', ppdf) histfit(samples, 50) 【实例截图】
【核心代码】
标签：
好例子网口号：伸出你的我的手 — 分享！
小贴士
感谢您为本站写下的评论，您的评论对其它用户来说具有重要的参考价值，所以请认真填写。
- 类似“顶”、“沙发”之类没有营养的文字，对勤劳贡献的楼主来说是令人沮丧的反馈信息。
- 相信您也不想看到一排文字/表情墙，所以请不要反馈意义不大的重复字符，也请尽量不要纯表情的回复。
- 提问之前请再仔细看一遍楼主的说明，或许是您遗漏了。
- 请勿到处挖坑绊人、招贴广告。既占空间让人厌烦，又没人会搭理，于人于己都无利。
关于好例子网
本站旨在为广大IT学习爱好者提供一个非营利性互相学习交流分享平台。本站所有资源都可以被免费获取学习研究。本站资源来自网友分享，对搜索内容的合法性不具有预见性、识别性、控制性，仅供学习研究，请务必在下载后24小时内给予删除，不得用于其他任何用途，否则后果自负。基于互联网的特殊性，平台无法对用户传输的作品、信息、内容的权属或合法性、安全性、合规性、真实性、科学性、完整权、有效性等进行实质审查；无论平台是否已进行审查，用户均应自行承担因其传输的作品、信息、内容而可能或已经产生的侵权或权属纠纷等法律责任。本站所有资源不代表本站的观点或立场，基于网友分享，根据中国法律《信息网络传播权保护条例》第二十二与二十三条之规定，若资源存在侵权或相关问题请联系本站客服人员，点此联系我们。关于更多版权及免责申明参见 版权及免责申明
网友评论
我要评论