【R言語】 C++ベースで早く動く MCMC Sampler {Stan}をRから動かす ~{RStan}パッケージのセットアップ方法

  • 13
    いいね
  • 0
    コメント
この記事は最終更新日から1年以上が経過しています。

PC 環境 】

OS : Windows 7 (64-bit)

R : version 3.1.0

1. RStudioR のインストール先を、パス名にスペースのないフォルダに移動

 → 「RStudio\bin」をCドライブ直下に移動した

2. {Rtools} パッケージを R にインストールする

http://cran.r-project.org/bin/windows/Rtools/ からダウンロード

※ 今回は、 Rtools31.exe をダウンロードした。

(ウェブサイトを見て、自分がいま使用している Rversion に適合するものを選ぶ)

Rtools.download.site.png

2-1. Rtools31.exe ファイルを、 C ドライブ直下に移動

2-2.以下を環境変数 Path を、パスの先頭に記述 ※パスの先頭にする

※ ③は、 お使いの R のバージョンに合わせたパス名を入力してください

   ① C:\Rtools\bin
   ② C:\Rtools\gcc-4.6.3\bin
   ③ C:\Program Files\R\R-3.1.0\bin
   ④ C:\RStudio\bin
   ⑤ C:\Windows\system32;

environment.var.setting.png

2-3. Rtools の動作チェック

(参考) (リンク先)Rtools 公式ページ 内の"Check if Rtools can be used in R"項目

Rtools.official.site.png

Rtools.official.site.2.png

RSys.getenv() 関数で、 R の中から、 Windows OS の環境変数 設定情報を表示させる

RStudio
Sys.getenv('PATH')

Rtools\bin と、 Rtools\gcc-4.6.3\bin が記載されていることをチェックする

(以下の実行結果は、自分が使用中の PC の環境変数の設定内容によって変わる)

RStudio
   [1] "C:\\Program Files\\R\\R-3.1.0\\bin\\x64;C:\\Rtools\\bin;C:\\Rtools\\gcc-4.6.3\\bin;C:\\Program Files\\R\\R-3.1.0\\bin;C:\\RStudio\\bin;C:\\Windows\\system32;C:\\Program Files\\Java\\jdk1.7.0_51\\bin;C:\\apache-maven-3.2.1\\bin;c:\\usr\\local\\bin;c:\\gs\\gs7.07\\bin;\\c:\\gs\\gs7.07\\lib;c:\\dviout;C:C:\\Program Files (x86)\\MeCab\\bin;C:\\Program Files (x86)\\MeCab;C:\\Python34;C:\\Python34\\Scripts;C:\\Rtools\\bin;C:\\Program Files\\R\\R-3.1.0\\bin;C:\\julia-e44b593905\\bin;C:\\Rtools\\bin;C:\\Rtools\\gcc-4.6.3\\bin;C:\\Program Files\\R\\R-3.1.0\\bin;C:\\RStudio\\bin;C:\\Windows\\system32;C:\\Program Files\\Java\\jdk1.7.0_51\\bin;C:\\apache-maven-3.2.1\\bin;c:\\usr\\local\\bin;c:\\gs\\gs7.07\\bin;\\c:\\gs\\gs7.07\\lib;c:\\dviout;C:C:\\Program Files (x86)\\MeCab\\bin;C:\\Program Files (x86)\\MeCab;C:\\Python34;C:\\Python34\\Scripts;C:\\Rtools\\bin;C:\\Program Files\\R\\R-3.1.0\\bin;C:\\julia-e44b593905\\bin;"

続けて、 C++ のコンパイラ g++Windows コンソールから呼び出せるかどうか確認

Rsystem() 関数で、 R の中から、 Windows OS コマンドプロンプトに命令指示

RStudio
system('g++ -v')

認識された g++ の情報が表示された。末尾に、 gcc version 4.6.3 とバージョン情報が表示

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

Rtools のダウンロードは無事に完了 】

3.【 次に、 R のパッケージ {inline}{Rcpp} をインストール 】

3-1. R を立ち上げて、 R のパッケージ {inline}{Rcpp}install.pacakges() で読み込み後、 library() でメモリに読み込む

RStudio
install.packages("inline")
install.packages("Rcpp")

library(Rcpp)
library(inline)

4.【 続いて、 {RStan} パッケージのインストール 】

(参考) (リンク先)RStan githubページ

RStan.install.site.png

※ 上記の 準備工程 Prerequistics の項目に記載されている下記は、すでに導入済み

(環境変数 PATH の設定も終わっている)

  • R 本体
  • Rtools

