RStan Getting Started (Japanese)

Kentaro Matsuura edited this page Feb 22, 2018 · 32 revisions

RStanをはじめよう

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

Latest Version: 2.17.3 (20 Jan 2018)

下の説明のほとんどはRStanの上記のバージョンに関するものです。

導入

RStanはStanをR から利用するためのインターフェースです。この文書では以下の説明をします:

Stanやその文法のさらなる情報が欲しい場合は以下の公式Webサイトを見てください:

RStanをインストールする前準備

R

Rのバージョンは3.0.2かそれ以降が必要です(Macの場合、もしXcodeをアップグレードしているとRのほか、様々なライブラリの再インストールが必要になるかもしれません)。しかし、RStanは3.2.0より前のバージョンだと何らかの点でうまく動かないことが知られています(ggplot2でのプロットが失敗するでしょう。また、rstan.package.skeleton関数が必要なファイルをダウンロードしないなど)。バージョン3.0.2というのは、最新のRが用意できないようなリモートサーバーで作業するために必要なバージョンです。RStanはRの最新の安定板とdevelop版でテストしています。もし、あなたが管理者権限を持っているなら、特にRStanのワークショップに参加する場合、Rの最新版を使わない理由はありません。最新版のRは以下から入手できます。

https://www.r-project.org/

