Rstan Tutorial
kosugitti
2015年9月4日
Stanを始めよう
環境の準備の準備
用意するもの。
- パソコン
- インターネット環境
環境の準備
統計パッケージRと，Rを使いやすくしてくれる統合環境RStudioを用意しよう。
Rの準備
Rのサイトはこちらから。
- 左のDownloadから，CRANを選び，最も近いミラーサイト(普通は日本にするよ)を選択。
- 自分のOSにあった所をクリックし，Rをインストールする。
- 最新バージョンはR3.2.2です。(2015.09.04段階)
- Windowsの人はbaseとrtoolsの両方をダウンロードする必要があります。
- Rtoolsを入れるWindowsユーザはRtools33を選びましょう。(Rのバージョンに合わせる)
- インストールする順番を間違えないで！RをいれてからRtoolsを入れましょう。
- Rtoolsのインストールの途中に出てくる，「edit the system path」は必ずチェックすること。
- MacユーザはXcodeをAppStoreからダウンロード・インストールしてください。
- Ubuntu系は自分で頑張れ。
Rstudioの準備
Rstudioはこちらから。
- OpenSourceEditionのDesktop版を選び，自分のOSにあったものをダウンロードする。
三つ全てインストールが終わったら，RStudioを起動してみよう。
ついでに
Rはパッケージを取り込むことで機能が増えていきます。 今後のことを考えて心理学関係のパッケージとベイズ関係のパッケージを一気に取り込んじゃおう！
install.packages("ctv") #install.packagesはパッケージをインストールする関数。
require(ctv) #requireあるいはlibraryはパッケージを装備する関数。
install.views("Psychometrics") #ctvパッケージにある，関係パッケージをまとめて入れる関数
install.views("Bayes") #Bayes関係のパッケージをじゃんじゃん入れてしまう！rstanのインストール
実は準備が整っていたら，ここで何の苦労もない。(バージョン2.7からの恩恵！)
install.packages("rstan")とっても簡単だね。
今更ながら，Stanプロジェクトのサイトを紹介しておきます。ここのInterfacesにRstan(stanのRインターフェイス)があります。
shiny-stanのインストール
もうひとつ。stanの結果を見やすくするためのアプリ，shinyStanというのも入れておきます。 入れるためのコードは次の通り。
install.packages("shinystan",repos='http://cran.rstudio.com')rstanのサンプルコード
インストールが終わったら，Rstanの「事始め」ページにある，次のコードを実行してみよう。
require(rstan)## Loading required package: rstan
## Loading required package: Rcpp
## Loading required package: inline
##
## Attaching package: 'inline'
##
## 以下のオブジェクトは 'package:Rcpp' からマスクされています:
##
## registerPlugin
##
## rstan (Version 2.7.0-1, packaged: 2015-07-17 18:12:01 UTC, GitRev: 05c3d0058b6a)
## For execution on a local, multicore CPU with excess RAM we recommend calling
## rstan_options(auto_write = TRUE)
## options(mc.cores = parallel::detectCores())rstan_options(auto_write = TRUE)
options(mc.cores = parallel::detectCores())これは何も起こらないと思いますが，「できるだけ早くしてね」という意味です。 実際に動かすのは次のコード。
stancode <- '
data {
int<lower=0> J; // number of schools
real y[J]; // estimated treatment effects
real<lower=0> sigma[J]; // s.e. of effect estimates
}
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 {
eta ~ normal(0, 1);
y ~ normal(theta, sigma);
}
'
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))
fit <- stan(model_code = stancode, data = schools_dat,
iter = 1000, chains = 4)## COMPILING THE C++ CODE FOR MODEL '19c33eca0c63b510805d5af1962613a7' NOW.print(fit)## Inference for Stan model: 19c33eca0c63b510805d5af1962613a7.
## 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 7.96 0.28 5.78 -4.48 4.75 7.86 11.18 20.22 434 1.00
## tau 7.14 0.42 6.31 0.20 2.54 5.56 10.07 22.53 222 1.01
## eta[1] 0.42 0.03 0.91 -1.38 -0.20 0.43 1.04 2.20 1087 1.00
## eta[2] -0.02 0.03 0.86 -1.80 -0.58 0.01 0.56 1.70 1162 1.00
## eta[3] -0.20 0.03 0.96 -2.09 -0.84 -0.19 0.43 1.70 1173 1.00
## eta[4] -0.04 0.03 0.90 -1.86 -0.62 0.00 0.55 1.70 1164 1.00
## eta[5] -0.34 0.03 0.87 -2.00 -0.89 -0.35 0.21 1.46 1157 1.00
## eta[6] -0.21 0.03 0.86 -1.91 -0.74 -0.20 0.35 1.54 1176 1.00
## eta[7] 0.29 0.03 0.87 -1.36 -0.33 0.32 0.88 1.93 1060 1.00
## eta[8] 0.06 0.02 0.90 -1.71 -0.55 0.07 0.67 1.76 1378 1.00
## theta[1] 11.92 0.32 8.52 -0.77 6.22 10.17 16.37 32.86 718 1.00
## theta[2] 8.08 0.16 6.21 -4.49 4.27 8.04 11.85 20.36 1523 1.00
## theta[3] 5.93 0.32 8.53 -14.03 1.73 6.65 10.79 22.24 721 1.00
## theta[4] 7.72 0.16 6.48 -5.25 4.13 7.61 11.65 20.86 1673 1.00
## theta[5] 5.17 0.17 6.39 -9.63 1.53 5.91 9.32 16.69 1395 1.00
## theta[6] 6.15 0.18 6.59 -8.71 2.45 6.53 10.06 18.76 1326 1.00
## theta[7] 10.65 0.20 6.84 -1.38 6.18 9.84 14.55 26.31 1189 1.00
## theta[8] 8.73 0.34 8.71 -7.98 4.05 8.32 13.01 27.84 655 1.00
## lp__ -4.77 0.14 2.81 -10.99 -6.49 -4.53 -2.79 0.19 375 1.01
##
## Samples were drawn using NUTS(diag_e) at Mon Sep 7 16:35:59 2015.
## 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)
これらの結果が何を表すか，はおいおい説明するとして，動くかどうかの確認をこれで済ませてください。
うまくいったら，Shinystanでの結果の確認の仕方も確認しておきましょう。
require(shinystan)
launch_shinystan(fit)Stanを使ってみよう
ちょっとその前に
RStudioはプロジェクト単位でソースコードを管理することができる。 今回の分析用に，新しくプロジェクトを作り，そこに様々なソースコードをいれておくようにしよう。
- FileのNew Projectを選ぶ
- New Directory(新しいフォルダを作る)を選ぶ。今までのフォルダに追加する場合はExisting Directoryを選ぶ
- New Directoryを選んだら，プロジェクトを置いておきたいところを選んでCreateを押す。
- Rstudioの画面で，上のタグに選んだプロジェクト名が入っていたら成功
- FileのNew FileからR Scriptを選ぶと真っ白いメモ帳がでてくるので，ここにコードを書く。
- 外部ファイルを読み込む(データファイルを読み込むなど)場合は，作ったプロジェクトフォルダの中に置くようにする。
- Rstudioの四つの窓のどこかにあるFilesタブを押すと，プロジェクトフォルダの中身が見られるよ。
世界一簡単なstanコード
データを正規分布に基づいて発生させ，そのデータがどんな正規分布から出てきたか推定させる，というコードです。 答えがわかっている問題を解かせるわけだから，実質的な意味はないんだけど，stanのコードがどういうものなのかを知るために。
require(rstan)
n <- 100
mu <- 50
sig <- 10
y <- rnorm(n,mu,sig)
hist(y)
stancode <- '
data{
int<lower=0> T;
real N[T]; // data
}
parameters {
real mu;
real<lower=0> s2;
}
model{
N~normal(mu,sqrt(s2));
s2~cauchy(0,5);
}
'
datastan <- list(N=y,T=n)
fit <- stan(model_code = stancode,data=datastan,iter=5000,chain=4)## COMPILING THE C++ CODE FOR MODEL '2252c37f415db79c9e3f3198c1fa1980' NOW.traceplot(fit,ask=T)
print(fit,digit=3)## Inference for Stan model: 2252c37f415db79c9e3f3198c1fa1980.
## 4 chains, each with iter=5000; warmup=2500; thin=1;
## post-warmup draws per chain=2500, total post-warmup draws=10000.
##
## mean se_mean sd 2.5% 25% 50% 75% 97.5%
## mu 48.990 0.012 0.911 47.207 48.368 48.997 49.598 50.798
## s2 85.285 0.160 12.078 65.066 76.709 84.161 92.679 112.005
## lp__ -274.039 0.018 0.983 -276.712 -274.425 -273.735 -273.345 -273.084
## n_eff Rhat
## mu 6146 1.000
## s2 5724 1.000
## lp__ 2918 1.001
##
## Samples were drawn using NUTS(diag_e) at Mon Sep 7 16:36:27 2015.
## 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).世界で二番目に簡単なstanコード
今度はもう少し複雑なコード。 二番目に簡単なrstanコードからの丸パクリです。
今回は回帰分析モデルにもとづいてデータを発生させ，それをstanに解かせるというもの。 これもコードの練習でしかないんだけど，「データ生成モデルから統計モデルを見る」という感覚に開眼してください。つまり，現実的には生成モデルはわからないデータだけがあるわけで，データの奥を見通す目を持つというか，そういうことをしているだという理解の仕方を。
N <- 1000
x <- rnorm(N, mean = 50, sd = 10)
y <- 10 + 0.8 * x + rnorm(N, mean =0, sd = 7)
plot(x,y)
abline(lm(y~x),col=2)
stancode <- '
data{
int<lower=0> N;
real x[N];
real y[N];
}
parameters {
real alpha;
real beta;
real<lower=0> s;
}
model{
for(i in 1:N)
y[i] ~ normal(alpha + beta * x[i], s);
alpha ~ normal(0, 100);
beta ~ normal(0, 100);
s ~ cauchy(0,5);
}
'
datastan <- list(N=N, x=x, y=y)
fit <- stan(model_code=stancode, data=datastan, iter=1000, chain=4)## COMPILING THE C++ CODE FOR MODEL '68cb52b4d2351820fa0138861c46e315' NOW.traceplot(fit, ask=T)
print(fit, digit=2)## Inference for Stan model: 68cb52b4d2351820fa0138861c46e315.
## 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%
## alpha 10.59 0.05 1.12 8.42 9.80 10.58 11.36 12.80
## beta 0.79 0.00 0.02 0.74 0.77 0.79 0.80 0.83
## s 7.36 0.01 0.17 7.03 7.24 7.36 7.46 7.69
## lp__ -2495.45 0.06 1.22 -2498.70 -2495.96 -2495.14 -2494.58 -2494.09
## n_eff Rhat
## alpha 582 1.00
## beta 582 1.00
## s 578 1.01
## lp__ 483 1.00
##
## Samples were drawn using NUTS(diag_e) at Mon Sep 7 16:36:56 2015.
## 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).基本はこういうスタイル。モデルとデータ生成分布を考えてstanに推定させてやる，という考え方です。 どういうモデルがいいのか(あるのか),どういう分布がいいのか(あるのか)というのはおいおいわかっていくと思います。
Enjoy!
追伸；RStanの導入と基本的な使い方についてのスライドも見ておいてね。