Pengenalan STAN di RStudio
Copy Right @ Prof Dr Suhartono M.Kom @ Magister Informatika UIN Maulana Malik Ibrahim@ (http://suhartono.lecturer.uin-malang.ac.id/riwayat-hidup/)
18 Januari 2021
library(rstan)## Loading required package: StanHeaders## Loading required package: ggplot2## rstan (Version 2.21.2, GitRev: 2e1f913d3ca3)## For execution on a local, multicore CPU with excess RAM we recommend calling
## options(mc.cores = parallel::detectCores()).
## To avoid recompilation of unchanged Stan programs, we recommend calling
## rstan_options(auto_write = TRUE)Stan adalah bahasa pemrograman probabilistik untuk menentukan model statistik. Stan memberikan inferensi Bayesian lengkap untuk model variabel kontinu melalui metode Markov Chain Monte Carlo seperti sampler No-U-Turn, bentuk adaptif dari pengambilan sampel Hamiltonian Monte Carlo. Perkiraan kemungkinan maksimum yang dikenai sanksi dihitung menggunakan metode pengoptimalan seperti algoritma Broyden-Fletcher-Goldfarb-Shanno memori terbatas. Stan dapat dipanggil melalui R menggunakan rstan paket, dan melalui Python menggunakan pystan paket. Kedua antarmuka mendukung pengambilan sampel dan inferensi berbasis pengoptimalan dengan diagnostik dan analisis posterior. Dalam pembicaraan ini diperlihatkan sekilas tentang sifat-sifat utama Stan. Juga ditunjukkan beberapa contoh: yang pertama terkait dengan model Bernoulli sederhana dan yang kedua, tentang model Lotka-Volterra berdasarkan persamaan diferensial biasa.
// saved as schools.stan
data { int<lower=0> J; // number of schools real y[J]; // estimated treatment effects real<lower=0> sigma[J]; // standard error of effect estimates } parameters { real mu; // population treatment effect real<lower=0> tau; // standard deviation in treatment effects vector[J] eta; // unscaled deviation from mu by school } transformed parameters { vector[J] theta = mu + tau * eta; // school treatment effects } model { target += normal_lpdf(eta | 0, 1); // prior log-density target += normal_lpdf(y | theta, sigma); // log-likelihood }
Pastikan program Stan Anda diakhiri dengan baris kosong tanpa karakter apa pun termasuk spasi dan komentar.
Dalam program Stan ini, kita membiarkan theta menjadi transformasi mu, eta, dan tau daripada mendeklarasikan theta di blok parameter, yang memungkinkan sampler akan berjalan lebih efisien (lihat penjelasan detail). Dalam kasus ini, kita akan mensimulasikan sampel acak untuk digunakan dalam contoh percobaan ini.
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))Untuk melihat data dapat dilakukan sebagai berikut :
print(schools_dat)## $J
## [1] 8
##
## $y
## [1] 28 8 -3 7 -1 1 18 12
##
## $sigma
## [1] 15 10 16 11 9 11 10 18Maka perintah stan dapat dilakukan dengan argumen ke file pada tempat file berada pada sistem file.
fit <- stan(file = 'schools.stan', data = schools_dat)## Trying to compile a simple C file## Running /usr/lib/R/bin/R CMD SHLIB foo.c
## gcc -std=gnu99 -I"/usr/share/R/include" -DNDEBUG -I"/home/suhartono/R/x86_64-pc-linux-gnu-library/3.6/Rcpp/include/" -I"/home/suhartono/R/x86_64-pc-linux-gnu-library/3.6/RcppEigen/include/" -I"/home/suhartono/R/x86_64-pc-linux-gnu-library/3.6/RcppEigen/include/unsupported" -I"/home/suhartono/R/x86_64-pc-linux-gnu-library/3.6/BH/include" -I"/home/suhartono/R/x86_64-pc-linux-gnu-library/3.6/StanHeaders/include/src/" -I"/home/suhartono/R/x86_64-pc-linux-gnu-library/3.6/StanHeaders/include/" -I"/home/suhartono/R/x86_64-pc-linux-gnu-library/3.6/RcppParallel/include/" -I"/home/suhartono/R/x86_64-pc-linux-gnu-library/3.6/rstan/include" -DEIGEN_NO_DEBUG -DBOOST_DISABLE_ASSERTS -DBOOST_PENDING_INTEGER_LOG2_HPP -DSTAN_THREADS -DBOOST_NO_AUTO_PTR -include '/home/suhartono/R/x86_64-pc-linux-gnu-library/3.6/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp' -D_REENTRANT -DRCPP_PARALLEL_USE_TBB=1 -fpic -g -O2 -fdebug-prefix-map=/build/r-base-V28x5H/r-base-3.6.3=. -fstack-protector-strong -Wformat -Werror=format-security -Wdate-time -D_FORTIFY_SOURCE=2 -g -c foo.c -o foo.o
## In file included from /home/suhartono/R/x86_64-pc-linux-gnu-library/3.6/RcppEigen/include/Eigen/Core:88:0,
## from /home/suhartono/R/x86_64-pc-linux-gnu-library/3.6/RcppEigen/include/Eigen/Dense:1,
## from /home/suhartono/R/x86_64-pc-linux-gnu-library/3.6/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp:13,
## from <command-line>:0:
## /home/suhartono/R/x86_64-pc-linux-gnu-library/3.6/RcppEigen/include/Eigen/src/Core/util/Macros.h:628:1: error: unknown type name ‘namespace’
## namespace Eigen {
## ^~~~~~~~~
## /home/suhartono/R/x86_64-pc-linux-gnu-library/3.6/RcppEigen/include/Eigen/src/Core/util/Macros.h:628:17: error: expected ‘=’, ‘,’, ‘;’, ‘asm’ or ‘__attribute__’ before ‘{’ token
## namespace Eigen {
## ^
## In file included from /home/suhartono/R/x86_64-pc-linux-gnu-library/3.6/RcppEigen/include/Eigen/Dense:1:0,
## from /home/suhartono/R/x86_64-pc-linux-gnu-library/3.6/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp:13,
## from <command-line>:0:
## /home/suhartono/R/x86_64-pc-linux-gnu-library/3.6/RcppEigen/include/Eigen/Core:96:10: fatal error: complex: No such file or directory
## #include <complex>
## ^~~~~~~~~
## compilation terminated.
## /usr/lib/R/etc/Makeconf:168: recipe for target 'foo.o' failed
## make: *** [foo.o] Error 1
##
## SAMPLING FOR MODEL 'schools' NOW (CHAIN 1).
## Chain 1:
## Chain 1: Gradient evaluation took 1.2e-05 seconds
## Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 0.12 seconds.
## Chain 1: Adjust your expectations accordingly!
## Chain 1:
## Chain 1:
## Chain 1: Iteration: 1 / 2000 [ 0%] (Warmup)
## Chain 1: Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 1: Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 1: Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 1: Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 1: Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 1: Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 1: Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 1: Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 1: Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 1: Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 1: Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 1:
## Chain 1: Elapsed Time: 0.06398 seconds (Warm-up)
## Chain 1: 0.069132 seconds (Sampling)
## Chain 1: 0.133112 seconds (Total)
## Chain 1:
##
## SAMPLING FOR MODEL 'schools' NOW (CHAIN 2).
## Chain 2:
## Chain 2: Gradient evaluation took 9e-06 seconds
## Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 0.09 seconds.
## Chain 2: Adjust your expectations accordingly!
## Chain 2:
## Chain 2:
## Chain 2: Iteration: 1 / 2000 [ 0%] (Warmup)
## Chain 2: Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 2: Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 2: Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 2: Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 2: Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 2: Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 2: Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 2: Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 2: Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 2: Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 2: Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 2:
## Chain 2: Elapsed Time: 0.066281 seconds (Warm-up)
## Chain 2: 0.084312 seconds (Sampling)
## Chain 2: 0.150593 seconds (Total)
## Chain 2:
##
## SAMPLING FOR MODEL 'schools' NOW (CHAIN 3).
## Chain 3:
## Chain 3: Gradient evaluation took 1.1e-05 seconds
## Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 0.11 seconds.
## Chain 3: Adjust your expectations accordingly!
## Chain 3:
## Chain 3:
## Chain 3: Iteration: 1 / 2000 [ 0%] (Warmup)
## Chain 3: Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 3: Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 3: Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 3: Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 3: Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 3: Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 3: Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 3: Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 3: Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 3: Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 3: Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 3:
## Chain 3: Elapsed Time: 0.068524 seconds (Warm-up)
## Chain 3: 0.05686 seconds (Sampling)
## Chain 3: 0.125384 seconds (Total)
## Chain 3:
##
## SAMPLING FOR MODEL 'schools' NOW (CHAIN 4).
## Chain 4:
## Chain 4: Gradient evaluation took 9e-06 seconds
## Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 0.09 seconds.
## Chain 4: Adjust your expectations accordingly!
## Chain 4:
## Chain 4:
## Chain 4: Iteration: 1 / 2000 [ 0%] (Warmup)
## Chain 4: Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 4: Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 4: Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 4: Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 4: Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 4: Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 4: Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 4: Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 4: Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 4: Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 4: Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 4:
## Chain 4: Elapsed Time: 0.075772 seconds (Warm-up)
## Chain 4: 0.073958 seconds (Sampling)
## Chain 4: 0.14973 seconds (Total)
## Chain 4:print(fit)## 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.77 0.12 5.06 -2.10 4.46 7.79 11.00 18.03 1646 1
## tau 6.76 0.14 5.54 0.30 2.58 5.44 9.40 21.65 1528 1
## eta[1] 0.41 0.02 0.95 -1.48 -0.20 0.43 1.05 2.28 3690 1
## eta[2] 0.01 0.01 0.85 -1.69 -0.55 0.01 0.56 1.69 3684 1
## eta[3] -0.20 0.02 0.94 -2.01 -0.83 -0.21 0.43 1.68 3508 1
## eta[4] -0.03 0.01 0.86 -1.72 -0.61 -0.02 0.54 1.68 3720 1
## eta[5] -0.36 0.01 0.85 -2.00 -0.91 -0.36 0.20 1.35 3629 1
## eta[6] -0.19 0.01 0.88 -1.94 -0.78 -0.18 0.40 1.58 3542 1
## eta[7] 0.36 0.01 0.88 -1.42 -0.21 0.38 0.93 2.08 3416 1
## eta[8] 0.06 0.02 0.92 -1.79 -0.56 0.07 0.68 1.89 3140 1
## theta[1] 11.59 0.18 8.61 -2.04 5.93 10.32 15.79 32.94 2406 1
## theta[2] 7.97 0.10 6.31 -4.47 3.96 7.94 11.93 21.03 4249 1
## theta[3] 5.92 0.14 8.13 -12.76 1.78 6.35 10.79 21.18 3480 1
## theta[4] 7.53 0.10 6.37 -5.63 3.61 7.55 11.50 20.44 4463 1
## theta[5] 4.87 0.10 6.31 -9.09 1.14 5.32 9.13 15.98 3984 1
## theta[6] 6.19 0.12 6.90 -8.46 2.14 6.54 10.70 19.27 3570 1
## theta[7] 10.78 0.12 6.82 -1.13 6.21 10.24 14.58 25.66 3496 1
## theta[8] 8.32 0.14 7.88 -8.02 3.85 8.14 12.66 25.47 3255 1
## lp__ -39.40 0.07 2.56 -44.94 -40.98 -39.13 -37.56 -35.17 1360 1
##
## Samples were drawn using NUTS(diag_e) at Mon Jan 18 05:47:49 2021.
## 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).plot(fit)## 'pars' not specified. Showing first 10 parameters by default.## ci_level: 0.8 (80% intervals)## outer_level: 0.95 (95% intervals)
pairs(fit, pars = c("mu", "tau", "lp__"))
la <- extract(fit, permuted = TRUE) # return a list of arrays
mu <- la$mu
### return an array of three dimensions: iterations, chains, parameters
a <- extract(fit, permuted = FALSE)
### use S3 functions on stanfit objects
a2 <- as.array(fit)
m <- as.matrix(fit)
d <- as.data.frame(fit)Daftar pustaka :