以下の記載指示に従い、 RStan をインストールする

RStan.install.site.3.png

RStudio コンソールに、下記を入力して実行

RStudio
Sys.setenv(MAKEFLAGS = "-j4") 

source('http://mc-stan.org/rstan/install.R', echo = TRUE, max.deparse.length = 2000)

install_rstan() 

※長い表示が続いた後、以下が末尾に表示される。

RStudio
 g++ -m64 -I"C:/PROGRA~1/R/R-31~1.0/include" -DNDEBUG -I"../inst/include/stansrc"   -I"../inst/include/stanlib/eigen_3.2.0"  -I"../inst/include/stanlib/eigen_3.2.0/unsupported"  -I"../inst/include/stanlib/boost_1.54.0" -I"../inst/include"  -DBOOST_RESULT_OF_USE_TR1 -DBOOST_NO_DECLTYPE -DBOOST_DISABLE_ASSERTS   -I"c:/usr/local/home/R/win-library/3.1/Rcpp/include" -I"d:/RCompile/CRANpkg/extralibs64/local/include"     -O3  -c gm__grammars__whitespace_grammar_inst.cpp -o gm__grammars__whitespace_grammar_inst.o
g++ -m64 -I"C:/PROGRA~1/R/R-31~1.0/include" -DNDEBUG -I"../inst/include/stansrc"   -I"../inst/include/stanlib/eigen_3.2.0"  -I"../inst/include/stanlib/eigen_3.2.0/unsupported"  -I"../inst/include/stanlib/boost_1.54.0" -I"../inst/include"  -DBOOST_RESULT_OF_USE_TR1 -DBOOST_NO_DECLTYPE -DBOOST_DISABLE_ASSERTS   -I"c:/usr/local/home/R/win-library/3.1/Rcpp/include" -I"d:/RCompile/CRANpkg/extralibs64/local/include"     -O3  -c agrad__rev__var_stack.cpp -o agrad__rev__var_stack.o
ar qc libstan.a agrad__rev__var_stack.o 
# -@if test ! -e ../inst/libstan/x64; then mkdir -p ../inst/libstan/x64; fi
# cp libstan.a ../inst/libstan/x64
# rm libstan.a 
g++ -m64 -shared -s -static-libgcc -o rstan.dll tmp.def stanc.o chains.o misc.o init.o gm__ast_def.o gm__grammars__bare_type_grammar_inst.o gm__grammars__expression07_grammar_inst.o gm__grammars__expression_grammar_inst.o gm__grammars__functions_grammar_inst.o gm__grammars__program_grammar_inst.o gm__grammars__statement_2_grammar_inst.o gm__grammars__statement_grammar_inst.o gm__grammars__term_grammar_inst.o gm__grammars__var_deccls_grammar_inst.o gm__grammars__whitespace_grammar_inst.o -Ld:/RCompile/CRANpkg/extralibs64/local/lib/x64 -Ld:/RCompile/CRANpkg/extralibs64/local/lib -LC:/PROGRA~1/R/R-31~1.0/bin/x64 -lR
installing via 'install.libs.R' to c:/usr/local/home/R/win-library/3.1/rstan
installing rstan libs to c:/usr/local/home/R/win-library/3.1/rstan/libs/x64
installing libstan.a to c:/usr/local/home/R/win-library/3.1/rstan/libstan/x64
** R
** inst
** preparing package for lazy loading
** help
*** installing help indices
** building package indices
** installing vignettes
** testing if installed package can be loaded
* DONE (rstan)

The downloaded source packages are in
    ‘C:\Users\HirofumiYashima\AppData\Local\Temp\Rtmpo59GaD\downloaded_packages’
rstan (Version 2.3.0, packaged: 2014-06-23 19:44:07 UTC, GitRev: 5dd1fe5a323a)
mode fast without NDEBUG for compiling C++ code has been set already

{RStan} パッケージのインストールが完了したので、 R のメモリに読み込む

RStudio
library(rstan) 

(実行結果)

RStudio
rstan (Version 2.3.0, packaged: 2014-06-23 19:44:07 UTC, GitRev: 5dd1fe5a323a)

5. Rstan を使ってみる

(参考) (リンク先)RStan githubページ内 "Example 1: Eight Schools"項目

RStan.install.site.4.examples.png

以下を、R (or RStudio)のコンソールに入力

RStudio
schools_code <- '
  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);
  }
'

RStudio

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

RStudio
fit <- stan(model_code = schools_code, data = schools_dat, 
            iter = 1000, chains = 4)

