RStan Getting Started (Japanese) - PrinceWangR/rstan GitHub Wiki

<wiki:toc max_depth="3" />

この文書はRStan Getting Startedを日本語に翻訳したものです。

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

RStanはStanをR (http://www.r-project.org/) から利用するためのインターフェースです。この文書では以下の説明をします:

• RStanをインストールする前準備
• RStanのインストールの仕方
• RStanの使い方

Stanやその文法のさらなる情報が欲しい時は以下を参照してください:

• Stan Home Page

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のバージョン0.99.441かそれ以降をインストールすることをオススメします。

https://www.rstudio.com/products/rstudio/

.stanファイル（Stanのモデルを書くファイル）にも完璧ではありませんが対応しています。

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

この小節に書いてあることは必須ではありません、やらなくてもRStanは動くからです。でも設定することを強くオススメします。もし個人のMakevarsファイルをまだ持っていないならこのリンク先に従って作りましょう。 以下の手順はRを起動して、そのMakevarsファイルをRで作る例です（R GUIでもターミナル上で起動したRでもオススメしたRStudioを起動して実行してもOKです）。

dotR <- file.path(Sys.getenv("HOME"), ".R")
if (!file.exists(dotR)) dir.create(dotR)
M <- file.path(dotR, "Makevars")
if (!file.exists(M)) file.create(M)
cat("\nCXXFLAGS=-O3 -mtune=native -march=native -Wno-unused-variable -Wno-unused-function",
    file = M, sep = "\n", append = TRUE)

ひとつ注意があります。Rのパッケージの中にはコンパイルの設定が平凡な場合だけをテストしているものがあって、上で設定したようにコンパイルの最適化レベルを3にすると、そういうパッケージをソースからインストールする時に不具合がでるかもしれません。

もしg++のバージョンが4.9かそれ以降ならば（Macでは珍しいですが）、以下をR上で実行することをオススメします。

cat("\nCXXFLAGS+=-flto -ffat-lto-objects",
    file = M, sep = "\n", append = TRUE)

さらに、OS Xの場合だけ （clang++はインストールしたとして）以下をR上で実行する必要があります。

cat("\nCC=clang", "CXX=clang++ -arch x86_64 -ftemplate-depth-256",
    file = M, sep = "\n", append = TRUE)

Rのバージョン3.3.xを使うと、C++11をサポートしているg++ 4.9.xに対応したWindows用Rtoolsをダウンロードして使うことができます。今のところStanではg++ 4.6とそれより前のバージョンではC++11をサポートしていません。しかし、g++の4.6よりあとのバージョンとclang++の最近のバージョンを使うとStanでもC++11は使えるとは思います。

C++11を使うかどうかに関わらず、もしRtools33（かそれ以降）を使うなら、以下を一度実行する必要があります。 （訳注：もしRtoolsをC:\Rtoolsにインストールした場合は不要です。Cドライブ直下でないところにインストールした場合、以下のように.RprofileでBINPREFを定義する、もしくはWindowsの環境変数を直接編集して設定する、もしくは(Rのインストールフォルダ)\etc\x64\Makeconfの19行目あたりのBINPREFを編集する必要があります）

cat('Sys.setenv(BINPREF = "C:/Rtools/mingw_$(WIN)/bin/")',
    file = file.path(Sys.getenv("HOME"), ".Rprofile"), 
    sep = "\n", append = TRUE)

もしg++のversion 6かそれ以上を使う場合、 Stanに関係がない警告が出力されるのを止めたいと思うかもしれません。それには次を実行します。

cat("\nCXXFLAGS += -Wno-ignored-attributes -Wno-deprecated-declarations", 
    file = M, sep = "\n", append = TRUE)

以下を実行することで、設定が正しいか確かめることができます。

cat(readLines(M), sep = "\n")

もし正しくなければ、ファイルのパスを以下で確かめて、そのファイルをテキストエディタで開いて編集することもできます。

cat(M)

• Rを起動します（R GUIでもターミナル上で起動したRでもオススメしたRStudioを起動して実行してもOKです）。

• ソースからビルドする場合だけ（OS XやWindowsでは珍しいと思いますが）、ビルドに何コア使うか指定します。例えば4コア使いたいならば、以下をR上で実行します。

Sys.setenv(MAKEFLAGS = "-j4")

• rstanパッケージの最新版とそれが依存しているパッケージを以下のように明示的に指定してインストールします。

# もしhttpsを使ったダウンロードに対応していないならば、httpsの「s」を取って実行してください
install.packages('rstan', repos = 'https://cloud.r-project.org/', dependencies = TRUE)

それでもうまくいかない場合、ソースからインストールすることを試してもよいかもしれません。

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」です、だから library(rstan) でrstanをloadする必要があります。

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

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

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

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

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

はじめに、このモデルを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として直接的に宣言する代わりに、muとetaを使って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))

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

