はじめに

ベイズ統計モデリングを行う確率的プログラミング言語にStanというものがあります。これをRで使うためのインターフェースとしてRStanがあります。これらを使うことで、統計モデリングにおいて、確率モデルをプログラミングコードで記述することで必要なパラメータの計算をほぼ自動で行うことができます。
私はベイズ統計モデリングの初心者なのもあり、「RStan: the R interface to Stan」by Stan Development Team(2019-07-09)」の内容を他の教科書もみながら試してみました。ということで、そのテキトーな訳も含めたメモになります。（ちょっと時間が無くなってしまい、中途半端に終わってしまったため、後日追記するか別記事に続きを書くなどしたいと思います。）

環境準備

Stanコードを動かすのにC++コンパイラが必要です。RStan Getting Startedを参考にインストールします。Macの場合はxcodeを入れれば十分のようです。
また、rstanはいくつか依存するパッケージがあるようですが、install.package("rstan', dependencies=TRUE)で入るので、特に意識する必要ありません。

実施フローについて

RStanを使ってベイズ統計モデリングを実施するにあたり、次のフローのように対応することが一般的とのことです。なお、下記の2から4はstan関数の呼び出しで内部的に連続的に実行するので、利用者は意識する必要は無いようです。

Stanを使って対数事後密度を書いて統計モデリングを表す。拡張子に.stanにてファイルを別途保存する。
stancでStanコードをC++コードに変換する。
変換したC++コードをRから呼び出し可能なDSO（DLLとも呼ぶ）にコンパイルする。
DSOを実行して事後分布からサンプリングする。
MCMCチェーンの収束状況を確認する。非収束している場合は対処する。
事後サンプルより、推定結果を解釈する。

Eight Schoolsの例をやってみる

Bayesian Data Analysis(Gelman et al. 2003)の5.5節にあるらしいEight Schoolsの例を紹介されているので、やってみます。
階層モデルを使用して、大学の入学試験に対するコーチングプログラムの効果をモデル化します。次の表にて、8つの高校での実験で得られた結果である推定値（平均スコア変化）と標準誤差（平均スコア変化の標準偏差）がデータとして得られているとします。

学校名	推定値($y_i$)	標準偏差 ($\sigma_i$)
A	28	15
B	8	1
C	-3	16
D	7	11
E	-1	9
F	1	11
G	18	10
H	12	18

それぞれのパラメータは次の確率分布で表わされるものと仮定します。

\begin{align}
&y_j \sim 𝖭𝗈𝗋𝗆𝖺𝗅(\theta_j, \sigma_j), \quad j=1\ldots8 \\
&\theta_j \sim 𝖭𝗈𝗋𝗆𝖺𝗅(\mu, \tau), \quad j=1\ldots8 \\
&p(\mu, \tau) \propto 1
\end{align}

Normal(μ, σ)は平均がμで標準偏差がσの正規分布を示します。今回、推定したいパラメータθはさらに別のパラメータを持つ確率分布より生成されるモデルをおいており、階層的になっていることに注意してください。このようなモデルを階層ベイズモデルと言います。これをStanで表現してみます。

Stanでのプログラミング

次のようにStanは、「dataブロック(データの宣言)」、「parameterブロック（サンプリングしたいパラメータの宣言）」、「modelブロック（尤度と事前分布の記述）」の３つで基本的には構成されます。下記例にはさらにtransformed parametersブロックもありますが、これはparameterブロックで宣言しているパラメータを使って変換を施すような別のパラメータを宣言するときに使います。

schools.stan

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;
  vector[J] eta;
}
transformed parameters {
  vector[J] theta;
  theta = mu + tau * eta;
}
model {
  target += normal_lpdf(eta | 0, 1);
  target += normal_lpdf(y | theta, sigma);
}

dataブロックは、先の表に示したデータをR側で代入できるように宣言しておきます。parameterブロックとtransformed parameterブロックで、modelブロックで使用するパラメータを宣言しておき、最後にmodelブロックで仮定する確率密度関数を定義しています。ちなみに、最後を空行で終わらせる必要がありますので、ご注意を。

R側でのプログラミングとモデルの確認

Stan側でモデルを定義できましたので、R側で入力するデータを作ってモデルを作成し、モデルの確認を行います。

データの作成とMCMCサンプリング