【実行結果】

(以下が表示されたあと、数十秒 待つ)

RStudio
TRANSLATING MODEL 'schools_code' FROM Stan CODE TO C++ CODE NOW.
COMPILING THE C++ CODE FOR MODEL 'schools_code' NOW.

(すると、以下のように、実行過程が表示される)

RStudio
cygwin warning:
  MS-DOS style path detected: C:/PROGRA~1/R/R-31~1.0/etc/x64/Makeconf
  Preferred POSIX equivalent is: /cygdrive/c/PROGRA~1/R/R-31~1.0/etc/x64/Makeconf
  CYGWIN environment variable option "nodosfilewarning" turns off this warning.
  Consult the user's guide for more details about POSIX paths:
    http://cygwin.com/cygwin-ug-net/using.html#using-pathnames
c:/usr/local/home/R/win-library/3.1/rstan/include//stansrc/stan/agrad/rev/var_stack.hpp:49:17: warning: 'void stan::agrad::free_memory()' defined but not used [-Wunused-function]

SAMPLING FOR MODEL 'schools_code' NOW (CHAIN 1).

Iteration:   1 / 1000 [  0%]  (Warmup)
Iteration: 100 / 1000 [ 10%]  (Warmup)
Iteration: 200 / 1000 [ 20%]  (Warmup)
Iteration: 300 / 1000 [ 30%]  (Warmup)
Iteration: 400 / 1000 [ 40%]  (Warmup)
Iteration: 500 / 1000 [ 50%]  (Warmup)
Iteration: 501 / 1000 [ 50%]  (Sampling)
Iteration: 600 / 1000 [ 60%]  (Sampling)
Iteration: 700 / 1000 [ 70%]  (Sampling)
Iteration: 800 / 1000 [ 80%]  (Sampling)
Iteration: 900 / 1000 [ 90%]  (Sampling)
Iteration: 1000 / 1000 [100%]  (Sampling)
#  Elapsed Time: 0.016 seconds (Warm-up)
#                0.016 seconds (Sampling)
#                0.032 seconds (Total)


SAMPLING FOR MODEL 'schools_code' NOW (CHAIN 2).

Iteration:   1 / 1000 [  0%]  (Warmup)
Iteration: 100 / 1000 [ 10%]  (Warmup)
Iteration: 200 / 1000 [ 20%]  (Warmup)
Iteration: 300 / 1000 [ 30%]  (Warmup)
Iteration: 400 / 1000 [ 40%]  (Warmup)
Iteration: 500 / 1000 [ 50%]  (Warmup)
Iteration: 501 / 1000 [ 50%]  (Sampling)
Iteration: 600 / 1000 [ 60%]  (Sampling)
Iteration: 700 / 1000 [ 70%]  (Sampling)
Iteration: 800 / 1000 [ 80%]  (Sampling)
Iteration: 900 / 1000 [ 90%]  (Sampling)
Iteration: 1000 / 1000 [100%]  (Sampling)
#  Elapsed Time: 0.015 seconds (Warm-up)
#                0.013 seconds (Sampling)
#                0.028 seconds (Total)


SAMPLING FOR MODEL 'schools_code' NOW (CHAIN 3).

Iteration:   1 / 1000 [  0%]  (Warmup)
Iteration: 100 / 1000 [ 10%]  (Warmup)
Iteration: 200 / 1000 [ 20%]  (Warmup)
Iteration: 300 / 1000 [ 30%]  (Warmup)
Iteration: 400 / 1000 [ 40%]  (Warmup)
Iteration: 500 / 1000 [ 50%]  (Warmup)
Iteration: 501 / 1000 [ 50%]  (Sampling)
Iteration: 600 / 1000 [ 60%]  (Sampling)
Iteration: 700 / 1000 [ 70%]  (Sampling)
Iteration: 800 / 1000 [ 80%]  (Sampling)
Iteration: 900 / 1000 [ 90%]  (Sampling)
Iteration: 1000 / 1000 [100%]  (Sampling)
#  Elapsed Time: 0.015 seconds (Warm-up)
#                0.012 seconds (Sampling)
#                0.027 seconds (Total)


SAMPLING FOR MODEL 'schools_code' NOW (CHAIN 4).

