R語言中的Stan概率編程MCMC采樣的貝葉斯模型

 

原文鏈接：http://tecdat.cn/?p=11161

 

 

概率編程使我們能夠實現統計模型，而不必擔心技術細節。這對於基於MCMC采樣的貝葉斯模型特別有用。 

stan簡介

Stan是用於貝葉斯推理的C ++庫。它基於No-U-Turn采樣器（NUTS），該采樣器用於根據用戶指定的模型和數據估計后驗分布。使用Stan執行分析涉及以下步驟：

1. 使用Stan建模語言指定統計模型。通常通過專用的.stan  文件完成此操作  。
2. 准備要提供給模型的數據。
3. 使用該stan 函數從后驗分布中采樣  。
4. 分析結果。

在本文中，我將通過兩個層次模型展示Stan的用法。我將使用第一個模型討論Stan的基本功能，並使用第二個示例演示更高級的應用程序。

 學校數據集

我們要使用的第一個數據集是  學校的數據集  。該數據集衡量了教練計划對大學入學考試（在美國使用的學業能力測驗（SAT））的影響。 數據集如下所示：

​

 正如我們所看到的，SAT並沒有兌現其諾言：對於八所學校中的大多數，短期教練的確提高了SAT分數 。對於此數據集，我們有興趣估算與每所學校相關的真實效果大小。可以使用兩種替代方法。首先，我們可以假設所有學校彼此獨立。但是，這將導致難以解釋的估計，因為學校的95％的后驗間隔由於高標准誤差而在很大程度上重疊。第二，假設所有學校的真實效果都相同，則可以匯總所有學校的數據。但是，這也是不合理的，因為應該有針對學校的效果（例如，不同的老師和學生）。

因此，需要另一個模型。分層模型的優點是可以合並來自所有八所學校（實驗）的信息，而無需假定它們具有共同的真實效果。我們可以通過以下方式指定層次貝葉斯模型：

 

​

根據該模型，教練的效果遵循正態分布，其均值是真實效果θjθj，其標准偏差為σjσj（從數據中得知）。真正的影響θjθj遵循參數μμ和ττ的正態分布。

定義Stan模型文件

在指定了要使用的模型之后，我們現在可以討論如何在Stan中指定此模型。在為上述模型定義Stan程序之前，讓我們看一下Stan建模語言的結構。

變量

在Stan中，可以通過以下方式定義變量：

int<lower=0> n; # lower bound is 0
int<upper=5> n; # upper bound is 5
int<lower=0,upper=5> n; # n is in [0,5]

注意，如果先驗已知變量，則應指定變量的上下邊界。

多維數據可以通過方括號指定：

vector[n] numbers; // a vector of length n
real[n] numbers;  // an array of floats with length n
matrix[n,n] matrix; // an n times n matrix

程序 

Stan中使用以下程序 ：

• data：用於指定以貝葉斯規則為條件的數據
• 轉換后的數據：用於預處理數據
• 參數  （必填）：用於指定模型的參數
• 轉換后的參數：用於計算后驗之前的參數處理
• 模型  （必填）：用於指定模型本身
• 生成數量：用於對結果進行后處理

對於  模型  程序塊，可以兩種等效方式指定分布。第一個，使用以下統計符號：

y ~ normal(mu, sigma); # y follows a normal distribution

第二種方法使用基於對數概率密度函數（lpdf）的程序化表示法：

target += normal_lpdf(y | mu, sigma); # increment the normal log density

Stan支持大量的概率分布。通過Stan指定模型時，該  lookup 函數會派上用場：它提供從R函數到Stan函數的映射。考慮以下示例：

library(rstan) # load stan package
lookup(rnorm)

##     StanFunction             Arguments ReturnType Page
## 355   normal_rng (real mu, real sigma)       real  494

在這里，我們看到rnorm R  中的等價  normal_rng 於Stan。

模型

現在，我們了解了Stan建模語言的基礎知識，我們可以定義模型，並將其存儲在一個名為的文件中  schools.stan：

注意，θ 永遠不會出現在參數中。這是因為我們沒有顯式地對θθ進行建模，而是對ηη（各個學校的標准化效果）進行了建模。然后， 根據μμ，ττ和ηη在變換后的參數部分構造θθ  。此參數化使采樣器更高效。

 

准備數據進行建模

在適合模型之前，我們需要將輸入數據編碼為一個列表，其參數應與Stan模型的數據部分中的條目相對應。對於學校數據，數據如下：

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 函數從后驗分布中采樣，該  函數執行以下三個步驟：

1. 它將模型規范轉換為C ++代碼。
2. 它將C ++代碼編譯為共享對象。
3. 它根據指定的模型，數據和設置從后驗分布中采樣。

如果  rstan_options(auto_write = TRUE) 已在活動的R會話中執行，則相同模型的后續調用將比第一次調用快得多，因為該  stan 函數隨后跳過了前兩個步驟（轉換和編譯模型）。此外，我們將設置要使用的內核數：

options(mc.cores = parallel::detectCores()) # parallelize
rstan_options(auto_write = TRUE)  # store compiled stan model

