R語言實(shí)現(xiàn)MCMC中的Metropolis–Hastings算法與吉布斯采樣|附代碼數(shù)據(jù)
原文鏈接：http://tecdat.cn/?p=3772
原文出處：拓端數(shù)據(jù)部落公眾號(hào)
?最近我們被客戶要求撰寫關(guān)于MCMC的研究報(bào)告，包括一些圖形和統(tǒng)計(jì)輸出。
創(chuàng)建測試數(shù)據(jù)
第一步，我們創(chuàng)建一些測試數(shù)據(jù)，用來擬合我們的模型。我們假設(shè)預(yù)測變量和因變量之間存在線性關(guān)系，所以我們用線性模型并添加一些噪音。
trueA <- 5trueB <- 0trueSd <- 10sampleSize <- 31?# 創(chuàng)建獨(dú)立的x值x <- (-(sampleSize-1)/2):((sampleSize-1)/2)# 根據(jù)ax + b + N（0，sd）創(chuàng)建因變量y <-? trueA * x + trueB + rnorm(n=sampleSize,mean=0,sd=trueSd)?plot(x,y, main="Test Data")
圖
定義統(tǒng)計(jì)模型
下一步是指定統(tǒng)計(jì)模型。我們已經(jīng)知道數(shù)據(jù)是用x和y之間的線性關(guān)系y = a * x + b和帶有標(biāo)準(zhǔn)差sd的正態(tài)誤差模型N（0，sd）創(chuàng)建的，所以讓我們使用相同的模型進(jìn)行擬合，看看如果我們可以檢索我們的原始參數(shù)值。
從模型中導(dǎo)出似然函數(shù)
為了估計(jì)貝葉斯分析中的參數(shù)，我們需要導(dǎo)出我們想要擬合的模型的似然函數(shù)。似然函數(shù)是我們期望觀察到的數(shù)據(jù)以我們所看到的模型的參數(shù)為條件發(fā)生的概率（密度）。因此，鑒于我們的線性模型y = b + a*x + N(0,sd)將參數(shù)（a, b, sd）作為輸入，我們必須返回在這個(gè)模型下獲得上述測試數(shù)據(jù)的概率（這聽起來比較復(fù)雜，正如你在代碼中看到的，我們只是計(jì)算預(yù)測值y = b + a*x與觀察到的y之間的差異，然后我們必須查找這種偏差發(fā)生的概率密度（使用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)))}
斜率參數(shù)對(duì)數(shù)似然曲線
作為說明，代碼的最后幾行繪制了斜率參數(shù)a的一系列參數(shù)值的似然函數(shù)。
為什么我們使用對(duì)數(shù)
您注意到結(jié)果是似然函數(shù)中概率的對(duì)數(shù)，這也是我對(duì)所有數(shù)據(jù)點(diǎn)的概率求和的原因（乘積的對(duì)數(shù)等于對(duì)數(shù)之和）。我們?yōu)槭裁匆鲞@個(gè)？ 因?yàn)楹芏嘈「怕食艘缘目赡苄院芸炀蜁?huì)變得非常?。ū热?0 ^ -34）。在某些階段，計(jì)算機(jī)程序存在數(shù)字四舍五入的問題。?
定義先驗(yàn)
第二步，與貝葉斯統(tǒng)計(jì)中一樣，我們必須為每個(gè)參數(shù)指定先驗(yàn)分布。為了方便起見，我對(duì)所有三個(gè)參數(shù)使用了均勻分布和正態(tài)分布。?無信息先驗(yàn)通常是1 / sigma（如果你想了解原因，請(qǐng)看這里）。
#先驗(yàn)分布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)?}
****后驗(yàn)
先驗(yàn)和似然性的乘積是MCMC將要處理的實(shí)際數(shù)量。這個(gè)函數(shù)被稱為后驗(yàn)?。同樣，在這里我們使用加總，因?yàn)槲覀兪褂脤?duì)數(shù)。
posterior <- function(param){???return (?(param) + prior(param))}
?
**MCMC
接下來是Metropolis-Hastings算法。該算法最常見的應(yīng)用之一（如本例所示）是從貝葉斯統(tǒng)計(jì)中的后驗(yàn)密度中提取樣本。然而，原則上，該算法可用于從任何可積函數(shù)中進(jìn)行采樣。因此，該算法的目的是在參數(shù)空間中跳轉(zhuǎn)，但是以某種方式使得在某一點(diǎn)上的概率與我們采樣的函數(shù)成比例（這通常稱為目標(biāo)函數(shù)）。在我們的例子中，這是上面定義的后驗(yàn)。
這是通過
從隨機(jī)參數(shù)值開始
根據(jù)稱為提議函數(shù)的某個(gè)概率密度，選擇接近舊值的新參數(shù)值
以概率p（新）/ p（舊）跳到這個(gè)新點(diǎn)，其中p是目標(biāo)函數(shù)，p> 1表示跳躍
當(dāng)我們運(yùn)行這個(gè)算法時(shí)，它訪問的參數(shù)的分布會(huì)收斂到目標(biāo)分布p。那么，讓我們?cè)赗中得到 ：
########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 = 5000acceptance = 1-mean(duplicated(chain[-(1:burnIn),]))
使用后驗(yàn)的對(duì)數(shù)可能在開始時(shí)有點(diǎn)混亂，特別是當(dāng)您查看計(jì)算接受概率的行時(shí)（probab = exp（后驗(yàn)分布（提議分布） - 后驗(yàn)（鏈[i，]）） ）。要理解我們?yōu)槭裁催@樣做，請(qǐng)注意p1 / p2 = exp [log（p1）-log（p2）]。
算法的第一步可能受初始值的偏差，因此通常被丟棄用于進(jìn)一步分析 。要看的一個(gè)有趣的輸出是接受率： 接受標(biāo)準(zhǔn)拒絕提議的頻率是多少？接受率可以受提議函數(shù)的影響：通常，提議越接近，接受率越大。然而，非常高的接受率通常是無益的：這意味著算法“停留”在同一點(diǎn) 。可以證明，20％到30％的接受率對(duì)于典型應(yīng)用來說是最佳的 。
?
###概要#######################?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]))??#進(jìn)行比較：summary(lm(y~x))
所得到的圖應(yīng)該看起來像下面的圖。你看到我們檢索到了或多或少用于創(chuàng)建數(shù)據(jù)的原始參數(shù)，你還看到我們?cè)谧罡吆篁?yàn)值周圍得到了一定的區(qū)域，這些后驗(yàn)值也有一些數(shù)據(jù)，這相當(dāng)于貝葉斯的置信區(qū)間。
圖：?上排顯示斜率（a）的后驗(yàn)估計(jì)，截距（b）和誤差的標(biāo)準(zhǔn)偏差（sd）。下一行顯示馬爾可夫鏈參數(shù)值。
?
還有問題嗎？請(qǐng)?jiān)谙旅媪粞裕?/p>
最受歡迎的見解
1.matlab使用貝葉斯優(yōu)化的深度學(xué)習(xí)
2.matlab貝葉斯隱馬爾可夫hmm模型實(shí)現(xiàn)
3.R語言Gibbs抽樣的貝葉斯簡單線性回歸仿真
4.R語言中的block Gibbs吉布斯采樣貝葉斯多元線性回歸
5.R語言中的Stan概率編程MCMC采樣的貝葉斯模型
6.Python用PyMC3實(shí)現(xiàn)貝葉斯線性回歸模型
7.R語言使用貝葉斯 層次模型進(jìn)行空間數(shù)據(jù)分析
8.R語言隨機(jī)搜索變量選擇SSVS估計(jì)貝葉斯向量自回歸（BVAR）模型
9.matlab貝葉斯隱馬爾可夫hmm模型實(shí)現(xiàn)