左のメニューの「Download」のところにある「CRAN」をクリックして、それからミラーサイトを選びます(http://cran.rstudio.com/ をオススメします。それは最も近くて信頼できるミラーサイトへリダイレクトしてくれるからです)。それからあなたのOS(Windows・Linux・Mac)にあわせてリンクをクリックします。Windowsだとダウンロードの前に「base」パッケージを選ぶもうひと手間があります。

LinuxとMacの場合はRのコマンドラインとGUIをインストールしましょう。デフォルトの設定のままで大丈夫です。

RStudio

必須ではありませんが、ほとんどのユーザにはRStudioをインストールすることを強くオススメします。.stanファイル(Stanのモデルを書くファイル)に完璧ではありませんが対応しており、(もしRStudioのバージョンが0.99.1259以降なら)Stan 2.10以降の文法にあわせてシンタックスハイライトもされます。

C++コンパイラ

RtoolsはRとうまく連携できるC++コンパイラを提供します。Rtoolsは主にGNU makeとGNU gccとUNIXライクなプラットフォームで一般的に使われるコマンド群を含みます。

以降ではWindowsの場合に、C++コンパイラをダウンロードしてインストールする方法を説明します。

  • Macの場合、RStan Mac OS X Prerequisite Installation Instructions を見てください。
  • Linuxの場合、使っているディストリビューションのパッケージマネージャーを使ってビルドに必要なものと最近のg++かclang++をインストールしてください。

Rtoolsのダウンロード

Rtoolsを http://cran.r-project.org/bin/windows/Rtools/ からインストールできます。もし、インストールしたRのバージョンが3.3.xならば、Rtools33.exeをダウンロードして実行してください。もし、バージョンが3.2.xかそれより下ならば、Rtools33.exeを使わないで、インストールしたRのバージョンにあったなかで最も新しい「frozen」になっているRtoolsを選んでください。もし、複数のバージョンのRを使いたいならば、複数のバージョンのRtoolsを入れて使うことも可能です。しかし、Stanのプログラムをコンパイルするために何らかの環境変数をいじる必要があるでしょう。

Rtoolsのインストール

Rtoolsをインストールする際に、環境変数PATHを編集するステップに注意してください。以下の画像で示したように、このオプションにチェックを入れましょう(何もタイプする必要はありません。ただチェックするだけです)。こうすることで、RからC++コンパイラを使えるようになります。

editpathrtools

RからRtoolsが使えるか確認する

  • 新しくRを起動します。
  • Sys.getenv("PATH")を実行して、環境変数PATHにRtoolsが含まれているかチェックします。もし、Rtoolsをc:\\Rtoolsにインストールした場合、以下のような出力が見えるはずです。
> Sys.getenv('PATH')
[1] "c:\\\\Rtools\\\\bin;c:\\\\Rtools\\\\gcc-4.6.3\\\\bin;...

上の出力の例では、はじめの方がc:\\Rtools\\bin and c:\\Rtools\\gcc-4.x-y\\binとなっています。もし、Rtoolsをc:\\Rtoolsにインストールして、これら2つが出てこなかったら、上のRtoolsのインストールの節に書いてあることにきちんと従ったか確かめてください。もし、うまくいかなかったら、選択肢は二つあります。一つ目はRtoolsを再インストールすることで、もう一つは環境変数PATHを手動で編集することです。簡単なので、前者の方をまずはオススメします。

  • Rから本当にg++を呼び出せるかチェックします。例えば、以下のようにしてRからgccのバージョンを見ることができます。
> system('g++ -v')
Using built-in specs.
COLLECT_GCC=c:\Rtools\GCC-46~1.3\bin\G__~1.EXE
COLLECT_LTO_WRAPPER=c:/rtools/gcc-46~1.3/bin/../libexec/gcc/i686-w64-mingw32/4.6.3/lto-wrapper.exe
Target: i686-w64-mingw32
Configured with: /data/gannet/ripley/Sources/mingw-test3/src/gcc/configure --host=i686-w64-mingw32 --build=x86_64-linux-gnu --target=i686-w64-mingw32 --with-sysroot=/data/gannet/ripley/Sources/mingw-test3/mingw32mingw32/mingw32 --prefix=/data/gannet/ripley/Sources/mingw-test3/mingw32mingw32/mingw32 --with-gmp=/data/gannet/ripley/Sources/mingw-test3/mingw32mingw32/prereq_install --with-mpfr=/data/gannet/ripley/Sources/mingw-test3/mingw32mingw32/prereq_install --with-mpc=/data/gannet/ripley/Sources/mingw-test3/mingw32mingw32/prereq_install --disable-shared --enable-static --enable-targets=all --enable-languages=c,c++,fortran --enable-libgomp --enable-sjlj-exceptions --enable-fully-dynamic-string --disable-nls --disable-werror --enable-checking=release --disable-win32-registry --disable-rpath --disable-werror CFLAGS='-O2 -mtune=core2 -fomit-frame-pointer' LDFLAGS=
Thread model: win32
gcc version 4.6.3 20111208 (prerelease) (GCC)

> system('where make')
c:\Rtools\bin\make.exe

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

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

  • Windowsの場合、ここを見てください。
  • MacかLinuxの場合、ここを見てください。

RStanのインストール方法

やっとRStanをインストールする準備ができました!以下のステップに順に従うとRStanをインストールできます。

  • Rを起動します(R GUIでもターミナル上で起動したRでもオススメしたRStudioを起動して実行してもOKです)。
  • rstanパッケージの最新版とそれが依存しているパッケージを以下のように明示的に指定してCRANからインストールします。
# もしhttpsを使ったダウンロードに対応していないならば、httpsの「s」を取って実行してください
install.packages('rstan', repos='https://cloud.r-project.org/', dependencies=TRUE)

dependencies=TRUEを省略しないでください。

  • Stanがちょうどバージョンアップをしている最中(StanのリリースはされたけどCRANのrstanのバイナリの最新版ができあがっていない状態)の場合には、上記のinstall.packages関数を実行する代わりにここを読んでインストールする必要があります。
  • それでもうまくいかない場合、ソースからインストールすることを試してもよいかもしれません。
# 注:以下の「4」をビルドに使いたいコア数に置き換えてください
Sys.setenv(MAKEFLAGS = "-j4")
install.packages("rstan", type = "source")
  • Rの再起動 インストールの後はたぶんRの再起動が必要です。 そしてrstanパッケージを呼び出す前に以前のバージョンのRStanによって作られた何らかのオブジェクトが勝手にloadされないか確かめてください。

  • Rで以下を実行して、C++が動くか確かめてください。

fx <- inline::cxxfunction( signature(x = "integer", y = "numeric" ) , '
  return ScalarReal( INTEGER(x)[0] * REAL(y)[0] ) ;
')
fx( 2L, 5 ) # 10になるはずです

RStanの使用方法

rstanのload

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

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

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

rstan_options(auto_write = TRUE)
options(mc.cores = parallel::detectCores())

これらのoptionはそれぞれ、再コンパイルをしなくてすむようにコンパイル後のStanプログラムをハードディスクに保存すること、複数のマルコフ連鎖を並列で実行することを指定しています。

例 1: Eight Schools

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

はじめに、モデルをStanプログラムで書いて8schools.stanという名前で保存します。中身は以下の通りです(ファイルはここで見ることができます)。

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

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 {
  target += normal_lpdf(eta | 0, 1);
  target += normal_lpdf(y | theta, sigma);
}

このモデルでは、thetaをparametersとして直接的に宣言する代わりに、muetaを使ってthetaを作りtransformed parametersにしています。このようにパラメータ化することでより効率的にサンプリングできます(詳しい説明はこちらを参照してください)。8schools.stanを作業ディレクトリに置いたとすると、以下の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コードでモデルをフィッティングさせることができます。

fit <- stan(file = '8schools.stan', data = schools_dat,
            iter = 1000, chains = 4)

stan関数のmodel_code引数を使うことでStanのモデルをRコード内の文字列で指定することもできます。しかしながら、これはオススメできません。なぜなら、ファイルを使えばprint文が使えたり、エラーメッセージの中の行番号を参考にすることができるからです。

stan関数から返されたfitオブジェクトはstanfitクラスのS4オブジェクトです。 printplotpairsのような関数はフィットした結果と関連しており、このあとのコードを使ってfitの中の結果を取り出すことができます。printはモデルのパラメータに加えてlp__という名前の対数事後確率(log-posterior)の要約を表示します(このあとの出力例を見てください)。他の関数やstanfitクラスの詳細が知りたい場合は、stanfitクラスのヘルプを見てください。

特に、MCMCサンプルを得るにはstanfitオブジェクトに対してextract関数を使います。 extractstanfitオブジェクトから、興味あるパラメータのarraysのlistもしくは1つのarrayとしてMCMCサンプルを取り出します。さらに、as.arrayas.matrixas.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)
> print(fit, digits = 1)
Inference for Stan model: 8schools.
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.2     0.2 5.4  -1.9   4.8   8.1  11.3  19.3   480    1
tau        6.8     0.3 6.2   0.3   2.5   5.2   9.2  21.7   425    1
eta[1]     0.4     0.0 1.0  -1.5  -0.3   0.4   1.0   2.2  2000    1
eta[2]     0.0     0.0 0.8  -1.7  -0.6   0.0   0.5   1.7  2000    1
eta[3]    -0.2     0.0 1.0  -2.1  -0.9  -0.2   0.4   1.7  2000    1
eta[4]    -0.1     0.0 0.9  -1.8  -0.7  -0.1   0.5   1.7  2000    1
eta[5]    -0.4     0.0 0.9  -2.1  -1.0  -0.4   0.2   1.4  2000    1
eta[6]    -0.2     0.0 0.9  -1.9  -0.8  -0.2   0.4   1.5  1731    1
eta[7]     0.3     0.0 0.9  -1.4  -0.2   0.4   0.9   2.0  1507    1
eta[8]     0.0     0.0 0.9  -1.9  -0.6   0.0   0.7   1.8  1988    1
theta[1]  11.5     0.3 8.8  -2.4   5.9  10.1  15.6  32.9   977    1
theta[2]   7.8     0.1 6.2  -4.7   4.1   7.9  11.6  20.3  2000    1
theta[3]   6.1     0.2 7.7 -11.2   2.1   6.4  10.5  20.2  2000    1
theta[4]   7.6     0.1 6.5  -4.9   3.8   7.8  11.4  21.3  2000    1
theta[5]   5.0     0.1 6.6  -9.3   1.2   5.6   9.3  16.7  2000    1
theta[6]   6.2     0.2 6.7  -8.2   2.2   6.5  10.5  18.5  2000    1
theta[7]  10.8     0.2 7.0  -1.3   6.1  10.1  15.1  26.8  2000    1
theta[8]   8.7     0.2 8.2  -7.3   3.9   8.4  12.8  27.2  1446    1
lp__     -39.5     0.1 2.6 -45.1 -41.2 -39.4 -37.7 -35.1   590    1

Samples were drawn using NUTS(diag_e) at Fri May  5 10:41:43 2017.
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).