Iteration:   1 / 1000 [  0%]  (Warmup)
Iteration: 100 / 1000 [ 10%]  (Warmup)
Iteration: 200 / 1000 [ 20%]  (Warmup)
Iteration: 300 / 1000 [ 30%]  (Warmup)
Iteration: 400 / 1000 [ 40%]  (Warmup)
Iteration: 500 / 1000 [ 50%]  (Warmup)
Iteration: 501 / 1000 [ 50%]  (Sampling)
Iteration: 600 / 1000 [ 60%]  (Sampling)
Iteration: 700 / 1000 [ 70%]  (Sampling)
Iteration: 800 / 1000 [ 80%]  (Sampling)
Iteration: 900 / 1000 [ 90%]  (Sampling)
Iteration: 1000 / 1000 [100%]  (Sampling)
#  Elapsed Time: 0.017 seconds (Warm-up)
#                0.011 seconds (Sampling)
#                0.028 seconds (Total)

> 

実行結果を表示させる

RStudio
print(fit);

実行結果が表示された

RStudio
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        8.00    0.16 4.88  -1.24  4.81  7.84 11.11 17.77   904 1.00
tau       6.14    0.20 4.93   0.19  2.37  5.10  8.69 18.71   625 1.01
eta[1]    0.41    0.03 0.94  -1.42 -0.22  0.39  1.05  2.24  1263 1.00
eta[2]   -0.02    0.02 0.90  -1.79 -0.63 -0.02  0.60  1.71  1426 1.00
eta[3]   -0.21    0.02 0.90  -1.98 -0.80 -0.24  0.39  1.62  1351 1.00
eta[4]   -0.03    0.02 0.90  -1.78 -0.63 -0.05  0.54  1.75  1506 1.00
eta[5]   -0.34    0.02 0.87  -2.05 -0.93 -0.37  0.24  1.39  1442 1.00
eta[6]   -0.20    0.02 0.95  -1.98 -0.81 -0.27  0.43  1.71  1578 1.00
eta[7]    0.33    0.02 0.86  -1.41 -0.22  0.35  0.91  1.92  1421 1.00
eta[8]    0.05    0.02 0.93  -1.86 -0.57  0.04  0.69  1.84  1525 1.00
theta[1] 11.28    0.21 7.92  -1.58  6.11 10.25 15.30 31.04  1364 1.00
theta[2]  7.88    0.14 6.25  -3.70  3.83  7.84 11.84 20.74  1932 1.00
theta[3]  6.13    0.20 7.50 -10.91  2.12  6.48 10.64 19.97  1365 1.00
theta[4]  7.55    0.16 6.45  -4.76  3.53  7.47 11.44 21.30  1645 1.00
theta[5]  5.41    0.15 6.07  -7.74  1.77  5.94  9.32 16.60  1545 1.00
theta[6]  6.35    0.16 6.49  -7.57  2.49  6.63 10.32 18.90  1597 1.00
theta[7] 10.47    0.18 6.77  -1.28  6.02  9.81 14.31 25.08  1487 1.00
theta[8]  8.57    0.21 7.47  -5.54  4.12  8.37 12.67 24.10  1230 1.00
lp__     -4.95    0.11 2.53 -10.76 -6.45 -4.72 -3.23 -0.51   566 1.01

Samples were drawn using NUTS(diag_e) at Mon Jun 30 16:36:51 2014.
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).

次に、以下の例を実行してみる。(※上記、RStan 公式ウェブサイトより)

RStan.install.site.5.examples.2.png

In RStan, we can also specify the Stan model using a file. For example, we can download 8schools.stan file to our working directory and use the following usage of function stan instead.

※ 公式ウェブサイトから、 8schools.stan ファイルをダウンロードして、以下を実行

※ 先ほど、 schools_code 変数に格納していた Stan スクリプトを、外部ファイル (ファイル名 8schools.stan )に保存して、 R から呼び出す方式を実行してみます

8schools.stan の青文字にマウスポインタを合わせて以下を実行
マウス右ボタンをクリック→「名前をつけてリンク先を保存」をクリック 
※すると、保存先のファイルを選択する画面が開く。ファイルの種類は、自動的に「STANファイル(.stan)」になっている
※ setwd()で、上記ファイルの格納先フォルダに、Rの作業フォルダをあわせる
⇒ list.files()で、上記ファイルが表示されるか確認

下記を実行

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

(実行過程の表示)

RStudio
TRANSLATING MODEL '8schools' FROM Stan CODE TO C++ CODE NOW.
COMPILING THE C++ CODE FOR MODEL '8schools' NOW.
cygwin warning:
  MS-DOS style path detected: C:/PROGRA~1/R/R-31~1.0/etc/x64/Makeconf
  Preferred POSIX equivalent is: /cygdrive/c/PROGRA~1/R/R-31~1.0/etc/x64/Makeconf
  CYGWIN environment variable option "nodosfilewarning" turns off this warning.
  Consult the user's guide for more details about POSIX paths:
    http://cygwin.com/cygwin-ug-net/using.html#using-pathnames
