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

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

一、Metropolis Hasting

1、算法理解

我們的目標(biāo)是對Target Distribution進行抽樣，首先，我們要引入一條具有平穩(wěn)分布的馬氏鏈，這條馬氏鏈?zhǔn)諗康钠椒€(wěn)分布我們稱為Proposal Distribution，而這條馬氏鏈的表現(xiàn)形式是概率轉(zhuǎn)移矩陣

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

如何根據(jù)這條馬氏鏈求得目標(biāo)分布呢？這里由馬氏鏈的細致平穩(wěn)性引入。

是目標(biāo)分布下的隨機變量， 是proposal distribution下的隨機變量。

（1.1）成立。（由馬氏鏈的細致平穩(wěn)性得到，表示i，j狀態(tài)之間的能量轉(zhuǎn)換相等）

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

（1.3）

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

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

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

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

2、proposal distribution的選取

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

3、共軛的正態(tài)分布示例

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

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

。

找到下一個狀態(tài) 。這里proposal distribution設(shè)為正態(tài)分布。生成 , .

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

接受率的概率實現(xiàn)。如果接受， ,否則， 。

import

二、Gibbs sampling

1、算法理解

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

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

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

2、示例

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

的泊松分布，每個雞蛋成功孵出小雞的比例為p。p未知，其先驗分布服從beta分布。 .參數(shù) 已知。我們的觀測數(shù)據(jù)只有每天孵出的小雞個數(shù) , 屬于隱變量，觀測不到。如何通過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

創(chuàng)作挑戰(zhàn)賽新人創(chuàng)作獎勵來咯，堅持創(chuàng)作打卡瓜分現(xiàn)金大獎

總結(jié)

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

如果覺得生活随笔網(wǎng)站內(nèi)容還不錯，歡迎將生活随笔推薦給好友。