まず、以下で、データの作成とMCMCサンプリングの実施を行い、その結果を保存します。下記ではデータを直で入れていますが、普通にread.csv等でデータファイルを別途作成し取り込んでも問題ありません。

ex_stan.R

library(rstan)

#dataブロックで宣言した変数に入れるデータを作成
schools_data <- list(
  J = 8,
  y = c(28,  8, -3,  7, -1,  1, 18, 12),
  sigma = c(15, 10, 16, 11,  9, 11, 10, 18)
)

#ベイズ統計モデリングの実施
fit1 <- stan(
  file = "schools.stan",  # Stan program
  data = schools_data,    # named list of data
  chains = 4,             # number of Markov chains
  warmup = 1000,          # number of warmup iterations per chain
  iter = 2000,            # total number of iterations per chain
  cores = 2,              # number of cores (could use one per chain)
  refresh = 0             # no progress shown
)

パラメータ収束の確認

次にパラメータの収束を確認します。次のようにモデルの中身の情報を確認します。

ex_stan.R(続き)

#パラメータの収束を確認する
print(fit1, pars=c("theta", "mu", "tau", "lp__"), probs=c(.1,.5,.9))

これを実行することで、次が出力されます。

Inference for Stan model: schools.
4 chains, each with iter=2000; warmup=1000; thin=1; 
post-warmup draws per chain=1000, total post-warmup draws=4000.

           mean se_mean   sd    10%    50%    90% n_eff Rhat
theta[1]  11.29    0.17 8.49   1.87  10.10  22.44  2544    1
theta[2]   7.69    0.09 6.29  -0.16   7.71  15.48  4847    1
theta[3]   6.07    0.15 7.57  -3.13   6.53  14.47  2628    1
theta[4]   7.49    0.11 6.56  -0.46   7.46  15.39  3871    1
theta[5]   5.14    0.12 6.32  -3.10   5.52  12.73  2814    1
theta[6]   6.03    0.11 6.66  -2.58   6.51  13.96  3702    1
theta[7]  10.51    0.11 6.65   2.93   9.85  19.25  3375    1
theta[8]   8.43    0.15 7.85  -0.10   8.14  17.81  2915    1
mu         7.69    0.13 5.01   1.61   7.71  13.98  1597    1
tau        6.39    0.15 5.35   0.90   5.19  13.15  1323    1
lp__     -39.63    0.08 2.72 -43.23 -39.34 -36.40  1202    1

Samples were drawn using NUTS(diag_e) at Wed Dec 25 22:20:51 2019.
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).

StanとRでベイズ統計モデリングによると、「chain数3以上かつ全パラメータのRhat<1.1で収束したとみなす」とのことなので、ひとまずはこれが守られているかを必要があります。この収束を確認するまで、次のステップに行くことはやめたほうが良さそうです。

traceplotを確認してみる

次にtraceplotを確認してみます。これは各MCMCチェーンの軌跡を表示させるもので、今回は4chainで実行していますので、4つ分の線グラフが見られます。ここで確認すべきは各chainで同じような軌跡になっているかです。そうなっていると収束していると判断できます。
次のように実行します。

ex_stan.R(続き)

traceplot(fit1, pars = c("mu", "tau"), inc_warmup = TRUE, nrow = 2)

すると、次の図が生成されます。

今回の場合、μもτも各chainが同じようにプロットされているのが分かるので、問題ないようです。

他の収束診断方法

ggmcmcでも収束診断のPDFファイルを出力することができます。何と次のように生成したモデルオブジェクトと出力ファイル名を指定するだけで、traceplotなど必要な図が載ったPDFが生成されます。

library(ggmcmc)
ggmcmc(ggs(fit1), file='ggmcmc_check.pdf')

詳しい説明はこちらUsing the ggmcmc packageを確認してください。
また、最近出たらしいbayesplotも良いらしいのですが、使い方がきちんと分かってないので割愛します。

最後に

ひとまず、チュートリアルに沿ってモデルの設計とモデル収束確認までさっとやってみました。また、収束確認については、ggmcmcも使ってみました。

アドカレ2019の最後に

今回は初の会社でのアドカレだったのですが、何とか形になりました。書いてくださった方々、お疲れ様でした。また、見てくださった方々、どうもありがとうございました。

実行環境

MacOS 10.14.6
R 3.6.2
RStudio 1.2.5033
rstan 2.19.2
ggmcmc 1.3

「RStan: the R interface to Stan」のチュートリアルを試してみる（収束確認まで）