c:/usr/local/home/R/win-library/3.1/rstan/include//stansrc/stan/agrad/rev/var_stack.hpp:49:17: warning: 'void stan::agrad::free_memory()' defined but not used [-Wunused-function]

SAMPLING FOR MODEL '8schools' NOW (CHAIN 1).

Iteration:   1 / 1000 [  0%]  (Warmup)
Iteration: 100 / 1000 [ 10%]  (Warmup)
Iteration: 200 / 1000 [ 20%]  (Warmup)
Iteration: 300 / 1000 [ 30%]  (Warmup)
Iteration: 400 / 1000 [ 40%]  (Warmup)
Iteration: 500 / 1000 [ 50%]  (Warmup)
Iteration: 501 / 1000 [ 50%]  (Sampling)
Iteration: 600 / 1000 [ 60%]  (Sampling)
Iteration: 700 / 1000 [ 70%]  (Sampling)
Iteration: 800 / 1000 [ 80%]  (Sampling)
Iteration: 900 / 1000 [ 90%]  (Sampling)
Iteration: 1000 / 1000 [100%]  (Sampling)
#  Elapsed Time: 0.017 seconds (Warm-up)
#                0.012 seconds (Sampling)
#                0.029 seconds (Total)


SAMPLING FOR MODEL '8schools' NOW (CHAIN 2).

Iteration:   1 / 1000 [  0%]  (Warmup)
Iteration: 100 / 1000 [ 10%]  (Warmup)
Iteration: 200 / 1000 [ 20%]  (Warmup)
Iteration: 300 / 1000 [ 30%]  (Warmup)
Iteration: 400 / 1000 [ 40%]  (Warmup)
Iteration: 500 / 1000 [ 50%]  (Warmup)
Iteration: 501 / 1000 [ 50%]  (Sampling)
Iteration: 600 / 1000 [ 60%]  (Sampling)
Iteration: 700 / 1000 [ 70%]  (Sampling)
Iteration: 800 / 1000 [ 80%]  (Sampling)
Iteration: 900 / 1000 [ 90%]  (Sampling)
Iteration: 1000 / 1000 [100%]  (Sampling)
#  Elapsed Time: 0.019 seconds (Warm-up)
#                0.015 seconds (Sampling)
#                0.034 seconds (Total)


SAMPLING FOR MODEL '8schools' NOW (CHAIN 3).

Iteration:   1 / 1000 [  0%]  (Warmup)
Iteration: 100 / 1000 [ 10%]  (Warmup)
Iteration: 200 / 1000 [ 20%]  (Warmup)
Iteration: 300 / 1000 [ 30%]  (Warmup)
Iteration: 400 / 1000 [ 40%]  (Warmup)
Iteration: 500 / 1000 [ 50%]  (Warmup)
Iteration: 501 / 1000 [ 50%]  (Sampling)
Iteration: 600 / 1000 [ 60%]  (Sampling)
Iteration: 700 / 1000 [ 70%]  (Sampling)
Iteration: 800 / 1000 [ 80%]  (Sampling)
Iteration: 900 / 1000 [ 90%]  (Sampling)
Iteration: 1000 / 1000 [100%]  (Sampling)
#  Elapsed Time: 0.015 seconds (Warm-up)
#                0.014 seconds (Sampling)
#                0.029 seconds (Total)


SAMPLING FOR MODEL '8schools' NOW (CHAIN 4).

Iteration:   1 / 1000 [  0%]  (Warmup)
Iteration: 100 / 1000 [ 10%]  (Warmup)
Iteration: 200 / 1000 [ 20%]  (Warmup)
Iteration: 300 / 1000 [ 30%]  (Warmup)
Iteration: 400 / 1000 [ 40%]  (Warmup)
Iteration: 500 / 1000 [ 50%]  (Warmup)
Iteration: 501 / 1000 [ 50%]  (Sampling)
Iteration: 600 / 1000 [ 60%]  (Sampling)
Iteration: 700 / 1000 [ 70%]  (Sampling)
Iteration: 800 / 1000 [ 80%]  (Sampling)
Iteration: 900 / 1000 [ 90%]  (Sampling)
Iteration: 1000 / 1000 [100%]  (Sampling)
#  Elapsed Time: 0.018 seconds (Warm-up)
#                0.019 seconds (Sampling)
#                0.037 seconds (Total)

