Skip to content

RStan Getting Started (繁體中文)

Jump to bottom
Andrew Johnson edited this page Feb 24, 2023 · 4 revisions

翻譯:謝男鴻(Nan-Hung Hsieh)


RStan 是以 R 作為操作界面以此來執行Stan。更多有關於Stan的資訊以及模擬用的程式語言都可參考Stan的官方網頁。

目前最新的Stan版本為： 2.17.3   (2018年1月20日)

以下所描述的安裝指引皆可以適用於較早的RStan版本。

安裝

首先，根據您所使用的作業系統並參考下列連結安裝Rstan：

假如您已經依上述連結安裝但卻無法安裝成功，您可以到Stan的論壇發表討論您所遇到問題：

如何使用RStan

假如您已經依照上述的安裝步驟成功的安裝RStan，接著請繼續參考以下的說明內容：

Loading the package

RStan套件的名稱為rstan (皆為小寫)，所以我們可以在R控制台中以library("rstan")來進行讀取

library("rstan") # 讀取後可以看到相關的啟動訊息

如啟動訊息所描述，假如您所使用的電腦系統具有多個中央處理器（CPU）並且有充足的記憶體（RAM）來進行平行運算模擬，則建議使用下列的設定步驟：

rstan_options(auto_write = TRUE)
options(mc.cores = parallel::detectCores())

這兩個步驟分別讓您可以自動儲存編譯後的Stan執行程式到硬碟中，以此可以免去重複的編譯程序，並且可以平行的運算執行多個馬可夫鍊。

範例一：八所學校（Eight Schools） - 教學成果分析

這是在Gelman等人於2003年的Bayesian Data Analysis這本書第五章第五節所提供的一個簡單範例，探討八所學校在不同的教學方式上對學術水準測驗考試(SAT, Scholastic Aptitude Test)所產生的不同影響。 在此，我們將此範例簡稱為八所學校（Eight Schools）。首先編寫一Stan程式的模擬模型並儲存及命名為8schools.stan：

// 存檔為 8schools.stan
data {
  int<lower=0> J; // 學校數目 
  real y[J]; // 測試效果的預測值
  real<lower=0> sigma[J]; // 測試效果的標準差 
}
parameters {
  real mu; 
  real<lower=0> tau;
  real eta[J];
}
transformed parameters {
  real theta[J];
  for (j in 1:J)
    theta[j] = mu + tau * eta[j];
}
model {
  target += normal_lpdf(eta | 0, 1);
  target += normal_lpdf(y | theta, sigma);
}

在此模式中，我們讓theta成為mu 及 eta的轉換參數，而不是直接的定義theta為一模式參數。 透過這樣的參數化過程可以更有效率的進行抽樣(細節解釋可參考這網頁說明)。接著我們可以定義所要探討的研究資料：

schools_dat <- list(J = 8, 
                    y = c(28,  8, -3,  7, -1,  1, 18, 12),
                    sigma = c(15, 10, 16, 11,  9, 11, 10, 18))

之後我可以藉由以下的R指令來得到擬合分析結果。值得注意的是，file = 必須能確實指示出8schools.stan這個檔案的所在位置，除非檔案已經確實放置在目前的R工作目錄中（在此條件下便可以執行下列的指令）。

fit <- stan(file = '8schools.stan', data = schools_dat, 
            iter = 1000, chains = 4)

另外，我們也可以使用stan內建的model_code功能來使用字元串來設定Stan模型。但我們建議使用上述的方法，因為使用定義的檔案在file = 這個功能可以讓您輸出狀態描述（print statement）並反參考(dereference)其中的錯誤訊息。

透過stan函式所獲得的輸出結果fit即為stanfit中的S4項目 。其他方法如print、plot及pairs皆是與擬合結果有關，因此我們可以透先前儲存的fit來檢試結果。 print提供模式參數以及對數化的後驗值（log-posterior）lp__的推估結果總覽（請參考以下的輸出範例）。關於stanfit的更多操作方法及細節可以進一步參考stanfit的說明頁面。

另外，我們可以使用extract這個函式對stanfit的項目來檢視分析內容的詳細結果。extract擷取了所有stanfit項目中的抽樣結果，並條列(list)的以陣列(array)方式列出參數的抽樣值。

