Skip to content

RStan Getting Started (Japanese)

Jump to bottom
Kentaro Matsuura edited this page Aug 26, 2020 · 45 revisions

RStanをはじめよう

この文書はRStan Getting Startedを日本語に翻訳したものです。一部構成を改変しています。

RStanはStanをRから使うためのインターフェースです。Stanについてもっと知りたい場合はStanのWebサイトを見てください。http://mc-stan.org/

Latest Version: 2.21.1 (2020年7月)

下の説明のほとんどは最新バージョンのRStanに関するものです。Rのバージョンは3.4.0かそれ以降が必要です。必要あれば、ここから最新版のRをインストールできます。

あと、ユーザにはRStudioのバージョン1.2.x以降の使用を強くオススメします。なぜならStanのコーディングを強力に支援する機能があるからです。

RStanのインストール

安全のため、以下の手順ですでに存在しているRStanパッケージを削除した方がよい場合があります。

remove.packages("rstan")
if (file.exists(".RData")) file.remove(".RData")

それからRを再起動します。

もしMacのCatalinaを使っているなら、ここを見てください。それ以外の場合、通常、以下を実行することでインストールできます（正確にこのまま入力してください）。

install.packages('rstan', repos='https://cloud.r-project.org/', dependencies=TRUE)

しかしながら、もしLinuxを使っているか、getOption("pkgType")を"source"に設定しているか、RがRStanの最新バージョンをソースからインストールしたいか聞いてくる場合には、対応するWindows, Mac, Linuxのページを読んで下さい。

C++コンパイラ

C++コンパイラがインストールされているかチェックするために、できればRStudio上で、またはRコンソール上で以下を実行します。

pkgbuild::has_build_tools(debug = TRUE)

それから、

  1. WindowsとRStudioを使用している場合（推奨）、もしRtoolsがない場合にはポップアップウィンドウが出てきてRtoolsをインストールするか聞いてくるでしょう。Yesをクリックして1分ぐらいインストールが終わるのを待って下さい。
  2. WindowsでRStudioを使用していない場合、Rコンソール上でRtoolsのインストールを指示するメッセージが出てくるのでそれに従って下さい。詳しくはRtoolsのダウンロードとインストールが参考になると思いますが、ふつうは"Installing RStan from source"（ソースからのインストール）の節で書かれていることが必要になることはありません。
  3. Macを使っている場合、リンクが出てくるかもしれません。しかし、それをクリックしないでこれに従って下さい。
  4. Linuxを使っている場合、これに従って下さい。

上の指示に従ってもインストールがうまくいかない場合、Stan Discourse Forumで聞くか、日本語がよければr-wakalangというslackのrstanチャンネルで聞くとよいでしょう。その場合、OSの情報、RStudioを使っているか否か、指示に従った結果の出力の情報をあわせて教えて下さい。

コンパイル設定のカスタマイズ（熟練者だけ）

以下は必須ではありません。やらなくてもRStanは動くからです。しかし、コンパイル設定のカスタマイズを行うと速度向上が期待できます。設定するにはR上で以下を一度だけ実行します。

dotR <- file.path(Sys.getenv("HOME"), ".R")
if (!file.exists(dotR)) dir.create(dotR)
M <- file.path(dotR, ifelse(.Platform$OS.type == "windows", "Makevars.win", "Makevars"))
if (!file.exists(M)) file.create(M)
cat(if( grepl("^darwin", R.version$os)) "\nCXX14FLAGS=-O3 -march=native -mtune=native -arch x86_64 -ftemplate-depth-256" else 
    if (.Platform$OS.type == "windows") "\nCXX14FLAGS=-O3 -mtune=native -mmmx -msse -msse2 -msse3 -mssse3 -msse4.1 -msse4.2" else
    "CXX14FLAGS = -fPIC",
    file = M, sep = "\n", append = TRUE)

しかしながら、最適化レベルをO3にすると他のパッケージで問題を起こすことがあります。また、-march=native -mtune=nativeを指定するとStanのプログラムが動かないことがあります。特にWindowsでは-march=nativeが問題を起こすことが知られています。もしC++コンパイラの設定を変えたい場合には、以下を実行してください。

M <- file.path(Sys.getenv("HOME"), ".R", ifelse(.Platform$OS.type == "windows", "Makevars.win", "Makevars"))
file.edit(M)

もしくは、~/.R以下にあるテキストファイルを直接編集しましょう。Windowsでは、この場所はたいていC:/Users/USERNAME/Documents/.Rです。最後に、Windowsでは今のところコンパイラフラグを指定する必要はないはずなので、問題がある場合はMakevars.winファイルを完全に削除してください。

RStanの使用方法

rstanのload

パッケージ名は「rstan」（すべて小文字）です、だから library(rstan) でrstanをloadする必要があります。

library(rstan) # スタートアップのメッセージが表示される

起動時のメッセージにあるように、もしマルチコアでメモリも十分にあるローカルマシンでrstanを使って並列で推定計算を実行させたい場合には、以下を実行します。

options(mc.cores = parallel::detectCores())

起動時の2番目のメッセージに従って

rstan_options(auto_write = TRUE)

を実行すると、（プログラムが変わらない限り）再コンパイルをしなくてすむように、コンパイル後のStanプログラムをハードディスクに保存します。

最後にWindowsの場合には、起動時の3番目のメッセージに従って次を実行するとよいでしょう。

Sys.setenv(LOCAL_CPPFLAGS = '-march=corei7')

もし前節でコンパイル設定のカスタマイズを行っている場合には、この実行は必要ありません。

例 1: Eight Schools

これはGelman et al (2003)の5.5節にある例です。8つの学校における教育の効果を研究したものです。簡単のため、この例を「eight schools」と呼びます。

はじめに、モデルをStanプログラムで書きます。RStudioを使っている場合は、メニューの[File]-[New File]-[Stan File]をクリックします。そうでない場合には、あなたのお気に入りのエディタを起動してください。そして、以下の内容をファイルに入力して8schools.stanという名前でRの作業ディレクトリに保存します。作業ディレクトリはgetwd()を実行することで確認できます。

data {
  int<lower=0> J;         // 学校の数
  real y[J];              // 推定されている教育の効果
  real<lower=0> sigma[J]; // 教育の効果の標準誤差
}

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 = '8schools.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)

例 2: Rats

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

例 3: 何でも！

開発チームが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に関する情報

References

  • 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
Clone this wiki locally