Stan은 R이나 파이썬 같은 프로그램에서 패키지형태로 제공되는 통계 패키지입니다.
제가 사용은 했지만 정확하게 어디까지 프로그램이고 어디까지 패키지이고 개념이 아직 헷갈리니 글 읽으시는 똑똑한 분들은 제가 용어를 잘못 쓰더라도 잘 이해해 주세요 ㅠㅠ
Stan은 베이지안 모델링을 하기 위한 아주 유용한 패키지인데요. 모델의 분포만 잘 정해주면 알아서 MCMC 샘플링을 돌려서 뚝딱 결과를 내주기 때문에 초보자도 코드를 조금만 이해하면 아주 유용하게 쓸 수 있습니다.
이번 포스팅에서는 Stan을 써보고 싶은데 어떻게 쓰는지 모르는 분들을 위해 R의 rstan 패키지를 이용하는 간단한 예시를 통해 설명해보겠습니다.
우선 R에서 패키지를 이용하려면 설치하고 라이브러리로 가져와야겠죠.
install.packages("rstan")로 rstan패키지를 설치해주고,
library(rstan)로 rstan 패키지를 라이브러리로 불러왔습니다.
라이브러리로 불러오니 샘플링에 멀티코어를 활용할 수 있게 해주는 명령어인 option(mc.cores=parallel::detectCores())와 같은 stan 모델을 반복해서 컴파일링하는것을 방지해주는 rstan_options(auto_write=TRUE)를 추천하네요.
추천하는대로 바로 적용해줍니다.
이제 간단한 베이지안 모델링을 해보겠습니다. 우리가 농작물 수확량의 사전 정보를 가지고 있다고 가정하고, 새롭게 주어진 정보를 이용해 농작물 수확량 평균을 추론하는 모델링을 하겠습니다.
Prior distribution과 Likelihood 모두 Normal distribution을 따른다고 가정하겠습니다.
농작물 수확량 평균의 Prior distribution은 Normal(3,1)로 주어진 것으로 가정합니다.
농작물 수확량의 Likelihood는 Normal(mu,1)로 가정합니다.
Prior distribution과 Likelihood는 예시를 간단하기 위해 임의로 설정한 것입니다. 모델링을 어떻게 하냐에 따라 이 분포는 얼마든지 달라질 수 있습니다.
데이터는 R에 기본적으로 내장되어있는 PlantGrowth 데이터셋 중에 Control 그룹의 데이터를 사용했습니다.
이제 Prior distribution, Likelihood, 데이터가 준비되었으니 모든 준비가 끝났네요.
Stan을 이용해 모델을 작성해 보겠습니다.
Stan은 모델을 R 코드와 별도로 작성해야하는 특징이 있습니다. 물론 R코드 안에 Stan모델을 집어넣을 수도 있지만 개인적으로 다른 파일로 작성하는 것이 깔끔해보인다고 생각합니다.
위와같에 모델을 작성했는데요. Stan모델을 크게 3개의 블록으로 나뉩니다.
먼저 Data 블록. Data 블록에는 측정되는 변수들의 정보를 입력해야합니다.
블록은 { 내용 } 형태로 범위가 정해지고, 각 줄 마지막에는 ;를 붙여줘야합니다.
각주는 //각주내용 으로 사용할 수 있습니다.
측정된 데이터 자체를 입력하는것은 아니고, 데이터가 어떤 형태인지를 작성해야합니다.
예시에서는 데이터가 수확량 하나밖에 없었습니다. 그래서 두번째줄인 vector[n] y; 만 들어가면 될것 같지만, 수확량 데이터의 갯수인 n역시 정의해주어야 합니다.
Data 블록에 입력가능한 변수의 형태로는 int(정수), real(실수), vector(벡터), matrix(행렬)형태가 있고, 벡터는 vector[a]의 형태로 길이가 a라는 것을 입력하고, 행렬은 matrix[b,c] 형태로 b행, c열의 행렬이라는 것을 표현합니다.
데이터의 최대값과 최소값을 설정하고 싶다면 형태 뒤에 <lower=d, upper=e>의 형태로 최소 d, 최대 e까지 범위를 제한할 수 있습니다. 마지막에 오는건 변수이름이겠죠.
예시 코드에서는 먼저 최소값이 1이며 정수인 관측값의 갯수 n을 정의하고, 관측값데이터를 벡터 y로 정의했습니다.
다음은 Parameter 블록인데요. 여기서는 우리가 추론하고싶은 변수를 정의해 줍니다.
예시에서는 수확량의 평균을 추론하고싶기 때문에 실수값인 평균을 real mu; 형태로 정의해주었습니다.
마지막으로 Model 블록입니다. 여기서는 Prior과 Likelihood를 정의해줍니다.
예시에서는 mu의 prior를 Normal(3,1)로 정해주었기 때문에 mu ~ normal(3,1)이라고 코드를 작성했습니다. stan에서는 normal의 두번째 파라메터를 분산이 아닌 표준편차로 인식하기때문에 주의가 필요합니다.
y의 likelihood를 Normal(mu,1)로 정의해 주었기 때문에 마찬가지로 y ~ normal(mu,1)이라고 코드를 작성했습니다.
워낙 간단한 예시여서 모델도 간단하게 나왔네요. 이제 코드를 stan 코드로 저장해줘야합니다.
Rstudio 기준으로 stan코드를 저장하는 방법은 간단합니다. 파일이름 확장자를 .stan으로 설정해주면 자동으로 R 코드가 아닌 stan 코드로 저장됩니다.
이제 stan 코드 작성을 마쳤으니 수확량 평균을 추론해 봅시다.
우선 R에 기본적으로 내장된 PlantGrowth 데이터셋을 data("PlantGrowth")로 불러왔습니다.
PlantGrowth 데이터셋은 이런 모습입니다.
우리는 이중에 yield<-PlantGrowth[1:10,1]을 이용해 처음 10개의 관측값만 사용합니다.
example<-stan_model("Example.stan") 코드를 이용해 먼저 만들어둔 stan코드를 example라는 이름으로 불러옵니다.
stan 모델의 데이터는 list형태로 작성되어야 하며, stan 코드에서 정의된 이름과 형태로 주어져야합니다.
data<-list(n=10,y=yield) 는 관측값의 갯수 n은 10으로, 관측값 y는 yield로 불러왔습니다.
이제 sample<-sampling(example,data)을 이용해 MCMC 샘플링을 해줍니다. sampling에서 첫번째 인자는 stan모델, 두번째 인자는 데이터가 저장된 list를 입력해줘야 합니다.
코드를 실행하면 위 그림과 같이 샘플링이 진행됩니다. 기본적으로 4개의 체인으로 2000번씩 샘플링을 진행합니다. 그리고 각 체인마다 처음 1000개의 샘플을 버려줍니다.
마지막으로 샘플링의 내용을 저장한 sample를 입력해봤습니다.
우리가 추론하고싶은 mu의 평균은 4.98 평균의 표준편차는 0.01로 나왔습니다. 표준편차가 아주 작아 평균을 추론값으로 사용해도 문제가 없겠네요.
추가로 각 분위수를 확인할 수 있고, n_eff는 유효한 샘플의 갯수입니다.
4개 체인이 각각 2000개의 샘플중 처음 1000개의 샘플을 버리고 나머지 1000개의 샘플을 남겼기 때문에 총 4000개의 샘플이 있는데 그중에 1555개의 샘플이 유효하다는 의미입니다.
샘플링 중 앞뒤 샘플간의 상관관계나 여러가지 이유로 유효샘플은 보통 추출한 샘플의 양보다 적습니다.
Rhat은 모델이 잘 수렴했는지를 알려주는 지표입니다. Rhat이 1에 가까울수록 잘 수렴했다는 뜻이 됩니다. 보통 1에 굉장히 가깝게 나오기 때문에 1.2 이하가 나오지 않는다면 모델에 문제가 있다는 뜻이 됩니다.
간단하게 rstan을 이용해 베이지안 추론을 하는 법을 알아봤습니다. 공부하시는데 작은 도움이 되는 포스팅이 됬으면 좋겠네요. 읽어주셔서 감사합니다.
'통계학' 카테고리의 다른 글
| [Latex] 이미지에 주석 달기 (2) | 2022.06.06 |
|---|---|
| [Latex]모르는 기호가 나왔을때 어떻게 찾죠? (0) | 2022.05.31 |
| [Latex]Latex 문서 작성할 때 유용한 사이트들 (표, 수식, 기호 도우미들) (1) | 2021.08.07 |
| [Latex]글씨 두껍게 이탤릭 밑줄 등 (0) | 2021.02.10 |