原文鏈接：http://tecdat.cn/?p=21545

原文出處：拓端數(shù)據(jù)部落公眾號

? 最近我們被客戶要求撰寫關于貝葉斯推斷的研究報告，包括一些圖形和統(tǒng)計輸出。

示例1：使用MCMC的指數(shù)分布采樣

任何MCMC方案的目標都是從“目標”分布產(chǎn)生樣本。在這種情況下，我們將使用平均值為1的指數(shù)分布作為我們的目標分布。所以我們從定義目標密度開始：

target = function(x){ ?if(x<0){ ? ?return(0)} ?else { ? ?return( exp(-x)) ?}}

定義了函數(shù)之后，我們現(xiàn)在可以用它來計算幾個值（只是為了說明函數(shù)的概念）：

target(1)

[1] 0.3678794

target(-1)

[1] 0

接下來，我們將規(guī)劃一個Metropolis-Hastings方案，從與目標成比例的分布中進行抽樣

x[1] = 3 ? ? #這只是一個起始值，我設置為3for(i in 2:1000){ ?A = target(proposedx)/target(currentx) ?if(runif(1)<A){ ? ?x[i] = proposedx ? ? ? # 接受概率min（1，a） ?} else { ? ?x[i] = currentx ? ? ? ?#否則“拒絕”行動，保持原樣 ?}

注意，x是馬爾可夫鏈的實現(xiàn)。我們可以畫幾個x的圖：

?

?

?

?

我們可以將其封裝在一個mcmc函數(shù)中，以使代碼更整潔，這樣更改起始值和提議分布更容易

?for(i in 2:niter){ ? ?currentx = x[i-1] ? ?proposedx = rnorm(1,mean=currentx,sd=proposalsd) ? ?A = target(proposedx)/target(currentx) ? ?if(runif(1)<A){ ? ? ?x[i] = proposedx ? ? ? # 接受概率min（1，a） ? ?} else { ? ? ?x[i] = currentx ? ? ? ?# 否則“拒絕”行動，保持原樣 ? ?}

現(xiàn)在我們將運行MCMC方案3次，看看結果有多相似：

z1=MCMC(1000,3,1)z2=MCMC(1000,3,1)z3=MCMC(1000,3,1)plot(z1,type="l")

?

?

par(mfcol=c(3,1)) #告訴R將3個圖形放在一個頁面上hist(z1,breaks=seq(0,maxz,length=20))

?

?

練習

使用函數(shù)easyMCMC了解以下內容：

1. 不同的起始值如何影響MCMC方案？

2. 較大/較小的提案標準差有什么影響？

3. 嘗試將目標函數(shù)更改為

target = function(x){ ? ?return((x>0 & x <1) + (x>2 & x<3))}

這個目標看起來像什么？如果提議sd太小怎么辦？（例如，嘗試1和0.1）

例2：估計等位基因頻率

在對雙等位基因座的基因型（如具有AA和AA等位基因的基因座）進行建模時，一個標準的假設是群體是“隨機”的。這意味著如果p是等位基因AA的頻率，那么基因型

和

將分別具有頻率

和

。

p一個簡單的先驗是假設它在[0,1]上是均勻的。 假設我們抽樣n個個體，觀察

基因型

、

基因型

和

基因型

。

下面的R代碼給出了一個簡短的MCMC例程，可以從p的后驗分布中進行采樣。請嘗試遍歷該代碼，看看它是如何工作的。

prior = function(p){ ?if((p<0) || (p>1)){ ?# || 這里意思是“或” ? ?return(0)} ?else{ ? ?return(1)}}likelihood = function(p, nAA, nAa, naa){ ?return(p^(2*nAA) * (2*p*(1-p))^nAa * (1-p)^(2*naa))}psampler = function(nAA, nAa, naa){ ?for(i in 2:niter){ ? ?if(runif(1)<A){ ? ? ?p[i] = newp ? ? ? # 接受概率min（1，a） ? ?} else { ? ? ?p[i] = currentp ? ? ? ?# 否則“拒絕”行動，保持原樣 ? ?}

對

運行此樣本。

現(xiàn)在用一些R代碼來比較后驗樣本和理論后驗樣本（在這種情況下可以通過分析獲得；因為我們觀察到121個As和79個as，在200個樣本中，p的后驗樣本是β（121+1,79+1）。

hist(z,prob=T)lines(x,dbeta(x,122, 80)) ?# 在直方圖上疊加β密度

您也可能希望將前5000 z的值丟棄為“burnin”（預燒期）。這里有一種方法，在R中僅選擇最后5000 z

hist(z[5001:10000])

?

?

練習

研究起點和提案標準偏差如何影響算法的收斂性。

例3：估計等位基因頻率和近交系數(shù)

一個復雜一點的替代方法是假設人們有一種傾向，即人們會與比“隨機”關系更密切的其他人近交（例如，在地理結構上的人口中可能會發(fā)生這種情況）。一個簡單的方法是引入一個額外的參數(shù)，即“近親繁殖系數(shù)”f，并假設

和

基因型有頻率

和

。

在大多數(shù)情況下，將f作為種群特征來對待是很自然的，因此假設f在各個位點上是恒定的。

請注意，f和p都被約束為介于0和1之間（包括0和1）。對于這兩個參數(shù)中的每一個，一個簡單的先驗是假設它們在[0,1]上是獨立的。 假設我們抽樣n個個體，觀察

基因型

、

基因型

和

基因型

。

練習：

• 編寫一個短的MCMC程序，從f和p的聯(lián)合分布中取樣。

sampler = function(){ ?f[1] = fstartval ?p[1] = pstartval ?for(i in 2:niter){ ? ?currentf = f[i-1] ? ?currentp = p[i-1] ? ?newf = currentf + ? ?newp = currentp + ? ?} ?return(list(f=f,p=p)) # 返回一個包含兩個名為f和p的元素的“l(fā)ist”}

• 使用此樣本獲得f和p的點估計（例如，使用后驗平均數(shù)）和f和p的區(qū)間估計（例如，90%后驗置信區(qū)間），數(shù)據(jù)：

?

附錄：GIBBS采樣

您也可以用Gibbs采樣器解決這個問題?

為此，您將想要使用以下“潛在變量”表示模型：

將zi相加得到與上述相同的模型:

?

最受歡迎的見解

1.使用R語言進行METROPLIS-IN-GIBBS采樣和MCMC運行

2.R語言中的Stan概率編程MCMC采樣的貝葉斯模型

3.R語言實現(xiàn)MCMC中的Metropolis–Hastings算法與吉布斯采樣

4.R語言BUGS JAGS貝葉斯分析 馬爾科夫鏈蒙特卡洛方法（MCMC）采樣

5.R語言中的block Gibbs吉布斯采樣貝葉斯多元線性回歸

6.R語言Gibbs抽樣的貝葉斯簡單線性回歸仿真分析

7.R語言用Rcpp加速Metropolis-Hastings抽樣估計貝葉斯邏輯回歸模型的參數(shù)

8.R語言使用Metropolis- Hasting抽樣算法進行邏輯回歸

9.R語言中基于混合數(shù)據(jù)抽樣(MIDAS)回歸的HAR-RV模型預測GDP增長