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

最近我們被客戶要求撰寫關(guān)于BUGS/JAGS貝葉斯分析的研究報告，包括一些圖形和統(tǒng)計輸出。

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

相關(guān)視頻：馬爾可夫鏈原理可視化解釋與R語言區(qū)制轉(zhuǎn)換Markov regime switching實例

馬爾可夫鏈原理可視化解釋與R語言區(qū)制轉(zhuǎn)換Markov?regime?switching實例

相關(guān)視頻

馬爾可夫鏈蒙特卡羅方法MCMC原理與R語言實現(xiàn)

，時長08:47

馬爾科夫鏈蒙特卡洛方法

MCMC的關(guān)鍵如下：

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

跳躍概率可以表征為：

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

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

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

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

令人驚訝的是，MCMC的核心并不難于描述或?qū)嵤?。讓我們看一個簡單的MCMC算法。

Metropolis-Hastings算法

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

MH算法可以表示為：

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

請注意，從本質(zhì)上講，這與“ Metropolis”模擬退火算法相同，后驗概率代替了概率，并且?k?參數(shù)設(shè)置為1。

二元正態(tài)例子

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

對于這個簡單的示例，我們不需要MCMC采樣器。一種實現(xiàn)方法是使用以下代碼，該代碼從具有相關(guān)參數(shù)ρ的雙變量標準正態(tài)分布中繪制并可視化任意數(shù)量的獨立樣本。

#################＃MCMC采樣的簡單示例########################### ＃首先，讓我們構(gòu)建一個從雙變量標準正態(tài)分布生成隨機數(shù)的函數(shù)rbvn<-function (n, rho) ? #用于從二元標準正態(tài)分布中提取任意數(shù)量的獨立樣本。{ ? ? ? ?x <- rnorm(n, 0, 1) ? ? ? ?y <- rnorm(n, rho * x, sqrt(1 - rho^2)) ? ? ? ?cbind(x, y)}########## 現(xiàn)在，從該分布圖中繪制隨機抽樣 bvn<-rbvn(10000,0.98)par(mfrow=c(3,2))plot(bvn,col=1:10000

################ ＃Metropolis-Hastings雙變量正態(tài)采樣器的實現(xiàn)...library(mvtnorm) ? ?# 加載一個包，該包使我們能夠計算mv正態(tài)分布的概率密度metropoli<- function (n, rho=0.98){ ? ?# 雙變量隨機數(shù)生成器的MCMC采樣器實現(xiàn) ? ?mat <- matrix(ncol = 2, nrow = n) ? # 用于存儲隨機樣本的矩陣 ? ?x <- 0 ? # 所有參數(shù)的初始值 ? ?prev <- dmvnorm(c(x,y),mean=c(0,0),sig # 起始位置分布的概率密度 ? ?mat[1, ] <- c(x, y) ? ? ? ?# 初始化馬爾可夫鏈 ? ? ?newx <- rnorm(1,x,0.5) ? ? # 進行跳轉(zhuǎn) ? ? ?newprob <- dmvnorm(c(newx,newy),sigma = ?# 評估跳轉(zhuǎn) ? ? ?ratio <- newprob/prev ? # 計算舊位置（跳出）和建議位置（跳到）的概率之比。 ? ? ?prob.accept <- min(1,ratio) ? ? # 決定接受新跳躍的概率！ ? ? ?if(rand<=prob.accept){ ? ? ? ?x=newx;y=newy ? ?# 將x和y設(shè)置為新位置 ? ? ? ?mat[counter,] <- c(x,y) ? ?# 將其存儲在存儲陣列中 ? ? ? ?prev <- newprob ? ?# 準備下一次迭代

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

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

點擊標題查閱往期內(nèi)容

R語言用WinBUGS 軟件對學(xué)術(shù)能力測驗建立層次（分層）貝葉斯模型

左右滑動查看更多

01

02

03

04

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

MCMC對粘液瘤病進行調(diào)查

#############黏液病示例的MCMC實現(xiàn)############

