-
-
Notifications
You must be signed in to change notification settings - Fork 266
RStan Getting Started (Japanese)
この文書はRStan Getting Startedを日本語に翻訳したものです。一部改変しています。
RStanはStanをRから使うためのインターフェースです。Stanについてもっと知りたい場合はStanのWebサイトを見てください。http://mc-stan.org/
Rのバージョンは4.2.0かそれ以降を強く推奨します。必要あれば、ここから最新版のRをインストールできます。 あと、ユーザにはRStudioのバージョン1.4.x以降の使用を強くオススメします。なぜならStanのコーディングを強力に支援する機能があるからです。
RStanの最新の開発版をインストールするには、以下を実行します。
# もしすでにrstanをインストールしているならば次の行を実行してください
# remove.packages(c("StanHeaders", "rstan"))
install.packages("rstan", repos = c("https://mc-stan.org/r-packages/", getOption("repos")))
このバージョンは公式ではリリースされていないものの、CRANの最新バージョンよりも新しいStan言語にアクセスできます。Stanコードをコンパイルするのに必要な環境を持っているか確認するには、「C++コンパイラ」と「インストールの確認」の節を見てください(下の「RStanのインストール」の節はスキップしてください)。
RStanのインストールの前に、RからC++コードをコンパイルできるようにする必要があります。OSに応じて以下のリンクに従ってください。
- Windows - Configuring C++ Toolchain
- Mac - Configuring C++ Toolchain
- Linux - Configuring C++ Toolchain
(訳注)上記のリンク先までは訳していなくてすみません。基本的に最新のRを使うと問題が少ないです。Windowsなら最新のRtoolsをインストールするだけで大丈夫です。Macならmacrtoolsパッケージを使います。上の指示に従ってもインストールがうまくいかない場合、Stan Discourse Forumで聞くか、日本語がよければr-wakalangというslackのrstanチャンネルで聞くとアドバイスがもらえるかもしれません。その場合、OSの情報、RStudioを使っているか否か、指示に従った結果の出力の情報をあわせて教えて下さい。
安全のため、以下の手順ですでに存在しているRStanパッケージを削除した方がよい場合があります。
remove.packages("rstan")
if (file.exists(".RData")) file.remove(".RData")
それからRを再起動します。
もしMac上でR 3.xを使っているなら、ここを見てください。それ以外の場合、通常、以下を実行することでインストールできます(正確にこのまま入力してください)。
Sys.setenv(DOWNLOAD_STATIC_LIBV8 = 1) # Linuxでnodejsライブラリがないときだけ必要
install.packages("rstan", repos = "https://cloud.r-project.org/", dependencies = TRUE)
無事インストールできたか確認するには、以下のRコードでRStanのexample/testモデルを実行します。
example(stan_model, package = "rstan", run.dontrun = TRUE)
モデルがコンパイルされてサンプリングが開始されるはずです。
パッケージ名は「rstan」(すべて小文字)です、だから library(rstan)
でrstanパッケージを読み込む必要があります。
library(rstan) # 起動時のメッセージが表示される
起動時のメッセージにあるように、もしマルチコアでメモリも十分にあるローカルマシンでrstanを使って並列で推定計算を実行させたい場合には、以下を実行します。
options(mc.cores = parallel::detectCores())
起動時の2番目のメッセージに従って
rstan_options(auto_write = TRUE)
を実行すると、(プログラムが変わらない限り)再コンパイルをしなくてすむように、コンパイル後のStanプログラムをハードディスクに保存します。これらのコマンドはrstan
パッケージを読み込むたびに実行する必要があるでしょう。
最後に、Windowsでは起動時の3番目のメッセージが表示され、--march=native
というコンパイラフラグを使わないように言ってくるでしょう。これまでのステップに従っていてMakevars.win
ファイルがこのフラグを含んでいない場合、このメッセージを無視してかまいません。
これはGelman et al (2003)の5.5節にある例です。8つの学校における教育の効果を研究したものです。簡単のため、この例を「eight schools」と呼びます。
はじめに、モデルをStanプログラムで書きます。RStudioを使っている場合は、メニューの[File]-[New File]-[Stan File]をクリックします。そうでない場合には、あなたのお気に入りのエディタを起動してください。そして、以下の内容をファイルに入力してschools.stan
という名前でRの作業ディレクトリに保存します。作業ディレクトリはgetwd()
を実行することで確認できます。
data {
int<lower=0> J; // 学校の数
array[J] real y; // 推定されている教育の効果
array[J] real<lower=0> sigma; // 教育の効果の標準誤差
}
parameters {
real mu; // 処置の効果(全体平均)
real<lower=0> tau; // 処置の効果の標準偏差
vector[J] eta; // 学校ごとのスケール前のバラつき
}
transformed parameters {
vector[J] theta = mu + tau * eta; // 学校ごとの処置の効果
}
model {
target += normal_lpdf(eta | 0, 1); // 事前分布の対数密度
target += normal_lpdf(y | theta, sigma); // 対数尤度
}
Stanプログラムの最後は(スペースやコメントなど一切の文字を含まない)空行で終わることに注意して下さい。
このモデルでは、theta
をparametersとして直接的に宣言する代わりに、mu
とeta
を使ってtheta
を作りtransformed parametersにしています。このようにパラメータ化することでより効率的にサンプリングできます(詳しい説明はこちらを参照してください)。
そして以下の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コードでモデルをフィッティングさせることができます。もし、作業ディレクトリにStanプログラムを作成しなかった場合には、引数のfile =
はその場所(ファイルのパス)を指定する必要があります。
fit <- stan(file = 'schools.stan', data = schools_dat)
stan
関数から返されたfit
オブジェクトはstanfit
クラスのS4オブジェクトです。print
・plot
・pairs
のような関数はフィットした結果と関連しており、このあとのコードを使ってfit
の中の結果を取り出すことができます。print
はモデルのパラメータに加えてlp__
という名前の対数事後確率(log-posterior)の要約を表示します(このあとの出力例を見てください)。他の関数やstanfit
クラスの詳細が知りたい場合は、stanfit
クラスのヘルプを見てください。
特に、MCMCサンプルを得るにはstanfit
オブジェクトに対してextract
関数を使います。 extract
はstanfit
オブジェクトから、興味あるパラメータのarraysのlistもしくは1つのarrayとしてMCMCサンプルを取り出します。さらに、as.array
・as.matrix
・as.data.frame
というS3関数が定義されています(R上でhelp("as.array.stanfit")
を使ってヘルプを見てください)。
print(fit)
plot(fit)
pairs(fit, pars = c("mu", "tau", "lp__"))
la <- extract(fit, permuted = TRUE) # arraysのlistを返す
mu <- la$mu
### iterations, chains, parametersの3次元arrayを返す
a <- extract(fit, permuted = FALSE)
### stanfitオブジェクトにS3関数を使う
a2 <- as.array(fit)
m <- as.matrix(fit)
d <- as.data.frame(fit)
Ratsの例も人気のある例です。 例えば、ここにOpenBUGSバージョンがあり、オリジナルはGelfand et al (1990)です。30匹のラットの体重を1週間ごと5週間にわたって記録したデータです。データは以下の表のようになっており、体重を記録した日付を表すのに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 <- 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', data = list(N=N, T=T, y=y, x=x, xbar=xbar))
開発チームがBUGS example(WinBUGSに付属している例)の大部分やその他のモデルをStanで作成しており、以下を実行することで試すことができます。
model <- stan_demo()
このあとはモデルの例をポップアップで出てきたリストから選んでいきます。はじめてstan_demo()
を呼ぶときは、これらの例をダウンロードするかどうかを聞いてきます。option 1を選んでrstanがインストールされたディレクトリに保存するようにしましょう。そうすると、またいつか実行する時に再ダウンロードしなくてもすみます。含まれるmodel
オブジェクトはstanfit
クラスのインスタンスなので、それらに対して、 print
・plot
・pairs
・extract
を呼び出すことができます。
RStanに関するさらなる詳細はrstanパッケージのvignetteに含まれるドキュメントで見ることができます。例えば、help(stan)
やhelp("stanfit-class")
を使うことでstan
関数やS4クラスであるstanfit
のヘルプを見ることができます。
StanのサンプラーやオプティマイザーやStanの文法に関する詳細を知りたい時にはStan's modeling language manualを見てください。
さらに、Stanの使い方を議論したい場合にはStan Discourse Forumに例を投稿したり、(R)Stanについて質問してもよいでしょう。助けが必要な時は、以下に関する十分な情報もあわせて記述することが重要です。
- Stanのモデルコード
- データ
- 実行するのに必要なRコード
-
stan
関数を呼ぶときにverbose=TRUE
かつcores=1
を指定して出力されたエラーメッセージのdump - C++コンパイラのバージョン。もし
gcc
を使っているならばg++ -v
で得ることができます。 - R上で
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