(実行が完了した)

(実行結果を表示)

RStudio
print(fit1);

(実行結果)

RStudio
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        7.95    0.23 4.98  -1.61  4.79  7.86 11.10 17.83   455 1.01
tau       6.48    0.32 5.76   0.27  2.34  5.05  8.85 19.88   328 1.01
eta[1]    0.38    0.03 0.97  -1.61 -0.25  0.39  1.04  2.24  1152 1.00
eta[2]    0.00    0.03 0.87  -1.67 -0.57 -0.02  0.58  1.77   923 1.00
eta[3]   -0.21    0.03 0.93  -2.01 -0.85 -0.22  0.41  1.69  1313 1.00
eta[4]   -0.03    0.02 0.90  -1.87 -0.62 -0.01  0.58  1.72  1462 1.00
eta[5]   -0.31    0.02 0.88  -1.92 -0.87 -0.34  0.22  1.54  1413 1.00
eta[6]   -0.18    0.02 0.87  -1.83 -0.77 -0.22  0.39  1.61  1456 1.00
eta[7]    0.35    0.03 0.90  -1.53 -0.24  0.38  0.95  2.05   991 1.00
eta[8]    0.05    0.03 1.00  -1.97 -0.62  0.06  0.72  1.97  1294 1.00
theta[1] 11.67    0.60 9.86  -2.52  5.90 10.09 15.47 33.96   274 1.01
theta[2]  7.93    0.16 6.16  -4.79  4.12  7.89 11.79 20.70  1487 1.00
theta[3]  6.01    0.22 7.63 -11.99  2.31  6.67 10.80 19.12  1186 1.00
theta[4]  7.63    0.16 6.48  -5.90  3.77  7.60 11.51 20.95  1555 1.00
theta[5]  5.23    0.15 6.19  -7.87  1.53  5.63  9.15 16.14  1618 1.00
theta[6]  6.27    0.18 6.68  -7.36  2.37  6.49 10.64 18.57  1418 1.00
theta[7] 10.75    0.19 6.80  -0.36  6.27  9.93 14.67 25.67  1305 1.00
theta[8]  8.35    0.22 7.87  -6.52  3.94  7.84 12.36 26.45  1240 1.00
lp__     -4.98    0.10 2.57 -10.64 -6.56 -4.72 -3.13 -0.75   693 1.00

Samples were drawn using NUTS(diag_e) at Mon Jun 30 16:37:52 2014.
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).

(再び、公式ウェブサイトに戻る)

RStan.install.site.5.examples.2.png

Once a model is fitted, we can use the fitted result as an input to fit the model with other data or settings. This would save us time of compiling the C++ code for the model. By specifying parameter fit for function stan, we can fit the model again. For example, if we want to sample more iterations, we can proceed as follows.

fit2 <- stan(fit = fit1, data = schools_dat, iter = 10000, chains = 4)

上記の記述に従い、
さきほど、生成したモデル「fit1」を入力引数にして、以下を実行する

RStudio
fit2 <- stan(fit = fit1, data = schools_dat, iter = 10000, chains = 4)

(実行過程)

RStudio
SAMPLING FOR MODEL '8schools' NOW (CHAIN 1).

Iteration:    1 / 10000 [  0%]  (Warmup)
Iteration: 1000 / 10000 [ 10%]  (Warmup)
Iteration: 2000 / 10000 [ 20%]  (Warmup)
Iteration: 3000 / 10000 [ 30%]  (Warmup)
Iteration: 4000 / 10000 [ 40%]  (Warmup)
Iteration: 5000 / 10000 [ 50%]  (Warmup)
Iteration: 5001 / 10000 [ 50%]  (Sampling)
Iteration: 6000 / 10000 [ 60%]  (Sampling)
Iteration: 7000 / 10000 [ 70%]  (Sampling)
Iteration: 8000 / 10000 [ 80%]  (Sampling)
Iteration: 9000 / 10000 [ 90%]  (Sampling)
Iteration: 10000 / 10000 [100%]  (Sampling)
#  Elapsed Time: 0.162 seconds (Warm-up)
#                0.14 seconds (Sampling)
#                0.302 seconds (Total)


SAMPLING FOR MODEL '8schools' NOW (CHAIN 2).

