R語言Rstan概率編程規(guī)劃MCMC采樣的貝葉斯模型簡介
原文http://tecdat.cn/?p=3234
?
概率編程使我們能夠實現(xiàn)統(tǒng)計模型，而無需擔心技術細節(jié)。它對基于MCMC采樣的貝葉斯模型特別有用。
?簡介
RStan是貝葉斯推理的C ++庫。它基于No-U-Turn采樣器（NUTS），用于根據(jù)用戶指定的模型和數(shù)據(jù)估計后驗分布。使用Stan執(zhí)行分析涉及以下步驟：
使用Stan建模語言指定統(tǒng)計模型。這通常通過專用的.stan文件完成。
準備要輸入模型的數(shù)據(jù)。
使用該
stan函數(shù)從后驗分布中取樣。分析結果。
在本文中，我將展示Stan使用兩個分層模型的用法。我將使用第一個模型來討論Stan的基本功能，并使用第二個示例來演示更高級的應用程序。
?數(shù)據(jù)集
?該數(shù)據(jù)集測量了教學計劃對大學入學考試的影響，即美國使用的學術能力測驗（SAT）。SAT的設計應該能夠抵抗直接針對提高測試成績的短期努力。相反，測試應反映在較長時間內獲得的知識。數(shù)據(jù)集如下：
學校估計效果?標準效果誤差?a2815b810C-316d711e-19f111g1810h1218
正如我們所看到的：對于八所學校中的大多數(shù)，短期確實增加了SAT成績。對于此數(shù)據(jù)集，我們有興趣估計與每所學校相關的真實效果大小。可以使用兩種替代方法。首先，我們可以假設所有學校都是相互獨立的。然而，這將導致難以解釋的估計，因為由于高標準誤差，學校的95％后驗間隔將在很大程度上重疊。其次，人們可以匯集所有學校的數(shù)據(jù)，假設所有學校的真實效果相同。然而，這也是不合理的，因為應該有針對學校的影響（例如不同的教師和學生）。
因此，需要另一種模型。分層模型的優(yōu)點是結合了所有八所學校（實驗）的信息，而沒有假設它們具有共同的真實效果。我們可以通過以下方式指定層次貝葉斯模型?
?
根據(jù)該模型，教學的效果遵循正態(tài)分布，其均值是真實效果，?θ? ，其標準差是?σ? ，從數(shù)據(jù)中已知。真正的效果，θ? ，遵循正態(tài)分布?μ? 和?τ 。
定義Stan模型文件
指定了我們將要使用的模型后，我們現(xiàn)在可以考慮如何在Stan中指定此模型。在為上面指定的模型定義Stan程序之前，讓我們先看看Stan建模語言的結構。
變量
在Stan中，變量可以通過以下方式定義：
int<lower=0> n; # 下限為0int<upper=5> n; # 上限為5int<lower=0,upper=5> n; # n在[0,5]中
注意，如果它們是先驗已知的，則應指定變量的上邊界和下邊界。
可以通過方括號指定多維數(shù)據(jù)：
vector[n] numbers; //長度為n的向量real[n] numbers; ?// 長度為n的浮點數(shù)組matrix[n,n] matrix; // ?n乘以n的矩陣
程序塊
Stan中使用了以下程序塊：
data：用于指定使用Bayes規(guī)則的條件
轉換數(shù)據(jù)：用于預處理數(shù)據(jù)
參數(shù)（必需）：用于指定模型的參數(shù)
變換后的參數(shù)：用于計算后驗之前的參數(shù)處理
model（必需）：用于指定模型?
生成數(shù)量：用于后處理結果
對于模型程序塊，可以以兩種等效方式指定分布。第一個，使用以下統(tǒng)計表示法：
y ~ normal(mu, sigma); # y服從正態(tài)分布
第二種方法使用基于對數(shù)概率密度函數(shù)（lpdf）的編程表示法：
?
target += normal_lpdf(y | mu, sigma); # 增加正常對數(shù)密度?
?學校的模型
現(xiàn)在我們已經了解了Stan建模langeu的基礎知識，我們可以定義我們的模型，我們將它存儲在一個名為的文件中schools.stan：
data {int<lower=0> n; //學校數(shù)}parameters {vector[n] eta; // 標準化的學校層的影響}transformed parameters {}model {}
?
?準備數(shù)據(jù)進行建模
在我們擬合模型之前，我們需要將輸入數(shù)據(jù)編碼為一個列表，其參數(shù)應該與Stan模型的數(shù)據(jù)部分中的條目相對應。?
從后驗分布中取樣
我們可以使用stan函數(shù)從后驗分布中進行采樣，執(zhí)行以下三個步驟：
它將模型規(guī)范轉換為C ++代碼。
它將C ++代碼編譯為共享對象。
它根據(jù)指定的模型，數(shù)據(jù)和設置從后驗分布中進行采樣。
?
現(xiàn)在，我們可以從后驗編譯模型和樣本。唯一需要的兩個參數(shù)stan是模型文件的位置和要輸入模型的數(shù)據(jù)。?
模型解釋
我們將首先對模型進行基本解釋，然后研究MCMC程序。
基本模型解釋
要使用擬合模型進行推理，我們可以使用該print函數(shù)。
## 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是其標準偏差。eta和theta的條目表示向量的估計η? 和?θ 。 列指示計算的值。百分比表示可信區(qū)間。例如，教練整體效果的95％可信區(qū)間，μ 是[?-?1.27?，18.26?] 。由于我們不太確定平均值，95％的可信區(qū)間為θ? 也很寬。例如，對于第一所學校，95％的可信區(qū)間是[?-?2.19?，32.33?] 。
我們可以使用plot函數(shù)可視化估算中的不確定性：
黑線表示95％的間隔，而紅線表示80％的間隔。圓圈表示平均值的估計值。
?
MCMC診斷
?通過繪制采樣程序的軌跡，我們可以確定采樣過程中是否出現(xiàn)任何問題。例如，如果鏈條在一個地方停留太長時間 。我們可以使用traceplot函數(shù)繪制模型中使用的四條鏈的痕跡：
traceplot(fit1, pars = c
?
從單個馬爾可夫鏈獲得樣本?
## ? ? ? ? ?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前端的包。?
分層回歸
現(xiàn)在我們對Stan有了基本的了解，我們可以深入了解更高級的應用程序：讓我們嘗試分層回歸。在傳統(tǒng)的回歸中，我們模擬了形式的關系
?
?=?β0+?X.β
?
該表示假定所有樣本具有相同的分布。如果存在一組樣本，那么我們就會遇到問題，因為組內和組之間的潛在差異將被忽略。
另一種方法是為每個組建立一個回歸模型。然而，在這種情況下，在估計單個模型時，小樣本量將是有問題的。
分層回歸是兩個極端之間的折衷。該模型假設這些組相似但仍然表現(xiàn)出差異。
假設每個樣本屬于其中一個?? 組。然后，分層回歸指定如下：
?
??=?α?+?X.?β（k?），????∈?{?1?，...?，?}?
?
大鼠數(shù)據(jù)集
分層回歸的典型示例是大鼠數(shù)據(jù)集。該縱向數(shù)據(jù)集包含測量5周的大鼠重量。讓我們加載數(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
?
數(shù)據(jù)顯示出不同大鼠的線性增長趨勢非常相似。然而，我們也看到大鼠具有不同的初始重量，不同的截距，以及不同的生長速率。因此，分層模型是合適的。
分層回歸模型的規(guī)范
?
我們現(xiàn)在可以指定模型并將其存儲在一個名為的文件中rats.stan：
data {int<lower=0> N; // 老鼠數(shù)int<lower=0> T; // 時間點數(shù)real y[N,T]; // 重量乘以時間的矩陣real xbar; // 時間序列中的天數(shù)的中位數(shù)}parameters {real alpha[N]; // 老鼠重量截距real mu_alpha; // 平均截距real mu_beta; //平均斜率real<lower=0> sigmasq_y;real<lower=0> sigmasq_alpha;}transformed parameters {real<lower=0> sigma_y; // 老鼠體重的標準差real<lower=0> sigma_alpha; // 截距分布的標準差sigma_y <- sqrt(sigmasq_y);sigma_alpha <- sqrt(sigmasq_alpha);sigma_beta <- sqrt(sigmasq_beta);}model {mu_alpha ~ normal(0, 100); // 非信息先驗mu_beta ~ normal(0, 100); // 非信息先驗sigmasq_alpha ~ inv_gamma(0.001, 0.001); // 普通的共軛先驗sigmasq_beta ~ inv_gamma(0.001, 0.001); // 普通的共軛先驗for (n in 1:N) // 對于每個樣本for (t in 1:T) ?// 每個時間點y[n,t] ~ normal(alpha[n] + beta[n] * (x[t] - xbar), sigma_y);}generated quantities {// 確定時間0（出生體重）的截距alpha0 <- mu_alpha - xbar * mu_beta;}
請注意，模型代碼估計方差（sigmasq變量）而不是標準偏差。此外，時間0的截距，即出生時大鼠的體重。我們還可以計算其他數(shù)量，例如，不同時間點的大鼠的估計重量。我們稍后會在R中執(zhí)行此操作。
數(shù)據(jù)準備
要為模型準備數(shù)據(jù)，我們首先將測量點提取為數(shù)值，然后在列表結構中對所有內容進行編碼：
擬合回歸模型
我們現(xiàn)在可以擬合大鼠體重數(shù)據(jù)集的貝葉斯分層回歸模型：
用層次回歸模型預測
確定了?α 和?β 對于每只大鼠，我們現(xiàn)在可以在任意時間點估計個體大鼠的體重。在這里，我們感興趣的是從第0天到第100天找到大鼠的體重。
?
與原始數(shù)據(jù)相比，模型的估計是平滑的，因為每條曲線都遵循線性模型。研究最后一個圖中顯示的置信區(qū)間，我們可以看出方差估計是合理的。離采樣區(qū)域越遠，不確定性越大。
?
?
非常感謝您閱讀本文，有任何問題請在下方留言！
