R語言實現MCMC中的Metropolis–Hastings算法與吉布斯採樣
原文連接：http://tecdat.cn/?p=3772
建立測試數據
第一步，咱們建立一些測試數據，用來擬合咱們的模型。咱們假設預測變量和因變量之間存在線性關係，因此咱們用線性模型並添加一些噪音。python
trueA <- 5 trueB <- 0 trueSd <- 10 sampleSize <- 31 # 建立獨立的x值 x <- (-(sampleSize-1)/2):((sampleSize-1)/2) # 根據ax + b + N（0，sd）建立因變量 y <- trueA * x + trueB + rnorm(n=sampleSize,mean=0,sd=trueSd) plot(x,y, main="Test Data")
圖算法
定義統計模型
下一步是指定統計模型。咱們已經知道數據是用x和y之間的線性關係y = a * x + b和帶有標準差sd的正態偏差模型N（0，sd）建立的，因此讓咱們使用相同的模型進行擬合，看看若是咱們能夠檢索咱們的原始參數值。編程
從模型中導出似然函數
爲了估計貝葉斯分析中的參數，咱們須要導出咱們想要擬合的模型的似然函數。似然函數是咱們指望觀察到的數據以咱們所看到的模型的參數爲條件發生的機率（密度）。所以，鑑於咱們的線性模型y = b + a*x + N(0,sd)將參數（a, b, sd）做爲輸入，咱們必須返回在這個模型下得到上述測試數據的機率（這聽起來比較複雜，正如你在代碼中看到的，咱們只是計算預測值y = b + a*x與觀察到的y之間的差別，而後咱們必須查找這種誤差發生的機率密度（使用dnorm）。函數
likelihood <- function(param){
a = param\[1\]
b = param\[2\]
sd = param\[3\]
pred = a*x + b
sumll = sum(singlelikelihoods)
(sumll)
}
slopevalues <- function(x){return(likelihood(c(x, trueB, trueSd)))}
斜率參數對數似然曲線post
做爲說明，代碼的最後幾行繪製了斜率參數a的一系列參數值的似然函數。學習
爲何咱們使用對數
您注意到結果是似然函數中機率的對數，這也是我對全部數據點的機率求和的緣由（乘積的對數等於對數之和）。咱們爲何要作這個？ 由於不少小几率乘以的可能性很快就會變得很是小（好比10 ^ -34）。在某些階段，計算機程序存在數字四捨五入的問題。 測試
定義先驗
第二步，與貝葉斯統計中同樣，咱們必須爲每一個參數指定先驗分佈。爲了方便起見，我對全部三個參數使用了均勻分佈和正態分佈。 _無信息先驗一般是1 / sigma（若是你想了解緣由，請看這裏）。 _優化
#先驗分佈
prior <- function(param){
a = param\[1\]
b = param\[2\]
sd = param\[3\]
aprior = (a, min=0, max=10, log = T)
bprior = dnorm(b, sd = 5, log = T)
}
後驗
先驗和似然性的乘積是MCMC將要處理的實際數量。這個函數被稱爲後驗 。一樣，在這裏咱們使用加總，由於咱們使用對數。spa
posterior <- function(param){
return ( (param) + prior(param))
}
MCMC
接下來是Metropolis-Hastings算法。該算法最多見的應用之一（如本例所示）是從貝葉斯統計中的後驗密度中提取樣本。然而，原則上，該算法可用於從任何可積函數中進行採樣。所以，該算法的目的是在參數空間中跳轉，可是以某種方式使得在某一點上的機率與咱們採樣的函數成比例（這一般稱爲目標函數）。在咱們的例子中，這是上面定義的後驗。code
這是經過
- 從隨機參數值開始
- 根據稱爲提議函數的某個機率密度，選擇接近舊值的新參數值
- 以機率p（新）/ p（舊）跳到這個新點，其中p是目標函數，p> 1表示跳躍
當咱們運行這個算法時，它訪問的參數的分佈會收斂到目標分佈p。那麼，讓咱們在R中獲得 ：
########Metropolis算法# ################
proposalfunction <- function(param){
return(rnorm(3,mean = param, sd= c(0.1,0.5,0.3)))
}
run\_metropolis\_MCMC <- function(startvalue, iterations){
for (i in 1:iterations){
if (runif(1) < probab){
chain\[i+1,\] = proposal
}else{
chain\[i+1,\] = chain\[i,\]
}
}
return(chain)
}
chain = run\_metropolis\_MCMC(startvalue, 10000)
burnIn = 5000
acceptance = 1-mean(duplicated(chain\[-(1:burnIn),\]))
使用後驗的對數可能在開始時有點混亂，特別是當您查看計算接受機率的行時（probab = exp（後驗分佈（提議分佈） - 後驗（鏈[i，]）） ）。要理解咱們爲何這樣作，請注意p1 / p2 = exp [log（p1）-log（p2）]。
算法的第一步可能受初始值的誤差，所以一般被丟棄用於進一步分析 。要看的一個有趣的輸出是接受率： 接受標準拒絕提議的頻率是多少？接受率能夠受提議函數的影響：一般，提議越接近，接受率越大。然而，很是高的接受率一般是無益的：這意味着算法「停留」在同一點 。能夠證實，20％到30％的接受率對於典型應用來講是最佳的 。
###概要####################### par(mfrow = c(2,3)) hist( \[-(1:burnIn),1\],nclass=30, , main="Posterior of a", xlab="True value = red line" ) abline(v = mean(chain\[-(1:burnIn),1\])) #進行比較： summary(lm(y~x))
生成的圖應該相似於下圖。您能夠看到咱們或多或少地檢索了用於建立數據的原始參數，而且您還看到咱們得到了一個圍繞最高後驗值的特定區域，這些區域也有一些數據支持，這是貝葉斯至關於置信度的區間。
所獲得的圖應該看起來像下面的圖。你看到咱們檢索到了或多或少用於建立數據的原始參數，你還看到咱們在最高後驗值周圍獲得了必定的區域，這些後驗值也有一些數據，這至關於貝葉斯的置信區間。
圖：上排顯示斜率（a）的後驗估計，截距（b）和偏差的標準誤差（sd）。下一行顯示馬爾可夫鏈參數值。
還有問題嗎？請在下面留言！
最受歡迎的看法
4.R語言中的block Gibbs吉布斯採樣貝葉斯多元線性迴歸