Iteration:    1 / 10000 [  0%]  (Warmup)
Iteration: 1000 / 10000 [ 10%]  (Warmup)
Iteration: 2000 / 10000 [ 20%]  (Warmup)
Iteration: 3000 / 10000 [ 30%]  (Warmup)
Iteration: 4000 / 10000 [ 40%]  (Warmup)
Iteration: 5000 / 10000 [ 50%]  (Warmup)
Iteration: 5001 / 10000 [ 50%]  (Sampling)
Iteration: 6000 / 10000 [ 60%]  (Sampling)
Iteration: 7000 / 10000 [ 70%]  (Sampling)
Iteration: 8000 / 10000 [ 80%]  (Sampling)
Iteration: 9000 / 10000 [ 90%]  (Sampling)
Iteration: 10000 / 10000 [100%]  (Sampling)
#  Elapsed Time: 0.132 seconds (Warm-up)
#                0.141 seconds (Sampling)
#                0.273 seconds (Total)


SAMPLING FOR MODEL '8schools' NOW (CHAIN 3).

Iteration:    1 / 10000 [  0%]  (Warmup)
Iteration: 1000 / 10000 [ 10%]  (Warmup)
Iteration: 2000 / 10000 [ 20%]  (Warmup)
Iteration: 3000 / 10000 [ 30%]  (Warmup)
Iteration: 4000 / 10000 [ 40%]  (Warmup)
Iteration: 5000 / 10000 [ 50%]  (Warmup)
Iteration: 5001 / 10000 [ 50%]  (Sampling)
Iteration: 6000 / 10000 [ 60%]  (Sampling)
Iteration: 7000 / 10000 [ 70%]  (Sampling)
Iteration: 8000 / 10000 [ 80%]  (Sampling)
Iteration: 9000 / 10000 [ 90%]  (Sampling)
Iteration: 10000 / 10000 [100%]  (Sampling)
#  Elapsed Time: 0.131 seconds (Warm-up)
#                0.162 seconds (Sampling)
#                0.293 seconds (Total)


SAMPLING FOR MODEL '8schools' NOW (CHAIN 4).

Iteration:    1 / 10000 [  0%]  (Warmup)
Iteration: 1000 / 10000 [ 10%]  (Warmup)
Iteration: 2000 / 10000 [ 20%]  (Warmup)
Iteration: 3000 / 10000 [ 30%]  (Warmup)
Iteration: 4000 / 10000 [ 40%]  (Warmup)
Iteration: 5000 / 10000 [ 50%]  (Warmup)
Iteration: 5001 / 10000 [ 50%]  (Sampling)
Iteration: 6000 / 10000 [ 60%]  (Sampling)
Iteration: 7000 / 10000 [ 70%]  (Sampling)
Iteration: 8000 / 10000 [ 80%]  (Sampling)
Iteration: 9000 / 10000 [ 90%]  (Sampling)
Iteration: 10000 / 10000 [100%]  (Sampling)
#  Elapsed Time: 0.14 seconds (Warm-up)
#                0.114 seconds (Sampling)
#                0.254 seconds (Total)

ウェブサイトの下記の記載に従い、結果を表示させる

RStan.install.site.5.examples.2.png

The object fit2, returned from function stan is an S4 object of class stanfit. Methods such as print and plot are associated with the fitted result so we can use the following code to check out the results in fit2. print provides a summary for the parameter of the model as well as the log-posterior with name lp__ (see the following example output).
For more methods and details of class stanfit, see the help of class stanfit.

In particular, we can use extract function on stanfit objects to obtain the samples. extract extracts samples from the stanfit object as a list of arrays for parameters of interest, or just an array. In addition, S3 functions as.array and as.matrix are defined for stanfit object (using help("as.array.stanfit") to check out the help document in R).

(実行結果を表示させる)

RStudio
print(fit2)

(実行結果)

RStudio
Inference for Stan model: 8schools.
4 chains, each with iter=10000; warmup=5000; thin=1; 
post-warmup draws per chain=5000, total post-warmup draws=20000.

          mean se_mean   sd   2.5%   25%   50%   75% 97.5% n_eff Rhat
