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インターフェイス)があります。
stanのプロモーション動画もCool!だぜ。
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 8.05 0.17 5.12 -1.94 4.70 7.87 11.35 18.44 904 1.00
## tau 6.23 0.22 5.32 0.22 2.21 4.98 8.91 19.47 608 1.00
## eta[1] 0.39 0.02 0.97 -1.54 -0.25 0.44 1.05 2.18 1559 1.00
## eta[2] 0.00 0.02 0.94 -1.81 -0.62 0.00 0.61 1.92 1615 1.00
## eta[3] -0.17 0.03 0.95 -2.03 -0.79 -0.18 0.43 1.74 1103 1.00
## eta[4] -0.06 0.02 0.90 -1.77 -0.65 -0.06 0.54 1.77 1402 1.00
## eta[5] -0.35 0.02 0.90 -2.04 -0.97 -0.39 0.24 1.56 1293 1.00
## eta[6] -0.20 0.02 0.92 -2.00 -0.78 -0.21 0.38 1.67 1502 1.00
## eta[7] 0.33 0.03 0.91 -1.51 -0.22 0.33 0.94 1.99 1297 1.00
## eta[8] 0.05 0.02 0.93 -1.84 -0.59 0.09 0.72 1.81 1417 1.00
## theta[1] 11.31 0.24 8.05 -2.56 6.12 10.45 15.31 30.26 1108 1.00
## theta[2] 8.06 0.14 6.40 -4.05 3.91 7.92 12.27 20.97 2000 1.00
## theta[3] 6.28 0.22 7.71 -11.29 2.37 6.95 11.02 20.16 1242 1.00
## theta[4] 7.79 0.17 6.65 -5.56 3.76 7.70 11.72 21.94 1495 1.00
## theta[5] 5.15 0.15 6.33 -8.36 1.25 5.62 9.31 16.37 1759 1.00
## theta[6] 6.33 0.17 6.54 -7.48 2.59 6.52 10.59 18.58 1409 1.00
## theta[7] 10.51 0.18 6.60 -1.55 6.20 10.14 14.19 25.80 1357 1.00
## theta[8] 8.59 0.21 7.56 -6.09 4.18 8.24 12.68 25.68 1274 1.00
## lp__ -5.12 0.12 2.78 -11.37 -6.88 -4.88 -3.14 -0.51 535 1.01
##
## Samples were drawn using NUTS(diag_e) at Tue Sep 8 11:03:43 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 49.218 0.013 0.994 47.252 48.550 49.226 49.888 51.147
## s2 100.487 0.193 14.420 75.928 90.452 99.258 109.442 131.934
## lp__ -282.320 0.018 0.984 -284.959 -282.715 -282.019 -281.609 -281.347
## n_eff Rhat
## mu 5918 1.001
## s2 5607 1.000
## lp__ 2932 1.001
##
## Samples were drawn using NUTS(diag_e) at Tue Sep 8 11:04:10 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.51 0.06 1.18 8.25 9.73 10.53 11.31 12.85
## beta 0.79 0.00 0.02 0.74 0.78 0.79 0.81 0.84
## s 6.89 0.01 0.15 6.59 6.78 6.88 6.99 7.20
## lp__ -2428.87 0.07 1.30 -2432.28 -2429.38 -2428.53 -2427.95 -2427.46
## n_eff Rhat
## alpha 393 1.00
## beta 389 1.00
## s 650 1.00
## lp__ 345 1.01
##
## Samples were drawn using NUTS(diag_e) at Tue Sep 8 11:04:39 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の導入と基本的な使い方についてのスライドも見ておいてね。