另外，S3的函式庫如as.array、as.matrix及as.data.frame可進一步定義stanfit的項目為陣列、矩陣及數據表（建議透過help("as.array.stanfit")檢視R的輔助文件）。

print(fit)
plot(fit)
pairs(fit, pars = c("mu", "tau", "lp__"))

la <- extract(fit, permuted = TRUE) # 回傳陣列列表 
mu <- la$mu 

### 回傳三維的陣列：疊代、馬可夫鍊及模式參數
a <- extract(fit, permuted = FALSE) 

### 對stanfit項目使用S3函式
a2 <- as.array(fit)
m <- as.matrix(fit)
d <- as.data.frame(fit)
> print(fit, digits = 1)
Inference for Stan model: 8schools.
4 chains, each with iter=1000; warmup=500; thin=1; 
post-warmup draws per chain=500, total post-warmup draws=2000.

          mean se_mean  sd  2.5%   25%   50%   75% 97.5% n_eff Rhat
mu         8.2     0.2 5.4  -1.9   4.8   8.1  11.3  19.3   480    1
tau        6.8     0.3 6.2   0.3   2.5   5.2   9.2  21.7   425    1
eta[1]     0.4     0.0 1.0  -1.5  -0.3   0.4   1.0   2.2  2000    1
eta[2]     0.0     0.0 0.8  -1.7  -0.6   0.0   0.5   1.7  2000    1
eta[3]    -0.2     0.0 1.0  -2.1  -0.9  -0.2   0.4   1.7  2000    1
eta[4]    -0.1     0.0 0.9  -1.8  -0.7  -0.1   0.5   1.7  2000    1
eta[5]    -0.4     0.0 0.9  -2.1  -1.0  -0.4   0.2   1.4  2000    1
eta[6]    -0.2     0.0 0.9  -1.9  -0.8  -0.2   0.4   1.5  1731    1
eta[7]     0.3     0.0 0.9  -1.4  -0.2   0.4   0.9   2.0  1507    1
eta[8]     0.0     0.0 0.9  -1.9  -0.6   0.0   0.7   1.8  1988    1
theta[1]  11.5     0.3 8.8  -2.4   5.9  10.1  15.6  32.9   977    1
theta[2]   7.8     0.1 6.2  -4.7   4.1   7.9  11.6  20.3  2000    1
theta[3]   6.1     0.2 7.7 -11.2   2.1   6.4  10.5  20.2  2000    1
theta[4]   7.6     0.1 6.5  -4.9   3.8   7.8  11.4  21.3  2000    1
theta[5]   5.0     0.1 6.6  -9.3   1.2   5.6   9.3  16.7  2000    1
theta[6]   6.2     0.2 6.7  -8.2   2.2   6.5  10.5  18.5  2000    1
theta[7]  10.8     0.2 7.0  -1.3   6.1  10.1  15.1  26.8  2000    1
theta[8]   8.7     0.2 8.2  -7.3   3.9   8.4  12.8  27.2  1446    1
lp__     -39.5     0.1 2.6 -45.1 -41.2 -39.4 -37.7 -35.1   590    1

Samples were drawn using NUTS(diag_e) at Fri May  5 10:41:43 2017.
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).

另外，如同於BUGS（或JAGS），在使用CmdStan（Stan的命令列界面）進行資料讀取時，必須先將所有的模式及參數資料儲存於以rdump為副檔名的檔案內。 假如我們已經建立了rdump檔案，我們就可以藉由read_rdump這個函式讀取所有資料到R列表中。在此範例中，假如我們擁有一個"8schools.rdump"的檔案，裡頭包括下列下列已經建立在工作目錄下的資料變數：

J <- 8
y <- c(28,  8, -3,  7, -1,  1, 18, 12)
sigma <- c(15, 10, 16, 11,  9, 11, 10, 18)

接著我們便可以用下列方式來讀取"8schools.rdump"中的資料：

schools_dat <- read_rdump('8schools.rdump')

實際上，透過source這個R內建函式也可以對rdump中的資料進行讀取。以下是透過著個方式，可以簡易略過對輸入資料data的宣告，進而直接讓stan搜尋在作業目錄下與8schools.stan的資料區具有相同命名的檔案：

