Rでstan使ってみる~導入編 for Mac~ - 医療職からデータサイエンティストへ

新年になったということで、今年はRでstanを使ったベイズモデリングに挑戦していきたいと思います。

今回は導入編ということで、こちら

RStan Getting Started · stan-dev/rstan Wiki · GitHub

を参考にしながら進めていきます。ちなみに私のPCはMacですので、Windowsユーザーの方はあしからず。

www.medi-08-data-06.work

環境のチェック

まずはstanの推奨環境の確認をしていきます。上記サイトを見てみると

R vesion 3.4.0以上
Rstudio version 1.2以上

以上が推奨環境のようです。私もRstudioのバージョンが1.2未満だったため最新バージョンにしました。ちなみにRstudioのバージョンはRstudio上部のAbout Rstudioで、Rのバージョンは

sessionInfo()

で確認できます。

Rや Rstudioの公式サイトはこちら

The Comprehensive R Archive Network

RStudio Preview - RStudio

さっそくインストール

まずは、念の為こちらのコードを実行後、リスタートします。

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

そしてrstanをインストールしていきます。

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

少し時間がかかります,,,,
ここで何か予期せぬことが起こった場合は、こちらを参考にしてください。

Installing RStan from source on a Mac · stan-dev/rstan Wiki · GitHub

さて、rstanをインストールしたら、こちらのコードを実行します。

> pkgbuild::has_build_tools(debug = TRUE)
[1] TRUE

これは、stanを動かしているC++言語がうまく利用できるかどうかみたいな確認ですかね....。

この辺りは、詳しくありませんがTrueが返ってこればOKみたいなので、大丈夫みたいです。

高速化するための設定

こちらは推奨設定ということでやらなくてもいいみたいなのですが、高速化のために設定しておきます。

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("\nCXX14FLAGS=-O3 -march=native -mtune=native",
    if( grepl("^darwin", R.version$os)) "CXX14FLAGS += -arch x86_64 -ftemplate-depth-256" else 
      if (.Platform$OS.type == "windows") "CXX11FLAGS=-O3 -march=native -mtune=native" else
        "CXX14FLAGS += -fPIC",
    file = M, sep = "\n", append = TRUE)

以上のコードをコピペするだけで OKです。やっていることは、HOMEディレクトリに.Rディレクトリ（隠しフォルダ）を作成して、その中にMakevarsというファイルを作成しています。

実際に作成できているかどうかは、Teminalから確認できます。

$ cat ~/.R/Makevars
CXX14FLAGS=-O3 -march=native -mtune=native
CXX14FLAGS += -arch x86_64 -ftemplate-depth-256

使ってみる

stanの使い方はおいおい学んでいくとして、実際に動作確認していきます。まずは、インストールしたstanを読み込みます。

> library("rstan")
 要求されたパッケージ ggplot2 をロード中です 
 要求されたパッケージ StanHeaders をロード中です 
rstan (Version 2.18.2, GitRev: 2e1f913d3ca3)
For execution on a local, multicore CPU with excess RAM we recommend calling
options(mc.cores = parallel::detectCores()).
To avoid recompilation of unchanged Stan programs, we recommend calling
rstan_options(auto_write = TRUE)

ここで、メッセージを読んでみるとPCがマルチコアマシンで十分なRAMがあるときは、以下を実行してくださいとある。

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

私のmacは高スペックではないですが、どうやら並列処理を行う設定らしいので、実行しておきます。

次のメッセージ、変更されていないstanモデルの再コンパイルを防ぐためには、以下のコードを実行することをお勧めします。とある。

モデルの自動保存のことらしいです。これもやっておきましょう。

rstan_options(auto_write = TRUE)

これで準備が整いました！

次に好きなテキストエディタ、もしくはRstudioのFile -> New File -> Stan Fileで新しいファイルを開いたら以下をコピペして、8schools.stanとしてファイルを保存します。

保存場所は、後から読み込むためgetwd()で確認できるRの作業ディレクトリにしておくと便利です。

よくわかりませんがとりあえずは動作確認ということで、指示通りに実行していきます。

data {
  int<lower=0> J;         // number of schools 
  real y[J];              // estimated treatment effects
  real<lower=0> sigma[J]; // standard error of effect estimates 
}
parameters {
  real mu;                // population treatment effect
  real<lower=0> tau;      // standard deviation in treatment effects
  vector[J] eta;          // unscaled deviation from mu by school
}
transformed parameters {
  vector[J] theta = mu + tau * eta;        // school treatment effects
}
model {
  target += normal_lpdf(eta | 0, 1);       // prior log-density
  target += normal_lpdf(y | theta, sigma); // log-likelihood
}

上記ファイルが保存できたら、今度は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))
#stanを実行
fit <- stan(file = '8schools.stan', data = schools_dat)

#少し時間がかかる...
#プロットしてみる
plot(fit)