視頻：R語言中的Stan概率編程MCMC采樣的貝葉斯模型
原文鏈接：http://tecdat.cn/?p=11161
概率編程使我們能夠?qū)崿F(xiàn)統(tǒng)計(jì)模型，而不必?fù)?dān)心技術(shù)細(xì)節(jié)。這對于基于MCMC采樣的貝葉斯模型特別有用。?
R語言中RStan貝葉斯層次模型分析示例
stan簡介
Stan是用于貝葉斯推理的C ++庫。它基于No-U-Turn采樣器（NUTS），該采樣器用于根據(jù)用戶指定的模型和數(shù)據(jù)估計(jì)后驗(yàn)分布。使用Stan執(zhí)行分析涉及以下步驟：
使用Stan建模語言指定統(tǒng)計(jì)模型。通過專用的.stan??文件完成此操作? 。
準(zhǔn)備要提供給模型的數(shù)據(jù)。
使用該
stan?函數(shù)從后驗(yàn)分布中采樣? 。分析結(jié)果。
在本文中，我將通過兩個(gè)層次模型展示Stan的用法。我將使用第一個(gè)模型討論Stan的基本功能，并使用第二個(gè)示例演示更高級的應(yīng)用。
?學(xué)校數(shù)據(jù)集
我們要使用的第一個(gè)數(shù)據(jù)集是??學(xué)校的數(shù)據(jù)集? 。該數(shù)據(jù)集衡量了教練計(jì)劃對大學(xué)入學(xué)考試（在美國使用的學(xué)業(yè)能力測驗(yàn)（SAT））的影響。 數(shù)據(jù)集如下所示：
正如我們所看到的：對于八所學(xué)校中的大多數(shù)，短期教練計(jì)劃的確提高了SAT分?jǐn)?shù) 。對于此數(shù)據(jù)集，我們有興趣估算與每所學(xué)校相關(guān)的真實(shí)教練計(jì)劃效果大小。我們考慮兩種替代方法。首先，我們可以假設(shè)所有學(xué)校彼此獨(dú)立。但是，這將難以解釋，因?yàn)閷W(xué)校的后驗(yàn)區(qū)間由于高標(biāo)準(zhǔn)差而在很大程度上重疊。第二，假設(shè)所有學(xué)校的真實(shí)效果都相同，則可以匯總所有學(xué)校的數(shù)據(jù)。但是，這也是不合理的，因?yàn)樵撚?jì)劃有針對學(xué)校的不同效果（例如，不同的老師和學(xué)生應(yīng)該有不同的計(jì)劃）。
因此，需要另一個(gè)模型。分層模型的優(yōu)點(diǎn)是可以合并來自所有八所學(xué)校的信息，而無需假定它們具有共同的真實(shí)效果。我們可以通過以下方式指定層次貝葉斯模型：
根據(jù)該模型，教練的效果遵循正態(tài)分布，其均值是真實(shí)效果θ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àn)已知變量，則應(yīng)指定變量的上下邊界。
多維數(shù)據(jù)可以通過方括號指定：
vector[n] numbers; // 長度為n的向量real[n] numbers; ?// 長度為n的浮點(diǎ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ù)：用于計(jì)算后驗(yàn)之前的參數(shù)處理
模型??（必填）：用于指定模型
生成數(shù)量：用于對結(jié)果進(jìn)行后處理
對于??模型??程序塊，可以兩種等效方式指定分布。第一個(gè)，使用以下統(tǒng)計(jì)符號：
y ~ normal(mu, sigma); # y 服從正態(tài)分布
第二種方法使用基于對數(shù)概率密度函數(shù)（lpdf）的程序化表示法：
target += normal_lpdf(y | mu, sigma); # 增加正態(tài)對數(shù)密度
Stan支持大量的概率分布。通過Stan指定模型時(shí)，該??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?等價(jià)于?Stan的?normal_rng?。
模型
現(xiàn)在，我們了解了Stan建模語言的基礎(chǔ)知識，我們可以定義模型，并將其存儲在一個(gè)名為的文件中??schools.stan：
注意，θ 永遠(yuǎn)不會出現(xiàn)在參數(shù)中。這是因?yàn)槲覀儧]有顯式地對θ進(jìn)行建模，而是對η（各個(gè)學(xué)校的標(biāo)準(zhǔn)化效果）進(jìn)行了建模。然后，?根據(jù)μ，τ和η在變換后的參數(shù)部分構(gòu)造θ? 。此參數(shù)化使采樣器更高效。
?
準(zhǔn)備數(shù)據(jù)進(jìn)行建模
在擬合模型之前，我們需要將輸入數(shù)據(jù)編碼為一個(gè)列表，其參數(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))
從后驗(yàn)分布抽樣
我們可以使用stan?函數(shù)從后驗(yàn)分布中采樣，函數(shù)執(zhí)行以下三個(gè)步驟：
它將模型規(guī)范轉(zhuǎn)換為C ++代碼。
它將C ++代碼編譯為共享對象。
它根據(jù)指定的模型，數(shù)據(jù)和設(shè)置從后驗(yàn)分布中采樣。
如果??rstan_options(auto_write = TRUE)，則相同模型的后續(xù)調(diào)用將比第一次調(diào)用快得多，因?yàn)樵??stan?函數(shù)隨后跳過了前兩個(gè)步驟（轉(zhuǎn)換和編譯模型）。此外，我們將設(shè)置要使用的內(nèi)核數(shù)：
options(mc.cores = parallel::detectCores()) # 并行化rstan_options(auto_write = TRUE) ?# 存儲編譯的stan模型
現(xiàn)在，我們可以從后驗(yàn)中編譯模型和樣本。
模型解釋
我們將首先對模型進(jì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).
在此，行名稱表示估計(jì)的參數(shù)：mu是后驗(yàn)分布的平均值，而tau是其標(biāo)準(zhǔn)偏差。eta和theta的條目分別表示矢量η和θ的估計(jì)值。這些列表示計(jì)算值。百分比表示置信區(qū)間。例如，教練計(jì)劃的總體效果的95％可信區(qū)間μ為[-1.27,18.26]。由于我們不確定平均值，因此θj的95％置信區(qū)間也很寬。例如，對于第一所學(xué)校，95％置信區(qū)間為[?2.19,32.33]。
我們可以使用以下plot?函數(shù)來可視化估計(jì)中的不確定性? ：
黑線表示95％的間隔，而紅線表示80％的間隔。圓圈表示平均值的估計(jì)。
我們可以使用以下extract?函數(shù)獲取生成的樣本? ：
# 獲取樣本samples <- extract(fit1, permuted = TRUE) # 每個(gè)參數(shù)1000個(gè)樣本
MCMC診斷
?通過繪制采樣過程的軌跡圖，我們可以確定采樣期間是否出了問題。例如，鏈條在一個(gè)位置停留的時(shí)間過長或在一個(gè)方向上走了太多步，就會有問題。我們可以使用traceplot?函數(shù)繪制模型中使用的四個(gè)鏈的軌跡? ：
# 診斷:
?
要從各個(gè)馬爾可夫鏈中獲取樣本，我們可以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
為了對采樣過程進(jìn)行更高級的分析，我們可以使用該??shinystan?軟件包 。使用該軟件包，可以通過以下方式啟動(dòng)Shiny應(yīng)用程序來分析擬合模型：
library(shinystan)launch_shinystan(fit1)
層次回歸
現(xiàn)在，我們對Stan有了基本的了解，我們可以深入研究更高級的應(yīng)用程序：讓我們嘗試一下層次回歸。在常規(guī)回歸中，我們對以下形式的關(guān)系進(jìn)行建模
此表示假設(shè)所有樣本都具有相同的分布。如果只存在一組樣本，那么我們就會遇到問題，因?yàn)閷⒑雎越M內(nèi)和組之間的潛在差異。
另一種選擇是為每個(gè)組建立一個(gè)回歸模型。但是，在這種情況下，估計(jì)單個(gè)模型時(shí)，小樣本量會帶來問題。
層次回歸是兩個(gè)極端之間的折衷。該模型假設(shè)組是相似的，但存在差異。
假設(shè)每個(gè)樣本都屬于K組之一。然后，層次回歸指定如下：
其中Yk是第k組的結(jié)果，αk是截距，Xk是特征，β（k）表示權(quán)重。層次模型不同于其中Yk分別擬合每個(gè)組的模型，因?yàn)榧俣▍?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個(gè)大鼠的截距由αi表示，斜率由βi表示。注意，測量時(shí)間的中心是x = 22，它是時(shí)間序列數(shù)據(jù)的中值測量值（第22天）。
現(xiàn)在，我們可以指定模型并將其存儲在名為?rats.stan的文件中?：
請注意，模型代碼估算的是方差（??sigmasq??變量）而不是標(biāo)準(zhǔn)差。?
資料準(zhǔn)備
為了準(zhǔn)備模型數(shù)據(jù)，我們首先將測量點(diǎn)提取為數(shù)值，然后將所有內(nèi)容編碼為列表結(jié)構(gòu)：
data <- list(N = nrow(df), T = ncol(df), x = days,y = df, xbar = median(days))
擬合回歸模型
現(xiàn)在，我們可以為老鼠體重?cái)?shù)據(jù)集擬合貝葉斯層次回歸模型：
# 模型包含截距（alpha）和斜率（beta）的估計(jì)
層次回歸模型的預(yù)測
在確定了每只大鼠的α和β之后，我們現(xiàn)在可以估計(jì)任意時(shí)間點(diǎn)單個(gè)大鼠的體重。在這里，我們尋找從第0天到第100天的大鼠體重。
?
ggplot(pred.df[pred.df$Rat %in% sel.rats, ],aes(x = Day, y = Weight, group = Rat,geom_line() ?+
與原始數(shù)據(jù)相比，該模型的估計(jì)是平滑的，因?yàn)槊織l曲線都遵循線性模型。研究最后一個(gè)圖中所示的置信區(qū)間，我們可以看到方差估計(jì)是合理的。我們對采樣時(shí)（第8至36天）的老鼠體重充滿信心，但是隨著離開采樣區(qū)域，不確定性會增加。
最受歡迎的見解
1.matlab使用貝葉斯優(yōu)化的深度學(xué)習(xí)
2.matlab貝葉斯隱馬爾可夫hmm模型實(shí)現(xiàn)
3.R語言Gibbs抽樣的貝葉斯簡單線性回歸仿真
4.R語言中的block Gibbs吉布斯采樣貝葉斯多元線性回歸
5.R語言中的Stan概率編程MCMC采樣的貝葉斯模型
6.Python用PyMC3實(shí)現(xiàn)貝葉斯線性回歸模型
7.R語言使用貝葉斯 層次模型進(jìn)行空間數(shù)據(jù)分析
8.R語言隨機(jī)搜索變量選擇SSVS估計(jì)貝葉斯向量自回歸（BVAR）模型
9.matlab貝葉斯隱馬爾可夫hmm模型實(shí)現(xiàn)
