RStan Getting Started (한국어)

RStan은 Stan의 R 인터페이스이다. Stan과 Stan을 사용한 모델링 언어에 관한 자세한 정보는 Stan 웹사이트에서 볼 수 있다.

mc-stan.org

최신 버전: 2.17.3 (2018년 1월 20일)

아래에 기술된 대부분의 설치 방법은 앞에서 언급한 RStan 버전에 적용된다.

설치

사용자의 OS에 맞는 RStan 설치 방법은 다음 링크를 참조 바란다.

자신에게 맞는 설치 방법을 다 따라하고도 제대로 설치가 되지 않는다면, Stan 사용자 메일링 리스트에서 도움을 받을 수 있다.

Stan 대화 포럼(영문)

RStan 사용법

이 문서의 나머지 부분은 사용자가 위 링크를 따라 RStan을 이미 설치했다고 가정하였다.

패키지 불러오기

패키지 이름은 rstan(모두 소문자)이고 library("rstan") 명령어로 패키지를 불러올 수 있다.

library("rstan") # 시작 메세지를 눈여겨보길 바란다

만약 사용자의 모형을 rstan을 이용해 멀티코어와 RAM이 많은 시스템에서 병렬로 계산하려면, 시작 메세지가 말해주는 지금 다음 명령어를 입력한다.

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

두 옵션은 자동으로 컴파일된 bare 버전 Stan 프로그램을 하드디스트에 저장하여 다중 말코프 체인을 병렬로 실행할 때마다 컴파일하지 않게 해준다.

예제 1: 여덟 개의 학교

이 예제는 겔만과 연구자들(2003)의 책 섹션 5.5 실려있는 예제로, 여덟 개의 학교에서 코칭 효과를 연구한 것이다. 단순하게 이 예제를 "여덟 개의 학교”라 부르자.

먼저 8schools.stan로 새 파일을 만들고 모델을 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의 변형된 매개변수로 정의한다. 이렇게 하므로써 표본을 더 효과적으로 추출할 수 있다. (자세한 설명(영문)). 데이터는 R에서 다음과 같이 준비한다.

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 명령어로 fit을 구할 수 있다. file = 의 인수는 R의 작업 디렉토리에 파일이 저장돼있지 않다면, 자신의 파일 시스템의 파일 위치로 지정해야한다.

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

Stan 모형을 지정할 때 stan 함수의 model_code 인수를 사용해 문자열로 지정할 수도 있다. 하지만 파일을 사용하면 print 문을 사용할 수 있고 오류 메시지의 줄 번호를 사용할 수 있으므로 문자열로 지정하는 것은 권장하지 않는다.

stan 함수에서 반환된 fit 객체는 stanfit 클래스의 S4 객체이다. print, plot, pairs와 같은 메소드는 fit 결과와 관련되어 있으므로 앞으로 나올 코드를 사용해 fit안의 결과를 확인할 수 있다. print는 모델의 매개변수 요약 뿐만 아니라 log-사후분포도 lp__라는 이름으로 보여준다 (다음 출력 예를 참고). 더 많은 메소드와 stanfit 클래스의 자세한 내용을 stanfit 클래스의 도움말을 읽어보길 바란다.

특히, 표본을 추출하기 위해 stanfit 객체의 extract 함수를 사용할 수 있다. extract 함수는 stanfit 객체에서 원하는 매개변수의 표본을 배열이나 배열 리스트로 추출한다. 참고로 stanfit 객체에 사용할 수 있는 S3 함수인 as.array, as.matrix, as.data.frame도 정의되어 있다(R에서 help("as.array.stanfit") 명령어로 도움말을 볼 수 있다).

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

la <- extract(fit, permuted = TRUE) # 배열 리스트를 반환한다 
mu <- la$mu 

### iterations, chains, parameters의 3차원 배열을 반환한다
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

2017년 5월 5일 금요일 10:41:43에 NUTS(diag_e)를 사용하여 표본을 추출하였다.
각 매개변수에 대해 n_eff는 유효표본크기의 대략적인 척도이고 Rhat은 나뉜 체인(수렴하고, Rhat=1일때)의
potential scale reduction factor이다.

덧붙여, BUGS (또는 JAGS)와 마찬가지로 CmdStan(Stan의 커맨드라인 인터페이스)를 사용하려면 R 덤프 파일의 모든 데이터가 필요하다. 이 파일을 가지고 있다면 rstan의 함수 read_rdump를 사용해 R 리스트로 모든 데이터를 읽을 수 있다. 예를 들어 다음 내용을 가진 "8schools.rdump" 파일을 작업 디렉토리에 가지고 있다고 하자. In addition, as in BUGS (or JAGS), CmdStan (the command line interface to Stan) needs all the data to be in an R dump file. In the case we have this file, rstan provides

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')

사실 R의 source 함수를 사용해 R 덤프 파일을 전역 환경 변수로 넣을 수 있다. 이때는 data 인수를 생략할 수 있고 stan은 8schools.stan 블록 데이터 안의 같은 이름을 가진 객체를 환경변수에서 검색한다.

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

예제 2: 쥐

쥐 예제도 유명한 예제이다. 예를 들어 OpenBUGS 버전(Gelfand와 연구자들(1990)에서 나옴). 데이터는 5주동안의 쥐 30마리의 성장을 매 주 관찰한 것이다. 다음 표는 데이터 리스트이고 x는 데이터를 수집한 날짜를 나타낸다. rats.txt 데이터와 모형 코드 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')

예제 3: 아무거나

많은 BUGS 예제와 Stan에 우리자 작성한 예제를 다음 명령어로 실행하여 팝업 메뉴에서 모델을 선택해 실행해 볼 수 있다.

model <- stan_demo()

처음 stan_demo()를 호출하면 다음 예제를 내려받을 것인지를 물어본다. rstan이 설치된 디렉토리에 내려받아 다음에 또다시 다운받지 않도록 하려면 옵션 1을 선택한다. 위의 model 객체는 stanfit 클래스의 인스턴스이므로 print, plot, pairs, extract 등을 사용할 수 있다.

더 많은 도움말

RStan의 더 자세한 설명은 rstan의 문서에서 찾을 수 있다. 예를 들어 help(stan)과 help("stanfit-class")를 사용해서 stan과 S4 클래스 stanfit 함수의 도움말을 볼 수 있다. 또한, Stan 모형화 언어 설명서(영문)를 참조하여 Stan 표본 추출기, 최적화기 그리고 Stan 모형화 언어에 대한 자세한 내용을 볼 수 있다.

덧붙여 Stan 사용자 메일링 리스트(영문) 에서도 Stan 사용법에 관해 이야기 할 수 있으모 예제를 올리거나 (R)Stan에 관해 질문할 수 있다. 도움이 필요할때는 다음 정보를 충분하게 함께 올려야 한다:

Stan 모형화 언어로 작성된 Stan 코드
데이터
필요한 R 코드
verbose=TRUE와 cores=1을 사용했을 때 나오는 오류 메시지 덤프
g++ -v를 입력해 출력되는 사용된 C++ 컴파일러 버전
sessionInfo() 함수를 사용해 출력되는 사용한 R에 관한 정보

참고 문헌

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

Provide feedback

Saved searches

Use saved searches to filter your results more quickly