原文連結：http://tecdat.cn/?p=3772
原文出處：拓端數據部落公眾號
最近我們被客戶要求撰寫關於MCMC的研究報告，包括一些圖形和統計輸出。
創建測試數據
第一步，我們創建一些測試數據，用來擬合我們的模型。我們假設預測變量和因變量之間存在線性關係，所以我們用線性模型並添加一些噪音。
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)))}
斜率參數對數似然曲線
作為說明，代碼的最後幾行繪製了斜率參數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將要處理的實際數量。這個函數被稱為後驗 。同樣，在這裡我們使用加總，因為我們使用對數。
posterior <- function(param){
return ( (param) + prior(param))
}
**MCMC
接下來是Metropolis-Hastings算法。該算法最常見的應用之一（如本例所示）是從貝葉斯統計中的後驗密度中提取樣本。然而，原則上，該算法可用於從任何可積函數中進行採樣。因此，該算法的目的是在參數空間中跳轉，但是以某種方式使得在某一點上的機率與我們採樣的函數成比例（這通常稱為目標函數）。在我們的例子中，這是上面定義的後驗。
這是通過
當我們運行這個算法時，它訪問的參數的分布會收斂到目標分布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）。下一行顯示馬爾可夫鏈參數值。
還有問題嗎？請在下面留言！
最受歡迎的見解
1.matlab使用貝葉斯優化的深度學習
2.matlab貝葉斯隱馬爾可夫hmm模型實現
3.R語言Gibbs抽樣的貝葉斯簡單線性回歸仿真
4.R語言中的block Gibbs吉布斯採樣貝葉斯多元線性回歸
5.R語言中的Stan機率編程MCMC採樣的貝葉斯模型
6.Python用PyMC3實現貝葉斯線性回歸模型
7.R語言使用貝葉斯 層次模型進行空間數據分析
8.R語言隨機搜索變量選擇SSVS估計貝葉斯向量自回歸（BVAR）模型
9.matlab貝葉斯隱馬爾可夫hmm模型實現