分布是一系列數字的規律組合。如果在收集了歷史中的幾百個數據后，我想知道這群數據背后的發射機制是什么，那么就得去尋找這個分布。當然這里的重點不是尋找分布，而是在已知分布的情況下，如何模擬這個機制發射出來的一系列數字呢？

MCMC（Markov Chain Monte Carlo）是馬爾科夫鏈下的蒙特卡洛方法，因為馬爾科夫鏈在滿足某些條件下具有平穩分布，如果能夠將平穩分布與目標分布聯系起來，那么就可以達到對目標分布進行抽樣的目的。這里主要介紹的是Metropolis Hasting 算法和Gibbs sampling 算法。

一、Metropolis Hasting

1、算法理解

我們的目標是對Target Distribution進行抽樣，首先，我們要引入一條具有平穩分布的馬氏鏈，這條馬氏鏈收斂的平穩分布我們稱為Proposal Distribution，而這條馬氏鏈的表現形式是概率轉移矩陣

，狀態空間 ，狀態空間也即是Proposal distribution的所有可能取值集合。

如何根據這條馬氏鏈求得目標分布呢？這里由馬氏鏈的細致平穩性引入。

是目標分布下的隨機變量， 是proposal distribution下的隨機變量。

（1.1）成立。（由馬氏鏈的細致平穩性得到，表示i，j狀態之間的能量轉換相等）

（1.2）（因為 與 是兩個不同的分布）

（1.3）

為了使（1.2）式成立，所以引入了接受率

。其中 即將不等式的左右兩邊互相相乘，即可得到式子（1.3）。接受率 表示是否決定抽取下一個樣本（i.e., 接受樣本j），因此我們需要將這個概率實現，因為在實際抽樣過程中，決定抽樣和不抽樣是一個二元過程，而不是說以多大的概率決定抽樣。這個概率實現可以用伯努利分布，也可以用均勻分布：當均勻分布下的數值小于接受率時，決定抽樣，反之不抽樣。

以上就是Metropolis抽樣方法的全部內容了，而Metropolis hasting 算法則對接受率做了一點改進。當接受率太小的時候，我們很難從當前的樣本值跳到其他狀態，所以對

進行了擴大。將 中的較大值擴充到1（即一定會抽取下個樣本），另外一個值等比例擴大。經過計算可以得到表達式 。在計算接受率的過程中，我們就會發現，目標分布的常數項被抵消了，也去除了歸一化的過程。

2、proposal distribution的選取

當proposal distribution與目標分布越靠近時，抽取的樣本也就越合理。但是proposal distribution下的馬氏鏈如何確定，兩個分布的距離如何衡量，這些也都是可以繼續探討也需要權衡的問題。

3、共軛的正態分布示例

已知， 未知，在貝葉斯統計下， 是一個隨機變量，其先驗分布為 已知。如何利用Metropolis-Hasting算法，在觀察數據Y下求得 后驗分布得期望和方差？

我們用M-H抽樣算法來檢驗上面得后驗分布是否準確。即在已知得各參數和觀測值y下抽出一系列的

。

找到下一個狀態 。這里proposal distribution設為正態分布。生成 , .

接受率 .其中， 為隨機變量 的概率密度。

接受率的概率實現。如果接受， ,否則， 。

import

二、Gibbs sampling

1、算法理解

Gibbs sampling適用于高維分布的抽樣問題。在M-H抽樣算法的基礎上，如果我們能夠比較容易的得到條件分布，那么就可以通過固定其他維度，一次只對一個維度上的條件分布抽樣的方法進行全局抽樣。

Gibbs sampling里的接受率恒為1。舉例說明，

兩個樣本點滿足馬氏鏈的細致平穩條件。因為 其中， 表示從A點轉移到B點的轉移概率。所以在二維的分布中，可以得到從任意一個點轉移到另外一個點都是平穩的，限制是每次變換只能轉移一個維度。二維轉移圖可如下所示。

2、示例

一只雞每天會下N個蛋,N服從參數為

的泊松分布，每個雞蛋成功孵出小雞的比例為p。p未知，其先驗分布服從beta分布。 .參數 已知。我們的觀測數據只有每天孵出的小雞個數 , 屬于隱變量，觀測不到。如何通過Gibbs Sampling 方法找到p的后驗期望呢？

在不引入隨機變量N的時候，后驗分布比較麻煩。引入N后，可得，

通過迭代，即可得到p,N的抽樣值。

x, lambda1, a, b = 7, 10, 1, 1 niter = 10000 p = [0 for i in range(niter)] N = [0 for i in range(niter)]#初始值 p[0] = 0.5 N[0] = 2*x for i in range(1,niter):p[i] = random.betavariate(x+a, N[i-1]-x+b)N[i] = x + np.random.poisson(lambda1*(1-p[i-1]))plt.hist(x = p, bins = 100,normed=True) plt.hist(x = N, bins = 100,normed=True) plt.show()

^[1][2][3]

參考

^https://blog.csdn.net/lin360580306/article/details/51240398

^https://www.jianshu.com/p/dc7624cf6adb

^Introduction to Probability--Joseph K.Blitzstein

創作挑戰賽新人創作獎勵來咯，堅持創作打卡瓜分現金大獎

總結

以上是生活随笔為你收集整理的t分布f分布与样本均值抽样分布_分布模拟1——MCMC抽样方法的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。