source('8schools.rdump') 
fit <- stan(file = '8schools.stan', iter = 1000, chains = 4)

範例二: 鼠群 （Rats） - 成長分析

鼠群也是一個著名的範例。這個例子同樣也在OpenBUGS的範例中使用過，其主要是Gelfand等人於1990年所發表的研究，針對30隻大鼠於五週內的成長紀錄。 以下的表格表示了老鼠的代號，x表示資料紀錄時的觀測天數。我們可以由連結內的資料rats.txt及Stan模式原始碼rats.stan來進行範例演練。

Rat x=8 x=15 x=22 x=29 x=36 Rat x=8 x=15 x=22 x=29 x=36
1 151 199 246 283 320 16 160 207 248 288 324
2 145 199 249 293 354 17 142 187 234 280 316
3 147 214 263 312 328 18 156 203 243 283 317
4 155 200 237 272 297 19 157 212 259 307 336
5 135 188 230 280 323 20 152 203 246 286 321
6 159 210 252 298 331 21 154 205 253 298 334
7 141 189 231 275 305 22 139 190 225 267 302
8 159 201 248 297 338 23 146 191 229 272 302
9 177 236 285 350 376 24 157 211 250 285 323
10 134 182 220 260 296 25 132 185 237 286 331
11 160 208 261 313 352 26 160 207 257 303 345
12 143 188 220 273 314 27 169 216 261 295 333
13 154 200 244 289 325 28 157 205 248 289 316
14 171 221 270 326 358 29 137 180 219 258 291
15 163 216 242 281 312 30 153 200 244 286 324 
y <- as.matrix(read.table('https://raw.github.com/wiki/stan-dev/rstan/rats.txt', header = TRUE))
x <- c(8, 15, 22, 29, 36)
xbar <- mean(x)
N <- nrow(y)
T <- ncol(y)
rats_fit <- stan(file = 'https://raw.githubusercontent.com/stan-dev/example-models/master/bugs_examples/vol1/rats/rats.stan')

範例三: 所有範例

您可以執行以下由Stan所內建的函式來演練所有的BUGS範例

model <- stan_demo()

從彈出的列表清單中選擇一個範例模型來練習。在第一次使用stan_demo()時，它會先詢問您是否要下載這些範例。您必須選擇選第一項來將這些範例儲存在rstan所安裝的目錄下，因此這些範例檔案在之後就可以直接讀取使用而不須要重複下載。 前面所描述的model會儲存所選擇的stanfit參考範例，模擬之後就可以接著使用print、plot、pairs以及extract等等函式來進行資料分析。

更多輔助內容

更多關於RStan的細節資料可以參考在rstan套件中的示範（vignette）文件。例如，使用help(stan)及help("stanfit-class")來檢視stan函式以及stanfit這個分類中的所有介紹。 同樣的，也可以參考 Stan's modeling language manual關於Stan抽樣、最佳化及Stan程式語言的詳細內容。 另外，Stan的社群論壇Stan User's Mailing list可以用來討論Stan的使用、張貼範例以及關於(R)Stan的問題詢問。 當您有需要時，建議能提供充足的資訊如下：

  • Stan模式的原始碼
  • 研究資料（data）
  • 必要的R原始碼
  • 當使用stan時所出現的所有錯誤訊息，可以由verbose=TRUE及cores=1來顯示
  • c++編譯器的版本。例如，使用g++ -v來測試gcc是否能正常使用
  • 關於所使用R的相關資訊，可以在R中執行sessionInfo()來檢視

參考資料

  • Gelman, A., Carlin, J. B., Stern, H. S., and Rubin, D. B. (2003). Bayesian Data Analysis, CRC Press, London, 2nd Edition.
  • Stan Development Team. Stan Modeling Language User's Guide and Reference Manual.
  • Gelfand, A. E., Hills S. E., Racine-Poon, A., and Smith A. F. M. (1990). "Illustration of Bayesian Inference in Normal Data Models Using Gibbs Sampling", Journal of the American Statistical Association, 85, 972-985.
  • Stan
  • R
  • BUGS
  • OpenBUGS
  • JAGS
  • Rcpp
Clone this wiki locally