拓端tecdat|R語言BUGS/JAGS貝葉斯分析: 馬爾科夫鏈蒙特卡洛方法（MCMC）采樣

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

 

 

馬爾科夫鏈蒙特卡洛方法

在許多情況下，我們沒有足夠的計算能力評估空間中所有n維像素的后驗概率 。在這些情況下，我們傾向於利用稱為Markov-Chain Monte Carlo 算法的程序 。此方法使用參數空間中的隨機跳躍來（最終）確定后驗分布。MCMC的關鍵如下：

跳躍概率的比例與后驗概率的比例成正比。

跳躍概率可以表征為：

概率（跳躍）*概率（接受）

從長遠來看，該鏈將花費大量時間在參數空間的高概率部分，從而實質上捕獲了后驗分布。有了足夠的跳躍，長期分布將與聯合后驗概率分布匹配。

MCMC本質上是一種特殊類型的隨機數生成器，旨在從難以描述（例如，多元，分層）的概率分布中采樣。在許多/大多數情況下，后驗分布是很難描述的概率分布。

MCMC使您可以從實際上不可能完全定義的概率分布中進行采樣！

令人驚訝的是，MCMC的核心並不難於描述或實施。讓我們看一個簡單的MCMC算法。

Metropolis-Hastings算法

該算法與模擬退火算法非常相似。

MH算法可以表示為：

Prob（acceptB | A）= min（1，Posterior（B）Posterior（A）⋅Prob（b→a）Prob（a→b））

請注意，從本質上講，這與“ Metropolis”模擬退火算法相同，后驗概率代替了概率，並且 k 參數設置為1。

 

二元正態例子

 

請記住，MCMC采樣器只是隨機數生成器的一種。我們可以使用Metropolis-Hastings采樣器來開發自己的隨機數生成器，生成進行簡單的已知分布。在此示例中，我們使用MH采樣器從標准雙變量正態概率分布生成隨機數。

對於這個簡單的示例，我們不需要MCMC采樣器。一種實現方法是使用以下代碼，該代碼從具有相關參數ρ的雙變量標准正態分布中繪制並可視化任意數量的獨立樣本。

1.  
   #################
2.  
   ＃MCMC采樣的簡單示例
3.  
   #################
4.  
    
5.  
    
6.  
   #########
7.  
   # ＃首先，讓我們構建一個從雙變量標准正態分布生成隨機數的函數
8.  
    
9.  
   rbvn<-function (n, rho) #用於從二元標准正態分布中提取任意數量的獨立樣本。
10.  
    {
11.  
    x <- rnorm(n, 0, 1)
12.  
    y <- rnorm(n, rho * x, sqrt(1 - rho^2))
13.  
    cbind(x, y)
14.  
    }
15.  
     
16.  
    #########
17.  
    # 現在，從該分布圖中繪制隨機抽樣
18.  
    bvn<-rbvn(10000,0.98)
19.  
    par(mfrow=c(3,2))