現在，我們可以從后驗中編譯模型和樣本。唯一需要的兩個參數  stan 是模型文件的位置和要饋送到模型的數據。此處顯示的其他參數僅用於概述可能的選項：

模型解釋

我們將首先對模型進行基本解釋，然后研究MCMC程序。

基本模型解釋

要使用擬合模型執行推斷，我們可以使用  print 函數。

print(fit1) # optional parameters: 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).

在此，行名稱表示估計的參數：mu是后驗分布的平均值，而tau是其標准偏差。eta和theta的條目分別表示矢量ηη和θθ的估計值。該  LP  條目顯示日志密度，這通常是不相關的。這些列表示計算值。百分比表示可靠的時間間隔。例如，指導的總體效果的95％可信區間μμ為[-1.27,18.26] [-1.27,18.26]。由於我們不確定平均值，因此θjθj的95％可信區間也很寬。例如，對於第一所學校，95％可信區間為[−2.19,32.33] [− 2.19,32.33]。

我們可以使用以下plot 函數來可視化估計中的不確定性  ：

​

黑線表示95％的間隔，而紅線表示80％的間隔。圓圈表示平均值的估計。

我們可以使用以下extract 函數獲取生成的樣本  ：

# retrieve the samples
samples <- extract(fit1, permuted = TRUE) # 1000 samples per parameter
mu <- samples$mu  # samples of mu only

MCMC診斷

 通過繪制采樣過程的軌跡圖，我們可以確定采樣期間是否出了問題。例如，如果鏈條在一個位置停留的時間過長或在一個方向上走了太多步，就是這種情況。我們可以使用以下traceplot 函數繪制模型中使用的四個鏈的軌跡  ：

# diagnostics:
traceplot(fit1, pars = c("mu", "tau"), inc_warmup = TRUE, nrow = 2)

​

 

要從各個馬爾可夫鏈中獲取樣本，我們可以extract 再次使用該  函數：

##          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應用程序來分析擬合模型：

library(shinystan)
launch_shinystan(fit1)

層次回歸

現在，我們對Stan有了基本的了解，我們可以深入研究更高級的應用程序：讓我們嘗試一下層次回歸。在常規回歸中，我們對以下形式的關系進行建模

 

​

此表示假設所有樣本都具有相同的分布。如果存在一組樣本，那么我們就會遇到問題，因為將忽略組內和組之間的潛在差異。

另一種選擇是為每個組建立一個回歸模型。但是，在這種情況下，估計單個模型時，小樣本量會帶來問題。

層次回歸是兩個極端之間的折衷。該模型假設組是相似的，但存在差異。

假設每個樣本都屬於KK組之一。然后，層次回歸指定如下：

​

其中YkYk是第kk組的結果，αkαk是截距，XkXk是特征，β（k）β（k）表示權重。層次模型不同於其中YkYk分別適合每個組的模型，因為假定參數αkαk和β（k）β（k）源自共同的分布。

 數據集

分層回歸的經典示例是 鼠數據集。該縱向數據集包含5周內測得的 鼠體重。讓我們加載數據：

##   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

讓我們調查數據：

library(ggplot2)
ggplot(ddf, aes(x = variable, y = value, group = Group)) + geom_line() + geom_point()

​

 數據顯示線性增長趨勢對於不同的大鼠非常相似。但是，我們還看到，大鼠的初始體重不同，需要不同的截距，並且生長速度也需要不同的斜率。因此，分層模型似乎是適當的。

層次回歸模型的規范

該模型可以指定如下：

​

第i個大鼠的截距由αiαi表示，斜率由βiβi表示。注意，測量時間的中心是x = 22x = 22，它是時間序列數據的中值測量值（第22天）。

現在，我們可以指定模型並將其存儲在名為的文件中  rats.stan：

請注意，模型代碼估算的是方差（  sigmasq  變量）而不是標准偏差。 

資料准備

為了准備模型數據，我們首先將測量點提取為數值，然后將所有內容編碼為列表結構：

days <- as.numeric(regmatches(colnames(df), regexpr("[0-9]*$", colnames(df))))
rat.data <- list(N = nrow(df), T = ncol(df), x = days,
                 y = df, xbar = median(days)) 

擬合回歸模型

現在，我們可以為老鼠體重數據集擬合貝葉斯層次回歸模型：

rat.model <- stan(
  file = "rats.stan",
  data = rat.data)
# model contains estimates for intercepts (alpha) and slopes (beta)

層次回歸模型的預測

在確定了每只大鼠的αα和ββ之后，我們現在可以估計任意時間點單個大鼠的體重。在這里，我們有興趣尋找從第0天到第100天的大鼠體重。

 ​

ggplot(pred.df[pred.df$Rat %in% sel.rats, ], 
       aes(x = Day, y = Weight, group = Rat, 
           ymin = Lwr_Weight, ymax = Upr_Weight)) +   
    geom_line()  +
    geom_errorbar(width=0.2, size=0.5, color="blue")

​

與原始數據相比，該模型的估計是平滑的，因為每條曲線都遵循線性模型。研究最后一個圖中所示的置信區間，我們可以看到方差估計是合理的。我們對采樣時（第8至36天）的老鼠體重充滿信心，但是隨着我們遠離采樣區域，不確定性會增加。