拓端tecdat|R語(yǔ)言貝葉斯推斷與MCMC：實(shí)現(xiàn)Metropolis
原文鏈接：http://tecdat.cn/?p=21545
原文出處：拓端數(shù)據(jù)部落公眾號(hào)
示例1：使用MCMC的指數(shù)分布采樣
任何MCMC方案的目標(biāo)都是從“目標(biāo)”分布產(chǎn)生樣本。在這種情況下，我們將使用平均值為1的指數(shù)分布作為我們的目標(biāo)分布。所以我們從定義目標(biāo)密度開(kāi)始：
target = function(x){if(x<0){return(0)}else {return( exp(-x))}}
定義了函數(shù)之后，我們現(xiàn)在可以用它來(lái)計(jì)算幾個(gè)值（只是為了說(shuō)明函數(shù)的概念）：
target(1)[1] 0.3678794target(-1)[1] 0
接下來(lái)，我們將規(guī)劃一個(gè)Metropolis-Hastings方案，從與目標(biāo)成比例的分布中進(jìn)行抽樣
x[1] = 3 ? ? #這只是一個(gè)起始值，我設(shè)置為3for(i in 2:1000){A = target(proposedx)/target(currentx)if(runif(1)<A){x[i] = proposedx ? ? ? # 接受概率min（1，a）} else {x[i] = currentx ? ? ? ?#否則“拒絕”行動(dòng)，保持原樣}
注意，x是馬爾可夫鏈的實(shí)現(xiàn)。我們可以畫(huà)幾個(gè)x的圖：
我們可以將其封裝在一個(gè)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 ? ? ? ?# 否則“拒絕”行動(dòng)，保持原樣}
現(xiàn)在我們將運(yùn)行MCMC方案3次，看看結(jié)果有多相似：
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個(gè)圖形放在一個(gè)頁(yè)面上hist(z1,breaks=seq(0,maxz,length=20))
練習(xí)
使用函數(shù)easyMCMC了解以下內(nèi)容：
不同的起始值如何影響MCMC方案？
較大/較小的提案標(biāo)準(zhǔn)差有什么影響？
嘗試將目標(biāo)函數(shù)更改為
target = function(x){return((x>0 & x <1) + (x>2 & x<3))}
這個(gè)目標(biāo)看起來(lái)像什么？如果提議sd太小怎么辦？（例如，嘗試1和0.1）
例2：估計(jì)等位基因頻率
在對(duì)雙等位基因座的基因型（如具有AA和AA等位基因的基因座）進(jìn)行建模時(shí)，一個(gè)標(biāo)準(zhǔn)的假設(shè)是群體是“隨機(jī)”的。這意味著如果p是等位基因AA的頻率，那么基因型
和
將分別具有頻率
和
。
p一個(gè)簡(jiǎn)單的先驗(yàn)是假設(shè)它在[0,1]上是均勻的。 假設(shè)我們抽樣n個(gè)個(gè)體，觀察
基因型
、
基因型
和
基因型
。
下面的R代碼給出了一個(gè)簡(jiǎn)短的MCMC例程，可以從p的后驗(yàn)分布中進(jìn)行采樣。請(qǐng)嘗試遍歷該代碼，看看它是如何工作的。
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 ? ? ? ?# 否則“拒絕”行動(dòng)，保持原樣}
對(duì)
運(yùn)行此樣本。
現(xiàn)在用一些R代碼來(lái)比較后驗(yàn)樣本和理論后驗(yàn)樣本（在這種情況下可以通過(guò)分析獲得；因?yàn)槲覀冇^察到121個(gè)As和79個(gè)as，在200個(gè)樣本中，p的后驗(yàn)樣本是β（121+1,79+1）。
hist(z,prob=T)lines(x,dbeta(x,122, 80)) ?# 在直方圖上疊加β密度
您也可能希望將前5000 z的值丟棄為“burnin”（預(yù)燒期）。這里有一種方法，在R中僅選擇最后5000 z
hist(z[5001:10000])
練習(xí)
研究起點(diǎn)和提案標(biāo)準(zhǔn)偏差如何影響算法的收斂性。
例3：估計(jì)等位基因頻率和近交系數(shù)
一個(gè)復(fù)雜一點(diǎn)的替代方法是假設(shè)人們有一種傾向，即人們會(huì)與比“隨機(jī)”關(guān)系更密切的其他人近交（例如，在地理結(jié)構(gòu)上的人口中可能會(huì)發(fā)生這種情況）。一個(gè)簡(jiǎn)單的方法是引入一個(gè)額外的參數(shù)，即“近親繁殖系數(shù)”f，并假設(shè)
和
基因型有頻率
和
。
在大多數(shù)情況下，將f作為種群特征來(lái)對(duì)待是很自然的，因此假設(shè)f在各個(gè)位點(diǎn)上是恒定的。
請(qǐng)注意，f和p都被約束為介于0和1之間（包括0和1）。對(duì)于這兩個(gè)參數(shù)中的每一個(gè)，一個(gè)簡(jiǎn)單的先驗(yàn)是假設(shè)它們?cè)赱0,1]上是獨(dú)立的。 假設(shè)我們抽樣n個(gè)個(gè)體，觀察
基因型
、
基因型
和
基因型
。
練習(xí)：
編寫(xiě)一個(gè)短的MCMC程序，從f和p的聯(lián)合分布中取樣。
sampler = function(){f[1] = fstartvalp[1] = pstartvalfor(i in 2:niter){currentf = f[i-1]currentp = p[i-1]newf = currentf +newp = currentp +}return(list(f=f,p=p)) # 返回一個(gè)包含兩個(gè)名為f和p的元素的“l(fā)ist”}
使用此樣本獲得f和p的點(diǎn)估計(jì)（例如，使用后驗(yàn)平均數(shù)）和f和p的區(qū)間估計(jì)（例如，90%后驗(yàn)置信區(qū)間），數(shù)據(jù)：
附錄：GIBBS采樣
您也可以用Gibbs采樣器解決這個(gè)問(wèn)題?
為此，您將想要使用以下“潛在變量”表示模型：
將zi相加得到與上述相同的模型:
?
最受歡迎的見(jiàn)解
1.使用R語(yǔ)言進(jìn)行METROPLIS-IN-GIBBS采樣和MCMC運(yùn)行
2.R語(yǔ)言中的Stan概率編程MCMC采樣的貝葉斯模型
3.R語(yǔ)言實(shí)現(xiàn)MCMC中的Metropolis–Hastings算法與吉布斯采樣
4.R語(yǔ)言BUGS JAGS貝葉斯分析 馬爾科夫鏈蒙特卡洛方法（MCMC）采樣
5.R語(yǔ)言中的block Gibbs吉布斯采樣貝葉斯多元線性回歸
6.R語(yǔ)言Gibbs抽樣的貝葉斯簡(jiǎn)單線性回歸仿真分析
7.R語(yǔ)言用Rcpp加速M(fèi)etropolis-Hastings抽樣估計(jì)貝葉斯邏輯回歸模型的參數(shù)
8.R語(yǔ)言使用Metropolis- Hasting抽樣算法進(jìn)行邏輯回歸
9.R語(yǔ)言中基于混合數(shù)據(jù)抽樣(MIDAS)回歸的HAR-RV模型預(yù)測(cè)GDP增長(zhǎng)