mu        8.01    0.07 5.09  -1.84  4.67  7.96 11.25 18.27  5540    1
tau       6.60    0.08 5.49   0.25  2.49  5.34  9.21 20.77  5017    1
eta[1]    0.38    0.01 0.94  -1.52 -0.22  0.40  1.02  2.20 10731    1
eta[2]    0.00    0.01 0.87  -1.72 -0.57  0.00  0.56  1.72 12225    1
eta[3]   -0.21    0.01 0.93  -2.03 -0.82 -0.21  0.40  1.67 13288    1
eta[4]   -0.04    0.01 0.89  -1.81 -0.62 -0.05  0.54  1.74 12446    1
eta[5]   -0.37    0.01 0.87  -2.06 -0.94 -0.39  0.19  1.45 11499    1
eta[6]   -0.21    0.01 0.89  -1.95 -0.80 -0.22  0.36  1.56 13001    1
eta[7]    0.33    0.01 0.90  -1.49 -0.25  0.35  0.93  2.08 11393    1
eta[8]    0.06    0.01 0.94  -1.83 -0.56  0.07  0.69  1.88 12202    1
theta[1] 11.47    0.09 8.41  -1.91  5.94 10.22 15.58 31.76  7993    1
theta[2]  7.94    0.05 6.25  -4.74  4.05  7.87 11.87 20.34 15143    1
theta[3]  6.12    0.08 7.76 -11.40  2.10  6.66 10.84 20.37 10693    1
theta[4]  7.69    0.05 6.58  -5.52  3.65  7.75 11.74 21.02 14920    1
theta[5]  5.09    0.05 6.34  -8.95  1.28  5.56  9.36 16.38 13846    1
theta[6]  6.20    0.06 6.73  -8.44  2.29  6.54 10.46 18.81 13474    1
theta[7] 10.74    0.06 6.83  -1.07  6.12 10.13 14.72 26.12 11121    1
theta[8]  8.57    0.07 7.99  -6.97  3.92  8.26 12.89 25.89 11535    1
lp__     -4.87    0.04 2.66 -10.82 -6.46 -4.61 -3.01 -0.39  5191    1

Samples were drawn using NUTS(diag_e) at Mon Jun 30 16:40:45 2014.
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).

(実行結果の描画)

RStudio
plot(fit2)

graph.png

(以下を実行)

RStudio
print(fit, digits = 1)

(実行結果)

RStudio
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        8.0     0.2 4.9  -1.2  4.8  7.8 11.1  17.8   904    1
tau       6.1     0.2 4.9   0.2  2.4  5.1  8.7  18.7   625    1
eta[1]    0.4     0.0 0.9  -1.4 -0.2  0.4  1.0   2.2  1263    1
eta[2]    0.0     0.0 0.9  -1.8 -0.6  0.0  0.6   1.7  1426    1
eta[3]   -0.2     0.0 0.9  -2.0 -0.8 -0.2  0.4   1.6  1351    1
eta[4]    0.0     0.0 0.9  -1.8 -0.6  0.0  0.5   1.8  1506    1
eta[5]   -0.3     0.0 0.9  -2.1 -0.9 -0.4  0.2   1.4  1442    1
eta[6]   -0.2     0.0 0.9  -2.0 -0.8 -0.3  0.4   1.7  1578    1
eta[7]    0.3     0.0 0.9  -1.4 -0.2  0.4  0.9   1.9  1421    1
eta[8]    0.1     0.0 0.9  -1.9 -0.6  0.0  0.7   1.8  1525    1
theta[1] 11.3     0.2 7.9  -1.6  6.1 10.3 15.3  31.0  1364    1
theta[2]  7.9     0.1 6.2  -3.7  3.8  7.8 11.8  20.7  1932    1
theta[3]  6.1     0.2 7.5 -10.9  2.1  6.5 10.6  20.0  1365    1
theta[4]  7.5     0.2 6.4  -4.8  3.5  7.5 11.4  21.3  1645    1
theta[5]  5.4     0.2 6.1  -7.7  1.8  5.9  9.3  16.6  1545    1
theta[6]  6.4     0.2 6.5  -7.6  2.5  6.6 10.3  18.9  1597    1
theta[7] 10.5     0.2 6.8  -1.3  6.0  9.8 14.3  25.1  1487    1
theta[8]  8.6     0.2 7.5  -5.5  4.1  8.4 12.7  24.1  1230    1
lp__     -5.0     0.1 2.5 -10.8 -6.4 -4.7 -3.2  -0.5   566    1

Samples were drawn using NUTS(diag_e) at Mon Jun 30 16:36:51 2014.
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).

RStan は、正常に動作した!) :laughing:  

【備考: (参考) MCMC Sampler RStan を使ったデータ解析】

Stan ファイルの記述文法や、 R から RStan 経由で Stan を動かして、ロジスティック回帰分析を行った事例は、下記が参考になります

(リンク)銀座で働くデータサイエンティストのブログ 「MCMCの計算にStanを使ってみた(超基礎・導入編)」

Ref.ginza.1.png

Ref.ginza.2.png

Ref.ginza.3.png

Ref.ginza.4.png

Ref.ginza.5.png