R語言BUGS/JAGS貝葉斯分析: 馬爾科夫鏈蒙特卡洛方法（MCMC）采樣
原文鏈接：http://tecdat.cn/?p=17884
原文出處：拓端數(shù)據(jù)部落公眾號
?
?
馬爾科夫鏈蒙特卡洛方法
在許多情況下，我們沒有足夠的計算能力評估空間中所有n維像素的后驗概率?。在這些情況下，我們傾向于利用稱為Markov-Chain Monte Carlo?算法的程序?。此方法使用參數(shù)空間中的隨機跳躍來（最終）確定后驗分布。MCMC的關鍵如下：
跳躍概率的比例與后驗概率的比例成正比。
跳躍概率可以表征為：
概率（跳躍）*概率（接受）
從長遠來看，該鏈將花費大量時間在參數(shù)空間的高概率部分，從而實質上捕獲了后驗分布。有了足夠的跳躍，長期分布將與聯(lián)合后驗概率分布匹配。
MCMC本質上是一種特殊類型的隨機數(shù)生成器，旨在從難以描述（例如，多元，分層）的概率分布中采樣。在許多/大多數(shù)情況下，后驗分布是很難描述的概率分布。
MCMC使您可以從實際上不可能完全定義的概率分布中進行采樣！
令人驚訝的是，MCMC的核心并不難于描述或實施。讓我們看一個簡單的MCMC算法。
Metropolis-Hastings算法
該算法與模擬退火算法非常相似。
MH算法可以表示為：
Prob（acceptB | A）= min（1，Posterior（B）Posterior（A）?Prob（b→a）Prob（a→b））
請注意，從本質上講，這與“ Metropolis”模擬退火算法相同，后驗概率代替了概率，并且?k?參數(shù)設置為1。
?
二元正態(tài)例子
?
請記住，MCMC采樣器只是隨機數(shù)生成器的一種。我們可以使用Metropolis-Hastings采樣器來開發(fā)自己的隨機數(shù)生成器，生成進行簡單的已知分布。在此示例中，我們使用MH采樣器從標準雙變量正態(tài)概率分布生成隨機數(shù)。
對于這個簡單的示例，我們不需要MCMC采樣器。一種實現(xiàn)方法是使用以下代碼，該代碼從具有相關參數(shù)ρ的雙變量標準正態(tài)分布中繪制并可視化任意數(shù)量的獨立樣本。
#################＃MCMC采樣的簡單示例########################### ＃首先，讓我們構建一個從雙變量標準正態(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) ? ? # 進行跳轉newprob <- dmvnorm(c(newx,newy),sigma =# 評估跳轉ratio <- newprob/prev ? # 計算舊位置（跳出）和建議位置（跳到）的概率之比。prob.accept <- min(1,ratio) ? ? # 決定接受新跳躍的概率！if(rand<=prob.accept){x=newx;y=newy ? ?# 將x和y設置為新位置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=
?
?
讓我們嘗試解決一個問題。
MCMC對粘液瘤病進行調查
#############黏液病示例的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ù)空間內此網(wǎng)格上的似然面##############GammaLogLikelihoodFunction <- function(par}surface2D <- matrix(nrow=length(shapevec),ncol=length(scalevec)) ? #初始化存儲變量newparams <- c(shasurface2D[i,j] <- GammaLogLikelihoodFunction(newparams)############# 可視化似然面############contour(x=shapevec,y=scalevec,z=surface2D,levels=c(-30,-40,-80,-500),add=T)
這是MH算法的實現(xiàn)，用于找到后驗分布！
首先，我們需要一個似然函數(shù)–這次，我們將返回真實概率–而不是對數(shù)轉換的概率
#############編寫非對數(shù)轉換的似然函數(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-22GammaLogLike(params)## [1] -49.58983
然后，我們需要預先分配參數(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-06prior2D <- 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)
我們假設形狀和比例?在先驗中是?獨立的（聯(lián)合先驗的乘法概率）。然而，并沒有對后驗參數(shù)相關性提出相同的假設，因為概率可以反映在后驗分布中。
然后，我們需要一個函數(shù)，該函數(shù)可以計算參數(shù)空間中任何給定跳轉的后驗概率比率。因為我們正在處理?后驗概率的?比率，所以?我們不需要計算歸一化常數(shù)。
無需歸一化常數(shù)，我們只需要計算加權似然比（即先驗加權的似然比）
############# 函數(shù)用于計算參數(shù)空間中任意兩點之間的后驗密度比PosteriorRatio <- function(oldguess,newguessoldLik <- max(1e-90,GammaLikelihoodFunction(oldguess)) ? # 計算舊猜測的可能性和先驗密度newLik <- GammaLikelihoodFunction(newguess) ? ? ? ? ? ? # 在新的猜測值下計算可能性和先驗密度return((newLik*newPrior)/(oldLik*oldPrior)) ? ? ? ? ?# 計算加權似然比}PosteriorRatio2 <- function(oldguess,newguess){oldLogLik <- GammaLogLikelihoodFunction(oldguess) ? # 計算舊猜測的可能性和先驗密度newLogLik <- GammaLogLikelihoodFunction(newguess) ? ? ? ? ? ? # 在新的猜測值下計算可能性和先驗密度return(exp((newLogLik+newLogPrior)-(oldLogLik+oldLogPrior))) ? ? ? ? ?# 計算加權似然比}
## [1] 0.01436301PosteriorRatio2(oldguess,newguess)## [1] 0.01436301
然后，我們需要一個函數(shù)進行新的推測或在參數(shù)空間中跳轉：
############# 為參數(shù)空間中的跳轉定義：使用正態(tài)分布函數(shù)進行新的推測newGuess <- function(oldguess)jump <- c(shape=rnorm(1,mean=0,sd=sdshapejump),scale=rnorm(1,0,sdscalejump))newguess <- abs(oldguess + ju}# 在原始推測附近設置新的推測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中設置起點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 <- 2post.rat <- PosteriorRatio2(oldguess,newguess)prob.accept <- min(1,postoldguess <- newguessguesses[coun#可視化contour(x=shapevec,y=scal
?
我們運行更長的時間?
########### 獲取更多MCMC示例chain.length <- 100oldgucounter <- 2while(counter <= chain.length){newguess <- newGuess(oldguess)post.rat <- Posteriorand <- runif(1)if(rand<=prob.accept){ewguesscounter=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ù)的后驗分布嗎？
首先，鏈的起點“記住”起始值，因此不是固定分布。我們需要刪除鏈的第一部分。
############# 刪除預燒期（允許MCMC有一段時間達到后驗概率）burn.in <- 100MCMCsamples <- guesses[-c(1:burn.in),]
?
?
但這看起來還不是很好。讓我們運行更長的時間，看看是否得到的東西看起來更像是隨機數(shù)生成器（白噪聲）
########### 再試一次-運行更長的時間chain.length <- 20000oldguess <- startingvo2(oldguess,newguess)prob.accept <- mi
?
讓我們首先刪除前5000個樣本作為預燒期
############## 使用更長的“預燒”burn.in <- 5000MCMCsamples <- guesses[-c(1:bur
現(xiàn)在，讓我們再次看一下鏈條
?
?
在評估這些跡線圖時，我們希望看到看起來像白噪聲的“平穩(wě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）實際上傾向于結合使用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)相當簡單的已知分布。在此示例中，我們使用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)在構造一個吉布斯采樣器gibbs<-function (n, rho){ ? ?# 雙變量隨機數(shù)生成器的gibbs采樣器實現(xiàn)mat <- matrix(ncol = 2, nrow = n) ? # 用于存儲隨機樣本的矩陣mat[1, ] <- c(x, y) ? ? ? ?# initialize the markov chainfor (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)
?
在這里，馬爾可夫鏈的樣本中有很多明顯的自相關。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）來調用JAGS并后處理MCMC的輸出。JAGS相對于WinBUGS/OpenBUGS的主要優(yōu)點在于平臺的獨立性，可以應用于各種操作系統(tǒng)，而WinBUGS/OpenBUGS只能應用于windows系統(tǒng)；JAGS也可以在64-bit平臺上以64-bit應用來進行編譯。
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ù)轉換為“ 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代碼中引用的所有相關數(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)在我們可以調用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) ? # 轉換為“ 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)
評估收斂
第一步是視覺檢查-我們尋找以下內容來評估收斂性：
當視為“軌跡圖”時，每個參數(shù)的鏈應看起來像白噪聲（模糊的毛毛蟲）或類似的噪聲
多個具有不同起始條件的鏈條看起來應該相同
我們可能在這里可以做得更好的一種方法是使鏈條運行更長的時間，并丟棄初始樣本我們還可以。
通過細化鏈來嘗試減少自相關，我們每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
# 轉換為“ 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)
從視覺上看，這看起來更好。現(xiàn)在我們可以使用更多的定量收斂指標。
Gelman-Rubin診斷
一種簡單而直觀的收斂診斷程序是?Gelman-Rubin診斷程序，該診斷程序基于簡單的蒙特卡洛誤差來評估鏈之間是否比應有的鏈距更大：
############### 運行收斂診斷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診斷。如果測試失敗，則應嘗試運行更長的鏈！
所以這個模型看起來不錯！
最受歡迎的見解
1.matlab使用貝葉斯優(yōu)化的深度學習
2.matlab貝葉斯隱馬爾可夫hmm模型實現(xiàn)
3.R語言Gibbs抽樣的貝葉斯簡單線性回歸仿真
4.R語言中的block Gibbs吉布斯采樣貝葉斯多元線性回歸
5.R語言中的Stan概率編程MCMC采樣的貝葉斯模型
6.Python用PyMC3實現(xiàn)貝葉斯線性回歸模型
7.R語言使用貝葉斯 層次模型進行空間數(shù)據(jù)分析
8.R語言隨機搜索變量選擇SSVS估計貝葉斯向量自回歸（BVAR）模型
9.matlab貝葉斯隱馬爾可夫hmm模型實現(xiàn)