subset(MyxDat,grade==1

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

選擇使用Gamma分布。這是經(jīng)驗分布：

############ 第100次可視化粘液病數(shù)據(jù)hist(Myx$titer,freq=FALSE)

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

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

二維（對數(shù)）似然面：

############### 定義二維參數(shù)空間##############shapevec <- seq(3,100,by=0.1) ? scalevec <- seq(0.01,0.5,by=0.001)############### ＃定義參數(shù)空間內(nèi)此網(wǎng)格上的似然面##############GammaLogLikelihoodFunction <- function(par}surface2D <- matrix(nrow=length(shapevec),ncol=length(scalevec)) ? #初始化存儲變量newparams <- c(sha ? ?surface2D[i,j] <- GammaLogLikelihoodFunction(newparams) ############# 可視化似然面############contour(x=shapevec,y=scalevec,z=surface2D,levels=c(-30,-40,-80,-500),add=T)

這是MH算法的實現(xiàn)，用于找到后驗分布！

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

#############編寫非對數(shù)轉(zhuǎn)換的似然函數(shù)GammaLike- function(params){ ?prod(dgamma(Myx$titer,shape=params['shape']params <- c(shape=40,

## shape scale ## 40.00 ?0.15

GammaLike(params)

## [1] 2.906766e-22

GammaLogLike(params)

## [1] -49.58983

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

############## 函數(shù)返回參數(shù)空間中任意點的先驗概率密度GammaPriorFunction <- function(params){ ?prior <- c(shape=NA,scale=NA)],3,100) ? ? ? ? ?# prior['scale'] <- dunif(params['GammaLogPriorFunction <- function(params){ ?prior <- c(shape=NA,scale=NA)'],shape=0.001,scale=1000,log=T) ?# prior['shape'] <- dunif(params['shape'],3,100) ? ? ? ? ?# prior['scale'] <- dunif(params['curve(dgamma(x,shape=0.01,scale=1000),3,100)

params

## shape scale ## 40.00 ?0.15

GammaPrior(params)

## [1] 1.104038e-06

prior2D <- matrix(nrow=length(shapevec),ncol=length(scalevec)) ? # 初始化存儲變量newparams <- c(shape=50,scale=0 ?for(j in 1:length(scalevec)){ ? ?newparams['scale'] <- sca ############# 可視化似然面############image(x=shapevec,y=scalevec,z=prior2D,zlim

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

我們假設(shè)形狀和比例?在先驗中是?獨立的（聯(lián)合先驗的乘法概率）。然而，并沒有對后驗參數(shù)相關(guān)性提出相同的假設(shè)，因為概率可以反映在后驗分布中。

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

無需歸一化常數(shù)，我們只需要計算加權(quán)似然比（即先驗加權(quán)的似然比）

############# 函數(shù)用于計算參數(shù)空間中任意兩點之間的后驗密度比PosteriorRatio <- function(oldguess,newguess ?oldLik <- max(1e-90,GammaLikelihoodFunction(oldguess)) ? # 計算舊猜測的可能性和先驗密度 ?newLik <- GammaLikelihoodFunction(newguess) ? ? ? ? ? ? # 在新的猜測值下計算可能性和先驗密度 ?return((newLik*newPrior)/(oldLik*oldPrior)) ? ? ? ? ?# 計算加權(quán)似然比}PosteriorRatio2 <- function(oldguess,newguess){ ?oldLogLik <- GammaLogLikelihoodFunction(oldguess) ? # 計算舊猜測的可能性和先驗密度 ?newLogLik <- GammaLogLikelihoodFunction(newguess) ? ? ? ? ? ? # 在新的猜測值下計算可能性和先驗密度 ?return(exp((newLogLik+newLogPrior)-(oldLogLik+oldLogPrior))) ? ? ? ? ?# 計算加權(quán)似然比}

## [1] 0.01436301

PosteriorRatio2(oldguess,newguess)

## [1] 0.01436301

然后，我們需要一個函數(shù)進行新的推測或在參數(shù)空間中跳轉(zhuǎn)：

############# 為參數(shù)空間中的跳轉(zhuǎn)定義：使用正態(tài)分布函數(shù)進行新的推測newGuess <- function(oldguess) ?jump <- c(shape=rnorm(1,mean=0,sd=sdshapejump),scale=rnorm(1,0,sdscalejump)) ?newguess <- abs(oldguess + ju} ?# 在原始推測附近設(shè)置新的推測newGuess(oldguess=params) ?

## ? ? ?shape ? ? ?scale ## 35.7132110 ?0.1576337

newGuess(oldguess=params)

## ? ? ?shape ? ? ?scale ## 45.1202345 ?0.2094243

newGuess(oldguess=params)

## ? ? ? shape ? ? ? scale ## 42.87840436 ?0.08152061

現(xiàn)在，我們準備實現(xiàn)Metropolis-Hastings MCMC算法：

我們需要一個初始點：

########### 在參數(shù)spacer中設(shè)置起點startingvals <- c(shape=75,scale=0.28) ? ?# 算法的起點

測試函數(shù)

############ 嘗試我們的新函數(shù)newguess <- newGuess(startingvals) ? ?# 在參數(shù)空間中跳躍newguess

## ? ? ?shape ? ? ?scale ## 73.9663949 ?0.3149796

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

## [1] 2.922783e-57

現(xiàn)在讓我們看一下Metropolis-Hastings：

################可視化Metropolis-Hastingschain.length <- 11gth,ncol=2)colnames(guesses) <- names(startingvals)guesses[1,] <- startingvalscounter <- 2 ?post.rat <- PosteriorRatio2(oldguess,newguess) ?prob.accept <- min(1,post ? ?oldguess <- newguess ? ?guesses[coun#可視化contour(x=shapevec,y=scal

我們運行更長的時間?

########### 獲取更多MCMC示例chain.length <- 100oldgucounter <- 2while(counter <= chain.length){ ?newguess <- newGuess(oldguess) ?post.rat <- Posterio ?rand <- runif(1) ?if(rand<=prob.accept){ewguess ? ?counter=counte#可視化image(x=shapevec,y=scalevec,z=suurface2D,levels=c(-30,-40,-80,-5lines(guesses,col="red")

更長的時間

#############更長的時間chain.length <- 1000oldguess <- startingvalschain.length,ncol=2)colnames(guesses) <- names(startingvals)guesses[1,] <- startingvaess) ?post.rat <- PosteriorRatio2(oldguess,newguess) ?prob.accept <- min(1,post.rat) ?rand <- runif(1) ? ?guesse#可視化image(x=shapevec,y=scalevec,rface2D,levels=c(-30,-40,-80,-500),a

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

現(xiàn)在，讓我們看一下“ shape”參數(shù)的鏈

############## 評估MCMC樣本的“軌跡圖” ...##### Shape 參數(shù)plot(1:chain.length,guesses[,'sha

對于比例參數(shù)?

###### 比例參數(shù)?plot(1:chain.length,guesses[,'scale'],type="l

我們可以說這些鏈已經(jīng)收斂于形狀參數(shù)的后驗分布嗎？

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

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

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

########### 再試一次-運行更長的時間chain.length <- 20000oldguess <- startingvo2(oldguess,newguess) ?prob.accept <- mi

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

############## 使用更長的“預(yù)燒”burn.in <- 5000MCMCsamples <- guesses[-c(1:bur

現(xiàn)在，讓我們再次看一下鏈條

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

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

現(xiàn)在我們可以檢查我們的后驗分布！

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

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

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

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

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

吉布斯Gibbs采樣器

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

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

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

二元正態(tài)例子

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

##############Gibbs采樣器的簡單示例###################### 首先，回顧一下我們簡單的雙變量正態(tài)采樣器rbvn<-function (n, rho){ ?#f函數(shù)用于從雙變量標準正態(tài)分布中提取任意數(shù)量的獨立樣本。 ? ? ? ?x <- rnorm(n, 0, 1) ? ? ? ?y <- rnorm(n, rho * x, sqrt(1 - rho^2))

############## 現(xiàn)在構(gòu)造一個吉布斯采樣器gibbs<-function (n, rho){ ? ?# 雙變量隨機數(shù)生成器的gibbs采樣器實現(xiàn) ? ?mat <- matrix(ncol = 2, nrow = n) ? # 用于存儲隨機樣本的矩陣 ? ?mat[1, ] <- c(x, y) ? ? ? ?# initialize the markov chain ? ?for (i in 2:n) { ? ? ? ? ? ?x <- rnorm(1, rho * y, sqrt(1 - rho^2)) ? ? ? ?# 以y為條件的x中的樣本 ? ? ? ? ? ?y <- rnorm(1, rho * x, sqrt(1 - rho^2)) ? ? ? ?# 以x為條件的y中的樣本 ? ? ? ? ? ?mat[i, ] <- c(x, y)

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

########### 測試吉布斯采樣器plot(ts(bvn[,2]))hist(bvn[,1],40)hist(bvn[,2],40)

在這里，馬爾可夫鏈的樣本中有很多明顯的自相關(guān)。Gibbs采樣器經(jīng)常有此問題。

示例

BUGS語言

最后，讓我們?yōu)槲覀冏钕矚g的粘瘤病示例創(chuàng)建一個Gibbs采樣器，為此，我們將使用BUGS語言（在JAGS中實現(xiàn)）來幫助我們！

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

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

