R語言中的Stan概率編程MCMC采樣的貝葉斯模型|附代碼數(shù)據(jù)
原文鏈接：http://tecdat.cn/?p=11161
最近我們被客戶要求撰寫關(guān)于貝葉斯模型的研究報告，包括一些圖形和統(tǒng)計輸出。
概率編程使我們能夠?qū)崿F(xiàn)統(tǒng)計模型，而不必?fù)?dān)心技術(shù)細(xì)節(jié)。這對于基于MCMC采樣的貝葉斯模型特別有用
R語言中RStan貝葉斯層次模型分析示例
stan簡介
Stan是用于貝葉斯推理的C ++庫。它基于No-U-Turn采樣器（NUTS），該采樣器用于根據(jù)用戶指定的模型和數(shù)據(jù)估計后驗分布。使用Stan執(zhí)行分析涉及以下步驟：
使用Stan建模語言指定統(tǒng)計模型。通過專用的_.stan_ ?文件完成此操作? 。
準(zhǔn)備要提供給模型的數(shù)據(jù)。
使用該
stan?函數(shù)從后驗分布中采樣? 。分析結(jié)果。
在本文中，我將通過兩個層次模型展示Stan的用法。我將使用第一個模型討論Stan的基本功能，并使用第二個示例演示更高級的應(yīng)用。
?學(xué)校數(shù)據(jù)集
我們要使用的第一個數(shù)據(jù)集是??學(xué)校的數(shù)據(jù)集? 。該數(shù)據(jù)集衡量了教練計劃對大學(xué)入學(xué)考試（在美國使用的學(xué)業(yè)能力測驗（SAT））的影響。數(shù)據(jù)集如下所示：
正如我們所看到的：對于八所學(xué)校中的大多數(shù)，短期教練計劃的確提高了SAT分?jǐn)?shù) 。對于此數(shù)據(jù)集，我們有興趣估算與每所學(xué)校相關(guān)的真實教練計劃效果大小。我們考慮兩種替代方法。首先，我們可以假設(shè)所有學(xué)校彼此獨立。但是，這將難以解釋，因為學(xué)校的后驗區(qū)間由于高標(biāo)準(zhǔn)差而在很大程度上重疊。第二，假設(shè)所有學(xué)校的真實效果都相同，則可以匯總所有學(xué)校的數(shù)據(jù)。但是，這也是不合理的，因為該計劃有針對學(xué)校的不同效果（例如，不同的老師和學(xué)生應(yīng)該有不同的計劃）。
因此，需要另一個模型。分層模型的優(yōu)點是可以合并來自所有八所學(xué)校的信息，而無需假定它們具有共同的真實效果。我們可以通過以下方式指定層次貝葉斯模型：
根據(jù)該模型，教練的效果遵循正態(tài)分布，其均值是真實效果θj，其標(biāo)準(zhǔn)偏差為σj（從數(shù)據(jù)中得知）。真正的影響θj遵循參數(shù)μ和τ的正態(tài)分布。
定義Stan模型文件
在指定了要使用的模型之后，我們現(xiàn)在可以討論如何在Stan中指定此模型。在為上述模型定義Stan程序之前，讓我們看一下Stan建模語言的結(jié)構(gòu)。
變量
在Stan中，可以通過以下方式定義變量：
int<lower=0>?n;?#?下界是0int<upper=5>?n;?#?上限是5int<lower=0,upper=5>?n;?#?n?的范圍是?[0,5]
注意，如果先驗已知變量，則應(yīng)指定變量的上下邊界。
多維數(shù)據(jù)可以通過方括號指定：
vector[n]?numbers;?//?長度為n的向量
real[n]?numbers;??//?長度為n的浮點數(shù)組
matrix[n,n]?matrix;?//?n乘n矩陣
程序?
Stan中使用以下程序 ：
data：用于指定以貝葉斯規(guī)則為條件的數(shù)據(jù)
轉(zhuǎn)換后的數(shù)據(jù)：用于預(yù)處理數(shù)據(jù)
參數(shù)??（必填）：用于指定模型的參數(shù)
轉(zhuǎn)換后的參數(shù)：用于計算后驗之前的參數(shù)處理
模型??（必填）：用于指定模型
生成數(shù)量：用于對結(jié)果進行后處理
點擊標(biāo)題查閱往期內(nèi)容
MCMC的rstan貝葉斯回歸模型和標(biāo)準(zhǔn)線性回歸模型比較
左右滑動查看更多
01
02
03
04
對于??模型??程序塊，可以兩種等效方式指定分布。第一個，使用以下統(tǒng)計符號：
y?~?normal(mu,?sigma);?#?y?服從正態(tài)分布
第二種方法使用基于對數(shù)概率密度函數(shù)（lpdf）的程序化表示法：
target?+=?normal_lpdf(y?|?mu,?sigma);?#?增加正態(tài)對數(shù)密度
Stan支持大量的概率分布。通過Stan指定模型時，該??lookup?函數(shù)會派上用場：它提供從R函數(shù)到Stan函數(shù)的映射?？紤]以下示例：
library(rstan)?#?加載stan包lookup(rnorm)##?????StanFunction?????????????Arguments?ReturnType?Page
##?355???normal_rng?(real?mu,?real?sigma)???????real??494
在這里，我們看到R中的rnorm?等價于?Stan的?normal_rng?。
模型
現(xiàn)在，我們了解了Stan建模語言的基礎(chǔ)知識，我們可以定義模型，并將其存儲在一個名為的文件中??schools.stan：
注意，θ 永遠(yuǎn)不會出現(xiàn)在參數(shù)中。這是因為我們沒有顯式地對θ進行建模，而是對η（各個學(xué)校的標(biāo)準(zhǔn)化效果）進行了建模。然后，?根據(jù)μ，τ和η在_變換后的參數(shù)_部分構(gòu)造θ? 。此參數(shù)化使采樣器更高效。
準(zhǔn)備數(shù)據(jù)進行建模
在擬合模型之前，我們需要將輸入數(shù)據(jù)編碼為一個列表，其參數(shù)應(yīng)與Stan模型的數(shù)據(jù)部分相對應(yīng)。對于學(xué)校數(shù)據(jù)，數(shù)據(jù)如下：
schools.data?<-?list(
??n?=?8,
??y?=?c(28,??8,?-3,??7,?-1,??1,?18,?12),
??sigma?=?c(15,?10,?16,?11,??9,?11,?10,?18)
)
從后驗分布抽樣
我們可以使用stan?函數(shù)從后驗分布中采樣，函數(shù)執(zhí)行以下三個步驟：
它將模型規(guī)范轉(zhuǎn)換為C ++代碼。
它將C ++代碼編譯為共享對象。
它根據(jù)指定的模型，數(shù)據(jù)和設(shè)置從后驗分布中采樣。
如果??rstan_options(auto_write = TRUE)，則相同模型的后續(xù)調(diào)用將比第一次調(diào)用快得多，因為該??stan?函數(shù)隨后跳過了前兩個步驟（轉(zhuǎn)換和編譯模型）。此外，我們將設(shè)置要使用的內(nèi)核數(shù)：
options(mc.cores?=?parallel::detectCores())?#?并行化rstan_options(auto_write?=?TRUE)??#?存儲編譯的stan模型
現(xiàn)在，我們可以從后驗中編譯模型和樣本。
模型解釋
我們將首先對模型進行基本解釋，然后研究MCMC程序。
基本模型解釋
要使用擬合模型執(zhí)行推斷，我們可以使用??print?函數(shù)。
print(fit1)?#?可選參數(shù)：pars，probs##?Inference?for?Stan?model:?schools.
##?4?chains,?each?with?iter=2000;?warmup=1000;?thin=1;?
##?post-warmup?draws?per?chain=1000,?total?post-warmup?draws=4000.##?
##????????????mean?se_mean???sd???2.5%????25%????50%????75%??97.5%?n_eff?Rhat
##?mu?????????7.67????0.15?5.14??-2.69???4.42???7.83??10.93??17.87??1185????1##?tau????????6.54????0.16?5.40???0.31???2.52???5.28???9.05??20.30??1157????1##?eta[1]?????0.42????0.01?0.92??-1.47??-0.18???0.44???1.03???2.18??4000????1##?eta[2]?????0.03????0.01?0.87??-1.74??-0.54???0.03???0.58???1.72??4000????1##?eta[3]????-0.18????0.02?0.92??-1.95??-0.81??-0.20???0.45???1.65??3690????1##?eta[4]????-0.03????0.01?0.92??-1.85??-0.64??-0.02???0.57???1.81??4000????1##?eta[5]????-0.33????0.01?0.86??-2.05??-0.89??-0.34???0.22???1.43??3318????1##?eta[6]????-0.20????0.01?0.87??-1.91??-0.80??-0.21???0.36???1.51??4000????1##?eta[7]?????0.37????0.02?0.87??-1.37??-0.23???0.37???0.96???2.02??3017????1##?eta[8]?????0.05????0.01?0.92??-1.77??-0.55???0.05???0.69???1.88??4000????1##?theta[1]??11.39????0.15?8.09??-2.21???6.14??10.30??15.56??30.22??2759????1##?theta[2]???7.92????0.10?6.25??-4.75???4.04???8.03??11.83??20.05??4000????1##?theta[3]???6.22????0.14?7.83?-11.41???2.03???6.64??10.80??20.97??3043????1##?theta[4]???7.58????0.10?6.54??-5.93???3.54???7.60??11.66??20.90??4000????1##?theta[5]???5.14????0.10?6.30??-8.68???1.40???5.63???9.50??16.12??4000????1##?theta[6]???6.08????0.10?6.62??-8.06???2.21???6.45??10.35??18.53??4000????1##?theta[7]??10.60????0.11?6.70??-0.94???6.15??10.01??14.48??25.75??4000????1##?theta[8]???8.19????0.14?8.18??-8.13???3.59???8.01??12.48??25.84??3361????1##?lp__?????-39.47????0.07?2.58?-45.21?-41.01?-39.28?-37.70?-34.99??1251????1##?
##?Samples?were?drawn?using?NUTS(diag_e)?at?Thu?Nov?29?11:17:50?2018.##?For?each?parameter,?n_eff?is?a?crude?measure?of?effective?sample?size,
##?and?Rhat?is?the?potential?scale?reduction?factor?on?split?chains?(at?
##?convergence,?Rhat=1).
在此，行名稱表示估計的參數(shù)：mu是后驗分布的平均值，而tau是其標(biāo)準(zhǔn)偏差。eta和theta的條目分別表示矢量η和θ的估計值。這些列表示計算值。百分比表示置信區(qū)間。例如，教練計劃的總體效果的95％可信區(qū)間μ為[-1.27,18.26]。由于我們不確定平均值，因此θj的95％置信區(qū)間也很寬。例如，對于第一所學(xué)校，95％置信區(qū)間為[?2.19,32.33]。
我們可以使用以下plot?函數(shù)來可視化估計中的不確定性? ：
黑線表示95％的間隔，而紅線表示80％的間隔。圓圈表示平均值的估計。
我們可以使用以下extract?函數(shù)獲取生成的樣本? ：
#?獲取樣本samples?<-?extract(fit1,?permuted?=?TRUE)?#?每個參數(shù)1000個樣本
MCMC診斷
通過繪制采樣過程的軌跡圖，我們可以確定采樣期間是否出了問題。例如，鏈條在一個位置停留的時間過長或在一個方向上走了太多步，就會有問題。我們可以使用traceplot?函數(shù)繪制模型中使用的四個鏈的軌跡? ：
#?診斷:
要從各個馬爾可夫鏈中獲取樣本，我們可以extract?再次使用函數(shù)：
##??????????parameters
##?chains???????????mu???????tau?????eta[1]?????eta[2]?????eta[3]?????eta[4]
##???chain:1??1.111120??2.729124?-0.1581242?-0.8498898??0.5025965?-1.9874554##???chain:2??3.633421??2.588945??1.2058772?-1.1173221??1.4830778??0.4838649##???chain:3?13.793056??3.144159??0.6023924?-1.1188243?-1.2393491?-0.6118482##???chain:4??3.673380?13.889267?-0.0869434??1.1900236?-0.0378830?-0.2687284##??????????parameters
##?chains????????eta[5]?????eta[6]?????eta[7]??????eta[8]???theta[1]
##???chain:1??0.3367602?-1.1940843??0.5834020?-0.08371249??0.6795797##???chain:2?-1.8057252??0.7429594??0.9517675??0.55907356??6.7553706##???chain:3?-1.5867789??0.6334288?-0.4613463?-1.44533007?15.6870727##???chain:4??0.1028605??0.3481214??0.9264762??0.45331024??2.4657999##??????????parameters
##?chains?????theta[2]?theta[3]????theta[4]??theta[5]??theta[6]??theta[7]
##???chain:1?-1.208335?2.482769?-4.31289292??2.030181?-2.147684??2.703297##???chain:2??0.740736?7.473028??4.88612054?-1.041502??5.556902??6.097494##???chain:3?10.275294?9.896345?11.86930758??8.803971?15.784656?12.342510##???chain:4?20.201935?3.147213?-0.05906019??5.102037??8.508530?16.541455##??????????parameters
##?chains?????theta[8]??????lp__
##???chain:1?0.8826584?-41.21499##???chain:2?5.0808317?-41.17178##???chain:3?9.2487083?-40.35351##???chain:4?9.9695268?-36.34043
為了對采樣過程進行更高級的分析，我們可以使用該??shinystan?軟件包 。使用該軟件包，可以通過以下方式啟動Shiny應(yīng)用程序來分析擬合模型：
library(shinystan)launch_shinystan(fit1)
層次回歸
現(xiàn)在，我們對Stan有了基本的了解，我們可以深入研究更高級的應(yīng)用程序：讓我們嘗試一下層次回歸。在常規(guī)回歸中，我們對以下形式的關(guān)系進行建模
此表示假設(shè)所有樣本都具有相同的分布。如果只存在一組樣本，那么我們就會遇到問題，因為將忽略組內(nèi)和組之間的潛在差異。
另一種選擇是為每個組建立一個回歸模型。但是，在這種情況下，估計單個模型時，小樣本量會帶來問題。
層次回歸是兩個極端之間的折衷。該模型假設(shè)組是相似的，但存在差異。
假設(shè)每個樣本都屬于K組之一。然后，層次回歸指定如下：
其中Yk是第k組的結(jié)果，αk是截距，Xk是特征，β（k）表示權(quán)重。層次模型不同于其中Yk分別擬合每個組的模型，因為假定參數(shù)αk和β（k）源自共同的分布。
?數(shù)據(jù)集
分層回歸的經(jīng)典示例是 老鼠數(shù)據(jù)集。該數(shù)據(jù)集包含5周內(nèi)測得的 鼠體重。讓我們加載數(shù)據(jù)：
##???day8?day15?day22?day29?day36
##?1??151???199???246???283???320##?2??145???199???249???293???354##?3??147???214???263???312???328##?4??155???200???237???272???297##?5??135???188???230???280???323##?6??159???210???252???298???331
讓我們調(diào)查數(shù)據(jù)：
library(ggplot2)ggplot(ddf,?aes(x?=?variable,?y?=?value,?group?=?Group))?+?geom_line()?+?geom_point()
數(shù)據(jù)顯示線性增長趨勢對于不同的大鼠非常相似。但是，我們還看到，大鼠的初始體重不同，需要不同的截距，并且生長速度也需要不同的斜率。因此，分層模型似乎是適當(dāng)?shù)摹?/p>
層次回歸模型的規(guī)范
該模型可以指定如下：
第i個大鼠的截距由αi表示，斜率由βi表示。注意，測量時間的中心是x = 22，它是時間序列數(shù)據(jù)的中值測量值（第22天）。
現(xiàn)在，我們可以指定模型并將其存儲在名為?rats.stan的文件中?：
請注意，模型代碼估算的是方差（??sigmasq??變量）而不是標(biāo)準(zhǔn)差。
資料準(zhǔn)備
為了準(zhǔn)備模型數(shù)據(jù)，我們首先將測量點提取為數(shù)值，然后將所有內(nèi)容編碼為列表結(jié)構(gòu)：
data?<-?list(N?=?nrow(df),?T?=?ncol(df),?x?=?days,
?????????????????y?=?df,?xbar?=?median(days))
擬合回歸模型
現(xiàn)在，我們可以為老鼠體重數(shù)據(jù)集擬合貝葉斯層次回歸模型：
#?模型包含截距（alpha）和斜率（beta）的估計
層次回歸模型的預(yù)測
在確定了每只大鼠的α和β之后，我們現(xiàn)在可以估計任意時間點單個大鼠的體重。在這里，我們尋找從第0天到第100天的大鼠體重。
ggplot(pred.df[pred.df$Rat?%in%?sel.rats,?],?
???????aes(x?=?Day,?y?=?Weight,?group?=?Rat,?
????geom_line()??+
與原始數(shù)據(jù)相比，該模型的估計是平滑的，因為每條曲線都遵循線性模型。研究最后一個圖中所示的置信區(qū)間，我們可以看到方差估計是合理的。我們對采樣時（第8至36天）的老鼠體重充滿信心，但是隨著離開采樣區(qū)域，不確定性會增加。
點擊文末?“閱讀原文”
獲取全文完整資料。
本文選自《R語言中的Stan概率編程MCMC采樣的貝葉斯模型》。
點擊標(biāo)題查閱往期內(nèi)容
R語言貝葉斯MCMC：用rstan建立線性回歸模型分析汽車數(shù)據(jù)和可視化診斷
【視頻】馬爾可夫鏈蒙特卡羅方法MCMC原理與R語言實現(xiàn)|數(shù)據(jù)分享
R語言實現(xiàn)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貝葉斯回歸分析住房負(fù)擔(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的貝葉斯分層混合模型的診斷準(zhǔn)確性研究
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采樣用于回歸的貝葉斯估計