20.  
    plot(bvn,col=1:10000

 

1.  
   ###############
2.  
   # ＃Metropolis-Hastings雙變量正態采樣器的實現...
3.  
    
4.  
   library(mvtnorm) # 加載一個包，該包使我們能夠計算mv正態分布的概率密度
5.  
    
6.  
   metropoli<- function (n, rho=0.98){ # 雙變量隨機數生成器的MCMC采樣器實現
7.  
   mat <- matrix(ncol = 2, nrow = n) # 用於存儲隨機樣本的矩陣
8.  
   x <- 0 # 所有參數的初始值
9.  
    
10.  
     
11.  
    prev <- dmvnorm(c(x,y),mean=c(0,0),sig
12.  
    # 起始位置分布的概率密度
13.  
    mat[1, ] <- c(x, y) # 初始化馬爾可夫鏈
14.  
     
15.  
    newx <- rnorm(1,x,0.5) # 進行跳轉
16.  
     
17.  
     
18.  
    newprob <- dmvnorm(c(newx,newy),sigma =
19.  
    # 評估跳轉
20.  
    ratio <- newprob/prev # 計算舊位置（跳出）和建議位置（跳到）的概率之比。
21.  
     
22.  
    prob.accept <- min(1,ratio) # 決定接受新跳躍的概率！
23.  
     
24.  
     
25.  
    if(rand<=prob.accept){
26.  
    x=newx;y=newy # 將x和y設置為新位置
27.  
    mat[counter,] <- c(x,y) # 將其存儲在存儲陣列中
28.  
     
29.  
     
30.  
    prev <- newprob # 准備下一次迭代
31.  
     

然后，我們可以使用MH采樣器從該已知分布中獲取隨機樣本…

1.  
   ###########
2.  
   # 測試新的M-H采樣器
3.  
    
4.  
   bvn<-metropoli(10000,0.98)
5.  
   par(mfrow=c(3,2))
6.  
   plot(bvn,col=1:10000)
7.  
   plot(bvn,type=

 

 

讓我們嘗試解決一個問題。

MCMC對粘液瘤病進行調查

1.  
   ############
2.  
   #黏液病示例的MCMC實現
3.  
   ############
4.  
    

1.  
    
2.  
    
3.  
   subset(MyxDat,grade==1
4.  
    

1.  
   ## grade day titer
2.  
   ## 1 1 2 5.207
3.  
   ## 2 1 2 5.734
4.  
   ## 3 1 2 6.613
5.  
   ## 4 1 3 5.997
6.  
   ## 5 1 3 6.612
7.  
   ## 6 1 3 6.810

選擇使用Gamma分布。這是經驗分布：

1.  
   ###########
2.  
   # 第100次可視化粘液病數據
3.  
   hist(Myx$titer,freq=FALSE)

我們需要估算最適合此經驗分布的伽馬速率和形狀參數。這是適合此分布的Gamma的一個示例：

1.  
   #########
2.  
   # ...覆蓋生成模型的數據（伽瑪分布）
3.  
    
4.  
    
5.  
   curve(dgamma(x,shape=40,scale=0.15),add=T,col="red")

 

二維（對數）似然面：

1.  
   ##############
2.  
   # 定義二維參數空間
3.  
   ##############
4.  
    
5.  
   shapevec <- seq(3,100,by=0.1)
6.  
   scalevec <- seq(0.01,0.5,by=0.001)
7.  
    
8.  
   ##############
9.  
   # ＃定義參數空間內此網格上的似然面
10.  
    ##############
11.  
     
12.  
    GammaLogLikelihoodFunction <- function(par
13.  
     
14.  
    }
15.  
    surface2D <- matrix(nrow=length(shapevec),ncol=length(scalevec)) #初始化存儲變量
16.  
     
17.  
    newparams <- c(sha
18.  
     
19.  
    surface2D[i,j] <- GammaLogLikelihoodFunction(newparams)
20.  
     
21.  
     
22.  
     
23.  
    ############
24.  
    # 可視化似然面
25.  
    ############
26.  
     
27.  
     
28.  
    contour(x=shapevec,y=scalevec,z=surface2D,levels=c(-30,-40,-80,-500),add=T)

這是MH算法的實現，用於找到后驗分布！

首先，我們需要一個似然函數–這次，我們將返回真實概率–而不是對數轉換的概率

1.  
   ############
2.  
   #編寫非對數轉換的似然函數
3.  
    
4.  
   GammaLike- function(params){
5.  
   prod(dgamma(Myx$titer,shape=params['shape']
6.  
    
7.  
   params <- c(shape=40,

1.  
   ## shape scale
2.  
   ## 40.00 0.15

GammaLike(params)

## [1] 2.906766e-22

GammaLogLike(params)

## [1] -49.58983

然后，我們需要預先分配參數！在這種情況下，我們分配gamma（shape = 0.01，scale = 100）和gamma（shape = 0.1，scale = 10）的分布（均值為1且方差略低）：

1.  
   #############
2.  
   # 函數返回參數空間中任意點的先驗概率密度
3.  
   GammaPriorFunction <- function(params){
4.  
   prior <- c(shape=NA,scale=NA)
5.  
   ],3,100)
6.  
   # prior['scale'] <- dunif(params['
7.  
    
8.  
   GammaLogPriorFunction <- function(params){
9.  
   prior <- c(shape=NA,scale=NA)
10.  
    '],shape=0.001,scale=1000,log=T)
11.  
    # prior['shape'] <- dunif(params['shape'],3,100)
12.  
    # prior['scale'] <- dunif(params['
13.  
     
14.  
    curve(dgamma(x,shape=0.01,scale=1000),3,100)

 

params

1.  
   ## shape scale
2.  
   ## 40.00 0.15

GammaPrior(params)

## [1] 1.104038e-06

1.  
   prior2D <- matrix(nrow=length(shapevec),ncol=length(scalevec)) # 初始化存儲變量
2.  
    
3.  
   newparams <- c(shape=50,scale=0
4.  
   for(j in 1:length(scalevec)){
5.  
   newparams['scale'] <- sca
6.  
    
7.  
    
8.  
    
9.  
   ############
10.  
    # 可視化似然面
11.  
    ############
12.  
     
13.  
    image(x=shapevec,y=scalevec,z=prior2D,zlim

 

contour(x=shapevec,y=scalevec,z=prior2D,levels=c(-30,-40,-80,-500),add=T)

我們假設形狀和比例 在先驗中是 獨立的（聯合先驗的乘法概率）。然而，並沒有對后驗參數相關性提出相同的假設，因為概率可以反映在后驗分布中。

然后，我們需要一個函數，該函數可以計算參數空間中任何給定跳轉的后驗概率比率。因為我們正在處理 后驗概率的 比率，所以 我們不需要計算歸一化常數。

無需歸一化常數，我們只需要計算加權似然比（即先驗加權的似然比）

1.  
   ############
2.  
   # 函數用於計算參數空間中任意兩點之間的后驗密度比
3.  
    
4.  
   PosteriorRatio <- function(oldguess,newguess
5.  
   oldLik <- max(1e-90,GammaLikelihoodFunction(oldguess)) # 計算舊猜測的可能性和先驗密度
6.  
    
7.  
   newLik <- GammaLikelihoodFunction(newguess) # 在新的猜測值下計算可能性和先驗密度
8.  
    
9.  
   return((newLik*newPrior)/(oldLik*oldPrior)) # 計算加權似然比
10.  
    }
11.  
     
12.  
    PosteriorRatio2 <- function(oldguess,newguess){
13.  
    oldLogLik <- GammaLogLikelihoodFunction(oldguess) # 計算舊猜測的可能性和先驗密度
14.  
     
15.  
    newLogLik <- GammaLogLikelihoodFunction(newguess) # 在新的猜測值下計算可能性和先驗密度
16.  
     
17.  
    return(exp((newLogLik+newLogPrior)-(oldLogLik+oldLogPrior))) # 計算加權似然比
18.  
    }
19.  
     
20.  
     

## [1] 0.01436301

PosteriorRatio2(oldguess,newguess)

## [1] 0.01436301

然后，我們需要一個函數進行新的推測或在參數空間中跳轉：

1.  
   ############
2.  
   # 為參數空間中的跳轉定義：使用正態分布函數進行新的推測
3.  
   newGuess <- function(oldguess)
4.  
    
5.  
   jump <- c(shape=rnorm(1,mean=0,sd=sdshapejump),scale=rnorm(1,0,sdscalejump))
6.  
   newguess <- abs(oldguess + ju
7.  
    
8.  
   }
9.  
   # 在原始推測附近設置新的推測
10.  
     
11.  
    newGuess(oldguess=params)

1.  
   ## shape scale
2.  
   ## 35.7132110 0.1576337

newGuess(oldguess=params)

1.  
   ## shape scale
2.  
   ## 45.1202345 0.2094243

newGuess(oldguess=params)

1.  
   ## shape scale
2.  
   ## 42.87840436 0.08152061

現在，我們准備實現Metropolis-Hastings MCMC算法：

我們需要一個初始點：

1.  
   ##########
2.  
   # 在參數spacer中設置起點
3.  
    
4.  
   startingvals <- c(shape=75,scale=0.28) # 算法的起點

測試函數

1.  
   ###########
2.  
   # 嘗試我們的新函數
3.  
    
4.  
   newguess <- newGuess(startingvals) # 在參數空間中跳躍
5.  
   newguess

1.  
   ## shape scale
2.  
   ## 73.9663949 0.3149796

PosteriorRatio2(startingvals,newguess)   # 后驗比例差異

## [1] 2.922783e-57

現在讓我們看一下Metropolis-Hastings：

1.  
   ###############
2.  
   #可視化Metropolis-Hastings
3.  
    
4.  
   chain.length <- 11
5.  
   gth,ncol=2)
6.  
   colnames(guesses) <- names(startingvals)
7.  
   guesses[1,] <- startingvals
8.  
   counter <- 2
9.  
    
10.  
    post.rat <- PosteriorRatio2(oldguess,newguess)
11.  
    prob.accept <- min(1,post
12.  
     
13.  
    oldguess <- newguess
14.  
    guesses[coun
15.  
     
16.  
     
17.  
     
18.  
    #可視化
19.  
     
20.  
    contour(x=shapevec,y=scal

 

我們運行更長的時間 

1.  
   ##########
2.  
   # 獲取更多MCMC示例
3.  
    
4.  
   chain.length <- 100
5.  
   oldgu
6.  
    
7.  
   counter <- 2
8.  
   while(counter <= chain.length){
9.  
   newguess <- newGuess(oldguess)
10.  
    post.rat <- Posterio
11.  
     
12.  
    rand <- runif(1)
13.  
    if(rand<=prob.accept){
14.  
    ewguess
15.  
    counter=counte
16.  
     
17.  
    #可視化
18.  
     
19.  
    image(x=shapevec,y=scalevec,z=su
20.  
    urface2D,levels=c(-30,-40,-80,-5
21.  
    lines(guesses,col="red")

更長的時間

1.  
   ############
2.  
   #更長的時間
3.  
    
4.  
    
5.  
    
6.  
   chain.length <- 1000
7.  
   oldguess <- startingvals
8.  
   chain.length,ncol=2)
9.  
   colnames(guesses) <- names(startingvals)
10.  
    guesses[1,] <- startingva
11.  
    ess)
12.  
    post.rat <- PosteriorRatio2(oldguess,newguess)
13.  
    prob.accept <- min(1,post.rat)
14.  
    rand <- runif(1)
15.  
     
16.  
    guesse
17.  
     
18.  
     
19.  
    #可視化
20.  
     
21.  
    image(x=shapevec,y=scalevec,
22.  
    rface2D,levels=c(-30,-40,-80,-500),a

看起來更好！搜索算法可以很好地找到參數空間的高似然部分！

現在，讓我們看一下“ shape”參數的鏈

1.  
   #############
2.  
   # 評估MCMC樣本的“軌跡圖” ...
3.  
    
4.  
   ##### Shape 參數
5.  
    
6.  
   plot(1:chain.length,guesses[,'sha

 

對於比例參數 

1.  
   ###### 比例參數 
2.  
    
3.  
   plot(1:chain.length,guesses[,'scale'],type="l

 

我們可以說這些鏈已經收斂於形狀參數的后驗分布嗎？

首先，鏈的起點“記住”起始值，因此不是固定分布。我們需要刪除鏈的第一部分。

1.  
   ############
2.  
   # 刪除預燒期（允許MCMC有一段時間達到后驗概率）
3.  
    
4.  
   burn.in <- 100
5.  
   MCMCsamples <- guesses[-c(1:burn.in),]
6.  
    

 

 

但這看起來還不是很好。讓我們運行更長的時間，看看是否得到的東西看起來更像是隨機數生成器（白噪聲）

1.  
   ##########
2.  
   # 再試一次-運行更長的時間
3.  
    
4.  
   chain.length <- 20000
5.  
   oldguess <- startingv
6.  
   o2(oldguess,newguess)
7.  
   prob.accept <- mi
8.  
    
9.  
    

 

讓我們首先刪除前5000個樣本作為預燒期

1.  
   #############
2.  
   # 使用更長的“預燒”
3.  
    
4.  
   burn.in <- 5000
5.  
   MCMCsamples <- guesses[-c(1:bur

現在，讓我們再次看一下鏈條

 

 

在評估這些跡線圖時，我們希望看到看起來像白噪聲的“平穩分布”。該軌跡圖看起來可能具有一些自相關。解決此問題的一種方法是稀疏MCMC樣本：

1.  
   ##########
2.  
   # “稀疏” MCMC樣本
3.  
    
4.  
   thinnedMCMC <- MCMCsamples[seq(1,chain.length,by=5),]

 

 

現在我們可以檢查我們的后驗分布！

1.  
   # 可視化后驗分布
2.  
    
3.  
   plot(density(thinnedMCMC[,'scale'])

 

 

我們可以像以前一樣可視化。

1.  
   #########
2.  
   # 更多后驗概率檢察
3.  
    
4.  
   par(mfrow=c(3,2))
5.  
   plot(thinnedMCMC,col=1:10000)
6.  
   plot(thinnedMCMC,type="l")

 

可以修改Metropolis-Hastings MCMC方法來擬合任意模型的任意數量的自由參數。但是，MH算法本身不一定是最有效和靈活的。在實驗中，我們使用吉布斯采樣，大多采用建模語言 BUGS 。

注意：BUGS實現（例如JAGS）實際上傾向於結合使用MH和Gibbs采樣，MH和Gibbs采樣器並不是唯一的MCMC例程。例如，“ stan”使用MH采樣的一種改進形式，稱為“ Hamiltonian Monte Carlo”。

吉布斯Gibbs采樣器

Gibbs采樣器非常簡單有效。基本上，該算法從完整的條件 概率分布（即， 在模型中所有其他參數的已知值作為條件的條件下，對任意參數i的后驗分布）中進行 連續采樣 。

在很多情況下，我們不能直接制定出我們的模型后驗分布，但我們 可以 分析出條件后驗分布。盡管如此，即使它在分析上不易處理，我們也可以使用單變量MH程序作為最后方法。

問：為什么Gibbs采樣器通常比純MH采樣器效率更高？

二元正態例子

MCMC采樣器只是隨機數生成器的一種。我們可以使用Gibbs采樣器來開發自己的隨機數生成器，以實現相當簡單的已知分布。在此示例中，我們使用Gibbs采樣器從標准雙變量正態概率分布生成隨機數。注意，吉布斯采樣器在許多方面都比MH算法更簡單明了。

1.  
   #############
2.  
   #Gibbs采樣器的簡單示例
3.  
   #############
4.  
    
5.  
   ########
6.  
   # 首先，回顧一下我們簡單的雙變量正態采樣器
7.  
   rbvn<-function (n, rho){ #f函數用於從雙變量標准正態分布中提取任意數量的獨立樣本。
8.  
   x <- rnorm(n, 0, 1)
9.  
   y <- rnorm(n, rho * x, sqrt(1 - rho^2))
10.  
     

 

1.  
   #############
2.  
   # 現在構造一個吉布斯采樣器
3.  
   gibbs<-function (n, rho){ # 雙變量隨機數生成器的gibbs采樣器實現
4.  
   mat <- matrix(ncol = 2, nrow = n) # 用於存儲隨機樣本的矩陣
5.  
    
6.  
   mat[1, ] <- c(x, y) # initialize the markov chain
7.  
   for (i in 2:n) {
8.  
   x <- rnorm(1, rho * y, sqrt(1 - rho^2)) # 以y為條件的x中的樣本
9.  
   y <- rnorm(1, rho * x, sqrt(1 - rho^2)) # 以x為條件的y中的樣本
10.  
    mat[i, ] <- c(x, y)
11.  
     

然后，我們可以使用Gibbs采樣器從該已知分布中獲取隨機樣本…

1.  
   ##########
2.  
   # 測試吉布斯采樣器
3.  
    
4.  
    
5.  
    
6.  
   plot(ts(bvn[,2]))
7.  
   hist(bvn[,1],40)
8.  
   hist(bvn[,2],40)

 

在這里，馬爾可夫鏈的樣本中有很多明顯的自相關。Gibbs采樣器經常有此問題。

示例

BUGS語言

最后，讓我們為我們最喜​​歡的粘瘤病示例創建一個Gibbs采樣器，為此，我們將使用BUGS語言（在JAGS中實現）來幫助我們！

JAGS，全稱是Just another Gibbs sampler，是基於BUGS語言開發的利用MCMC來進行貝葉斯建模的軟件包。它沒有提供建模所用的GUI以及MCMC抽樣的后處理，這些要在其它的程序軟件上來處理，比如說利用R包（rjags）來調用JAGS並后處理MCMC的輸出。JAGS相對於WinBUGS/OpenBUGS的主要優點在於平台的獨立性，可以應用於各種操作系統，而WinBUGS/OpenBUGS只能應用於windows系統；JAGS也可以在64-bit平台上以64-bit應用來進行編譯。

BUGS語言看起來與R類似，但是有幾個主要區別：

• 首先，BUGS是一種編譯語言，因此代碼中的操作順序並不重要
• BUGS不是矢量化的-您需要使用FOR循環
• 在BUGS中，幾個概率分布的參數差異很大。值得注意的是，正態分布通過均值和精度（1 / Variance ）進行參數化。

這是用BUGS語言實現的粘液病示例：

 

1.  
    
2.  
   model {
3.  
    
4.  
   #############
5.  
   # 似然
6.  
   ############
7.  
   for(obs in 1:n.observations){
8.  
   titer[obs] ~ dgamma(shape,rate
9.  
    
10.  
    #############
11.  
    # 先驗
12.  
    ############
13.  
     
14.  
     
15.  
     
16.  
    rate <- 1/scale # 將BUGS的scale參數轉換為“ rate”
17.  
    }

我們可以使用R中的“ cat”函數將此模型寫到您的工作目錄中的文本文件中：

1.  
   ###########
2.  
   # BUGS建模語言中的粘液瘤示例
3.  
    
4.  
   ##########
5.  
   # 將BUGS模型寫入文件
6.  
    
7.  
   cat("
8.  
   model {
9.  
    
10.  
    #############
11.  
    # 似然
12.  
    ############
13.  
    for(obs in 1:n.observations){
14.  
    titer[obs] ~ dgamma(shape,rate)
15.  
     
16.  
     
17.  
     
18.  
    #############
19.  
    # 先驗
20.  
    ############
21.  
    shape ~ dgamma(0.001,0.001
22.  
     
23.  
    file="BUGSmodel.txt")

現在我們已經將BUGS模型打包為文本文件，我們將數據捆綁到一個列表對象中，該列表對象包含BUGS代碼中引用的所有相關數據：

1.  
   ############
2.  
   # 將數據封裝到單個“列表”對象中
3.  
    
4.  
   myx.data <- list(
5.  
    
6.  
   n.observations = length(Myx$titer

1.  
   ## $titer
2.  
   ## [1] 5.207 5.734 6.613 5.997 6.612 6.810 5.930 6.501 7.182 7.292 7.819
3.  
   ## [12] 7.489 6.918 6.808 6.235 6.916 4.196 7.682 8.189 7.707 7.597 7.112
4.  
   ## [23] 7.354 7.158 7.466 7.927 8.499
5.  
   ##
6.  
   ## $n.observations
7.  
   ## [1] 27

然后，我們需要為所有參數定義初始值。將其定義為一個函數很方便，因此可以使用不同的起始值來初始化每個MCMC鏈。 

1.  
   ###########
2.  
   # 用於為所有自由參數生成隨機初始值的函數
3.  
    
4.  
    
5.  
    
6.  
    
7.  
   init.vals.for.bugs()

1.  
   ## $shape
2.  
   ## [1] 51.80414
3.  
   ##
4.  
   ## $scale
5.  
   ## [1] 0.2202549

init.vals.for.bugs()

1.  
   ## $shape
2.  
   ## [1] 89.85618
3.  
   ##
4.  
   ## $scale
5.  
   ## [1] 0.2678733

init.vals.for.bugs()

1.  
   ## $shape
2.  
   ## [1] 69.22457
3.  
   ##
4.  
   ## $scale
5.  
   ## [1] 0.1728908

現在我們可以調用JAGS！

1.  
   ###########
2.  
   # 運行 JAGS
3.  
   ##########

## Loading required package: rjags

1.  
   ## The following object is masked from 'package:coda':
2.  
   ##
3.  
   ## traceplot

1.  
    
2.  
    
3.  
   jags.fit <- run.jags(model="BUGSmodel.txt",
4.  
    

1.  
   ## Compiling rjags model...
2.  
   ## Calling the simulation using the rjags method...
3.  
   ## Adapting the model for 100 iterations...
4.  
   ## Running the model for 5000 iterations...
5.  
   ## Simulation complete
6.  
   ## Finished running the simulation

1.  
   jags.fit) # 轉換為“ MCMC”對象（CODA包）
2.  
    
3.  
    

1.  
   ##
2.  
   ## Iterations = 101:5100
3.  
   ## Thinning interval = 1
4.  
   ## Number of chains = 3
5.  
   ## Sample size per chain = 5000
6.  
   ##
7.  
   ## 1. Empirical mean and standard deviation for each variable,
8.  
   ## plus standard error of the mean:
9.  
   ##
10.  
    ## Mean SD Naive SE Time-series SE
11.  
    ## shape 51.6013 14.36711 0.1173070 2.216657
12.  
    ## scale 0.1454 0.04405 0.0003597 0.006722
13.  
    ##
14.  
    ## 2. Quantiles for each variable:
15.  
    ##
16.  
    ## 2.5% 25% 50% 75% 97.5%
17.  
    ## shape 28.76333 40.9574 50.1722 60.2463 82.7532
18.  
    ## scale 0.08346 0.1148 0.1378 0.1687 0.2389

plot(jagsfit.mcmc)

評估收斂

第一步是視覺檢查-我們尋找以下內容來評估收斂性：

• 當視為“軌跡圖”時，每個參數的鏈應看起來像白噪聲（模糊的毛毛蟲）或類似的噪聲
• 多個具有不同起始條件的鏈條看起來應該相同

我們可能在這里可以做得更好的一種方法是使鏈條運行更長的時間，並丟棄初始樣本我們還可以。

通過細化鏈來嘗試減少自相關，我們每20個樣本中僅保留1個。

1.  
   ################
2.  
   #運行鏈更長時間
3.  
    
4.  
   jags.fit <-
5.  
   sample = 10000,n.chains = 3,adapt = 1000,burnin = 1000,
6.  
   summarise = FALSE,thin=20,method = "parallel" )

1.  
   ## Calling 3 simulations using the parallel method...
2.  
   ## Following the progress of chain 1 (the program will wait for all
3.  
   ## chains to finish before continuing):
4.  
   ## Welcome to JAGS 4.2.0 on Wed Oct 23 11:33:35 2019
5.  
   ## JAGS is free software and comes with ABSOLUTELY NO WARRANTY
6.  
   ## Loading module: basemod: ok
7.  
   ## Loading module: bugs: ok
8.  
   ## . . Reading data file data.txt
9.  
   ## . Compiling model graph
10.  
    ## Resolving undeclared variables
11.  
    ## Allocating nodes
12.  
    ## Graph information:
13.  
    ## Observed stochastic nodes: 27
14.  
    ## Unobserved stochastic nodes: 2
15.  
    ## Total graph size: 37
16.  
    ## . Reading parameter file inits1.txt
17.  
    ## . Initializing model
18.  
    ## . Adapting 1000
19.  
    ## -------------------------------------------------| 1000
20.  
    ## ++++++++++++++++++++++++++++++++++++++++++++++++++ 100%
21.  
    ## Adaptation successful
22.  
    ## . Updating 1000
23.  
    ## -------------------------------------------------| 1000
24.  
    ## ************************************************** 100%
25.  
    ## . . . Updating 200000
26.  
    ## -------------------------------------------------| 200000
27.  
    ## ************************************************** 100%
28.  
    ## . . . . Updating 0
29.  
    ## . Deleting model
30.  
    ## .
31.  
    ## All chains have finished
32.  
    ## Simulation complete. Reading coda files...
33.  
    ## Coda files loaded successfully
34.  
    ## Finished running the simulation

1.  
   jagsfit.mcmc <- as.mcmc.list
2.  
   # 轉換為“ MCMC”對象（CODA包）
3.  
    

1.  
   ##
2.  
   ## Iterations = 2001:201981
3.  
   ## Thinning interval = 20
4.  
   ## Number of chains = 3
5.  
   ## Sample size per chain = 10000
6.  
   ##
7.  
   ## 1. Empirical mean and standard deviation for each variable,
8.  
   ## plus standard error of the mean:
9.  
   ##
10.  
    ## Mean SD Naive SE Time-series SE
11.  
    ## shape 47.1460 12.9801 0.0749404 0.292218
12.  
    ## scale 0.1591 0.0482 0.0002783 0.001075
13.  
    ##
14.  
    ## 2. Quantiles for each variable:
15.  
    ##
16.  
    ## 2.5% 25% 50% 75% 97.5%
17.  
    ## shape 25.14757 37.9988 45.9142 55.1436 75.5850
18.  
    ## scale 0.09147 0.1256 0.1509 0.1825 0.2773

plot(jagsfit.mcmc)

從視覺上看，這看起來更好。現在我們可以使用更多的定量收斂指標。

Gelman-Rubin診斷

一種簡單而直觀的收斂診斷程序是 Gelman-Rubin診斷程序，該診斷程序基於簡單的蒙特卡洛誤差來評估鏈之間是否比應有的鏈距更大：

1.  
   ##############
2.  
   # 運行收斂診斷
3.  
    
4.  
   gelman(jagsfit.mcmc)

1.  
   ## Potential scale reduction factors:
2.  
   ##
3.  
   ## Point est. Upper C.I.
4.  
   ## shape 1 1.00
5.  
   ## scale 1 1.01
6.  
   ##
7.  
   ## Multivariate psrf
8.  
   ##
9.  
   ## 1

通常，1.1或更高的值被認為收斂不佳。為模型中的所有可用參數計算GR診斷。如果測試失敗，則應嘗試運行更長的鏈！

所以這個模型看起來不錯！

---

最受歡迎的見解

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模型實現