• 首先，BUGS是一種編譯語言，因此代碼中的操作順序并不重要

• BUGS不是矢量化的-您需要使用FOR循環(huán)

• 在BUGS中，幾個概率分布的參數(shù)差異很大。值得注意的是，正態(tài)分布通過均值和精度（1 / Variance ）進行參數(shù)化。

這是用BUGS語言實現(xiàn)的粘液病示例：

model { ?############# ?# 似然 ?############ ?for(obs in 1:n.observations){ ? ?titer[obs] ~ dgamma(shape,rate ?############# ?# 先驗 ?############ ?rate <- 1/scale ? # 將BUGS的scale參數(shù)轉(zhuǎn)換為“ rate”}

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

############ BUGS建模語言中的粘液瘤示例########### 將BUGS模型寫入文件cat(" ?model { ? ?############# ? ?# 似然 ? ?############ ? ?for(obs in 1:n.observations){ ? ? ?titer[obs] ~ dgamma(shape,rate) ? ?############# ? ?# 先驗 ? ?############ ? ?shape ~ dgamma(0.001,0.001file="BUGSmodel.txt")

現(xiàn)在我們已經(jīng)將BUGS模型打包為文本文件，我們將數(shù)據(jù)捆綁到一個列表對象中，該列表對象包含BUGS代碼中引用的所有相關(guān)數(shù)據(jù)：

############# 將數(shù)據(jù)封裝到單個“列表”對象中myx.data <- list( ?n.observations = length(Myx$titer

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

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

############ 用于為所有自由參數(shù)生成隨機初始值的函數(shù)init.vals.for.bugs()

## $shape## [1] 51.80414## ## $scale## [1] 0.2202549

init.vals.for.bugs()

## $shape## [1] 89.85618## ## $scale## [1] 0.2678733

init.vals.for.bugs()

## $shape## [1] 69.22457## ## $scale## [1] 0.1728908

現(xiàn)在我們可以調(diào)用JAGS！

############ 運行 JAGS ##########

## Loading required package: rjags

## The following object is masked from 'package:coda':## ## ? ? traceplot

jags.fit <- run.jags(model="BUGSmodel.txt",

## Compiling rjags model...## Calling the simulation using the rjags method...## Adapting the model for 100 iterations...## Running the model for 5000 iterations...## Simulation complete## Finished running the simulation

jags.fit) ? # 轉(zhuǎn)換為“ MCMC”對象（CODA包）

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

plot(jagsfit.mcmc)

評估收斂

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

• 當(dāng)視為“軌跡圖”時，每個參數(shù)的鏈應(yīng)看起來像白噪聲（模糊的毛毛蟲）或類似的噪聲

• 多個具有不同起始條件的鏈條看起來應(yīng)該相同

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

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

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

## Calling 3 simulations using the parallel method...## Following the progress of chain 1 (the program will wait for all## chains to finish before continuing):## Welcome to JAGS 4.2.0 on Wed Oct 23 11:33:35 2019## JAGS is free software and comes with ABSOLUTELY NO WARRANTY## Loading module: basemod: ok## Loading module: bugs: ok## . . Reading data file data.txt## . Compiling model graph## ? ?Resolving undeclared variables## ? ?Allocating nodes## Graph information:## ? ?Observed stochastic nodes: 27## ? ?Unobserved stochastic nodes: 2## ? ?Total graph size: 37## . Reading parameter file inits1.txt## . Initializing model## . Adapting 1000## -------------------------------------------------| 1000## ++++++++++++++++++++++++++++++++++++++++++++++++++ 100%## Adaptation successful## . Updating 1000## -------------------------------------------------| 1000## ************************************************** 100%## . . . Updating 200000## -------------------------------------------------| 200000## ************************************************** 100%## . . . . Updating 0## . Deleting model## . ## All chains have finished## Simulation complete. ?Reading coda files...## Coda files loaded successfully## Finished running the simulation

jagsfit.mcmc <- as.mcmc.list ?# 轉(zhuǎn)換為“ MCMC”對象（CODA包）

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

plot(jagsfit.mcmc)

從視覺上看，這看起來更好?，F(xiàn)在我們可以使用更多的定量收斂指標。

Gelman-Rubin診斷

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

############### 運行收斂診斷gelman(jagsfit.mcmc)

## Potential scale reduction factors:## ## ? ? ? Point est. Upper C.I.## shape ? ? ? ? ?1 ? ? ? 1.00## scale ? ? ? ? ?1 ? ? ? 1.01## ## Multivariate psrf## ## 1

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

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

點擊文末?“閱讀原文”

獲取全文完整代碼數(shù)據(jù)資料。

本文選自《R語言BUGS/JAGS貝葉斯分析: 馬爾科夫鏈蒙特卡洛方法（MCMC）采樣》。

點擊標題查閱往期內(nèi)容

R語言MCMC:Metropolis-Hastings采樣用于回歸的貝葉斯估計
Python用MCMC馬爾科夫鏈蒙特卡洛、拒絕抽樣和Metropolis-Hastings采樣算法
R語言貝葉斯METROPOLIS-HASTINGS GIBBS 吉布斯采樣器估計變點指數(shù)分布分析泊松過程車站等待時間
R語言馬爾可夫MCMC中的METROPOLIS HASTINGS，MH算法抽樣（采樣）法可視化實例
python貝葉斯隨機過程：馬爾可夫鏈Markov-Chain，MC和Metropolis-Hastings，MH采樣算法可視化
Python貝葉斯推斷Metropolis-Hastings（M-H）MCMC采樣算法的實現(xiàn)
Metropolis Hastings采樣和貝葉斯泊松回歸Poisson模型
Matlab用BUGS馬爾可夫區(qū)制轉(zhuǎn)換Markov switching隨機波動率模型、序列蒙特卡羅SMC、M H采樣分析時間序列R語言RSTAN MCMC：NUTS采樣算法用LASSO 構(gòu)建貝葉斯線性回歸模型分析職業(yè)聲望數(shù)據(jù)
R語言BUGS序列蒙特卡羅SMC、馬爾可夫轉(zhuǎn)換隨機波動率SV模型、粒子濾波、Metropolis Hasting采樣時間序列分析
R語言Metropolis Hastings采樣和貝葉斯泊松回歸Poisson模型
R語言貝葉斯MCMC：用rstan建立線性回歸模型分析汽車數(shù)據(jù)和可視化診斷
R語言貝葉斯MCMC：GLM邏輯回歸、Rstan線性回歸、Metropolis Hastings與Gibbs采樣算法實例
R語言貝葉斯Poisson泊松-正態(tài)分布模型分析職業(yè)足球比賽進球數(shù)
R語言用Rcpp加速Metropolis-Hastings抽樣估計貝葉斯邏輯回歸模型的參數(shù)
R語言邏輯回歸、Naive Bayes貝葉斯、決策樹、隨機森林算法預(yù)測心臟病
R語言中貝葉斯網(wǎng)絡(luò)（BN）、動態(tài)貝葉斯網(wǎng)絡(luò)、線性模型分析錯頜畸形數(shù)據(jù)
R語言中的block Gibbs吉布斯采樣貝葉斯多元線性回歸
Python貝葉斯回歸分析住房負擔(dān)能力數(shù)據(jù)集
R語言實現(xiàn)貝葉斯分位數(shù)回歸、lasso和自適應(yīng)lasso貝葉斯分位數(shù)回歸分析
Python用PyMC3實現(xiàn)貝葉斯線性回歸模型
R語言用WinBUGS 軟件對學(xué)術(shù)能力測驗建立層次（分層）貝葉斯模型
R語言Gibbs抽樣的貝葉斯簡單線性回歸仿真分析
R語言和STAN,JAGS：用RSTAN,RJAG建立貝葉斯多元線性回歸預(yù)測選舉數(shù)據(jù)
R語言基于copula的貝葉斯分層混合模型的診斷準確性研究
R語言貝葉斯線性回歸和多元線性回歸構(gòu)建工資預(yù)測模型
R語言貝葉斯推斷與MCMC：實現(xiàn)Metropolis-Hastings 采樣算法示例
R語言stan進行基于貝葉斯推斷的回歸模型
R語言中RStan貝葉斯層次模型分析示例
R語言使用Metropolis-Hastings采樣算法自適應(yīng)貝葉斯估計與可視化
R語言隨機搜索變量選擇SSVS估計貝葉斯向量自回歸（BVAR）模型
WinBUGS對多元隨機波動率模型：貝葉斯估計與模型比較
R語言實現(xiàn)MCMC中的Metropolis–Hastings算法與吉布斯采樣
R語言貝葉斯推斷與MCMC：實現(xiàn)Metropolis-Hastings 采樣算法示例
R語言使用Metropolis-Hastings采樣算法自適應(yīng)貝葉斯估計與可視化
視頻：R語言中的Stan概率編程MCMC采樣的貝葉斯模型
R語言MCMC:Metropolis-Hastings采樣用于回歸的貝葉斯估計