さらに、BUGS(やJAGS)と同様に、CmdStan(Stanのコマンドラインインターフェース)を使うにはすべてのデータをRのdumpファイルの形式で持っておく必要があります。dumpファイルがある場合、rstanにはread_rdumpという関数があり、これを使って全てのデータを1つのRの名前付きリストとして読み込むことができます。例えば、作業ディレクトリに以下の中身を持つ "8schools.rdump" という名前のファイルがあるとします。

J <- 8
y <- c(28,  8, -3,  7, -1,  1, 18, 12)
sigma <- c(15, 10, 16, 11,  9, 11, 10, 18)

そうすると以下のように "8schools.rdump" からデータを読み込むことができます。

schools_dat <- read_rdump('8schools.rdump')

dumpファイルはRのsource関数を使って読み込むこともでき、dumpファイル内の各変数はグローバル環境内のオブジェクトとなります。この場合、Stanを実行する時のstan関数のdata引数を省略することができます。stan関数は 8schools.stan のdataブロックにある変数名と同じ名前のオブジェクトを呼び出された環境から探すからです。すなわち、

source('8schools.rdump')
fit <- stan(file = '8schools.stan', iter = 1000, chains = 4)

例 2: Rats

Ratsの例も人気のある例です。 例えば、ここOpenBUGSバージョンがあり、オリジナルはGelfand et al (1990)です。30匹のラットの体重を1週間ごと5週間にわたって記録したデータです。データは以下の表のようになっており、体重を記録した日付を表すのにxを使っています。この例を試すのにリンク先のデータとモデルのコードを使うことができます(rats.txtrats.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クラスのインスタンスなので、それらに対して、 printplotpairsextractを呼び出すことができます。

さらなるヘルプ

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
You can’t perform that action at this time.
You signed in with another tab or window. Reload to refresh your session. You signed out in another tab or window. Reload to refresh your session.
Press h to open a hovercard with more details.