stan関数のmodel_code引数を使うことでStanのモデルをRコード内の文字列で指定することもできます。しかしながら、これはオススメできません。

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という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関数のas.array（やas.matrix）を使う
a2 <- as.array(fit)
m <- as.matrix(fit)

> print(fit, digits = 1)
Inference for Stan model: schools_code.
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.9     0.2 4.9  -2.1  4.5  7.9 11.0  17.8   422    1
tau       6.3     0.3 5.0   0.2  2.5  5.2  8.9  18.7   214    1
eta[1]    0.4     0.0 0.9  -1.5 -0.2  0.4  1.0   2.1   928    1
eta[2]    0.0     0.0 0.9  -1.8 -0.6  0.0  0.5   1.8  1640    1
eta[3]   -0.2     0.0 1.0  -2.1 -0.8 -0.2  0.4   1.8  1243    1
eta[4]    0.0     0.0 0.9  -1.7 -0.6  0.0  0.6   1.7  1421    1
eta[5]   -0.3     0.0 0.9  -2.0 -1.0 -0.4  0.3   1.5   883    1
eta[6]   -0.2     0.0 0.9  -2.0 -0.8 -0.2  0.4   1.6   926    1
eta[7]    0.4     0.0 0.9  -1.4 -0.2  0.4  0.9   2.1   969    1
eta[8]    0.1     0.0 1.0  -1.8 -0.6  0.1  0.7   2.0  1365    1
theta[1] 11.4     0.3 8.1  -1.4  5.9 10.3 15.2  30.6   574    1
theta[2]  7.7     0.2 6.1  -3.7  3.9  7.8 11.4  19.5   762    1
theta[3]  5.8     0.3 7.9 -12.1  1.8  6.5 10.5  19.9   715    1
theta[4]  8.0     0.2 6.5  -5.4  3.9  8.1 12.3  20.2   977    1
theta[5]  5.0     0.3 6.7 -10.3  1.3  5.7  9.5  16.5   667    1
theta[6]  6.0     0.2 6.6  -8.4  2.0  6.2 10.2  18.6   976    1
theta[7] 10.8     0.3 6.8  -1.1  6.2 10.2 14.9  26.0   596    1
theta[8]  8.6     0.3 7.9  -6.1  4.0  8.1 12.6  27.7   629    1
lp__     -5.0     0.1 2.6 -10.7 -6.6 -4.8 -3.1  -0.5   367    1

Samples were drawn using NUTS2 at Fri Apr 12 22:09:54 2013.
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_y <- 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)

Ratsの例も人気のある例です。 例えば、ここにOpenBUGSバージョンがあり、オリジナルはGelfand et al (1990)です。30匹のラットの体重を1週間ごと5週間にわたって記録したデータです。データは以下の表のようになっており、体重を記録した日付を表すのにxを使っています。この例を試すのにリンク先のデータとモデルのコードを使うことができます（rats.txtとrats.stan）。

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

開発チームが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 User's Mailing listに例を投稿したり、(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.
• The Stan Development Team (2015). 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