RStan: the R interface to Stan

Stan Development Team

2022-11-26

この vignette では、Stanに対するRのインターフェースであるRStanを紹介する。StanはNo-U-Turnサンプラー（ハミルトンモンテカルロの変種）を用いたベイズ推論や、最適化による頻度論的推論のためのC++ライブラリである。GelmanCarlinSternRubin:2003の例を通して、RStanのフィーチャを説明する。

イントロダクション

Stanはベイズモデリングと推論のためのC++ライブラリで、主にNo-U-Turn sampler (NUTS) (Hoffman and Gelman 2012)、ユーザが指定したモデルとデータが与えられたときの事後シミュレーションを得るために使用する。また、LBFGS最適化アルゴリズムを利用して、対数尤度のような目的関数を最大化することも可能である。Rパッケージの__rstan__は、StanへのRインターフェースであるRStanを提供する。rstan__パッケージは、R (R Core Team 2014) からStanモデルを簡便に適合させ、事後推論や対数事後密度やその勾配の評価などの中間量を含む出力にアクセスすることができるようにするものである。

この vignette では、rstan パッケージに含まれる機能を簡潔に紹介している。Stanのウェブサイト mc-stan.orgには、さらなる詳細と、StanとRStanを含む多くのインターフェースの両方を操作する方法に関する最新情報が掲載されている。例えば、RStan Getting Started (The Stan Development Team 2014) を参照。

Prerequisites

Stan にはモデリング言語があり、これはベイズグラフィカルモデリングパッケージ BUGS (Lunn et al. 2000) のものと似ているが同一ではない。パーサーは Stan 言語で表現されたモデルを C++ コードに変換し、実行可能なプログラムにコンパイルして R の動的共有オブジェクト（DSO）としてロードし、ユーザーが呼び出すことができるようにする。

C++ コンパイラ、例えば g++ または clang++のような C++ コンパイラがこの処理に必要である。RStan で使用する C++ コンパイラのインストール方法については、 RStan-Getting-Startedを参照。

また、rstan パッケージは、他のいくつかの R パッケージに大きく依存している。

これらの依存関係は、従来のメカニズムで rstan パッケージをインストールした場合、自動的にインストールされるはずである。

典型的なワークフロー

以下は、ベイズ推定にRStan経由でStanを使用する場合の典型的なワークフローである。

  1. 統計モデルを、その対数事後密度（モデルの未知パラメータに依存しない正規化定数まで）を、モデリング言語 Stan を用いて記述することで表現する。拡張子 .stan の別ファイルの使用を推奨するが、R内の文字列でも可能。
  2. Stan プログラムを stanc 関数で C++ コードに変換する。
  3. C++ コードをコンパイルして、R がロードできる DSO（ダイナミックリンクライブラリ（DLL）とも呼ばれる）を作成する。
  4. DSOを実行して事後分布からサンプリングする。
  5. MCMC連鎖の非収束を診断する。
  6. 事後標本（事後分布からのMCMCドロー）に基づく推論を行う。

便利なことに、上記のステップ2、3、4はすべて、stan 関数を1回呼び出すだけで暗黙のうちに実行される。

例

この vignette の残りの部分では、:2003の5.5節で説明されている階層的メタ分析モデルを使用する。大学入学試験におけるコーチングプログラムの効果をモデル化するために、階層モデルが使用される。下の表に示すデータは、8つの高校で行われた実験結果をまとめたもので、それぞれについて標準誤差を推定している。これらのデータとモデルは、完全ベイズ推論の例として歴史的に興味深いものである (Rubin 1981) 。略して、これを_Eight Schools_の例と呼ぶ。

School Estimate (\(y_j\)) Standard Error (\(\sigma_j\))
A 28 15
B 8 10
C -3 16
D 7 11
E -1 9
F 1 11
G 18 10
H 12 18

ここで8校の例を使うのは、単純ではあるが、この研究でもともと関心のあるパラメータ（8校のそれぞれにおけるコーチングの効果）と、モデル化された集団におけるこれらの効果の変動を表すハイパーパラメータの間に依存性があるという、自明ではないマルコフ連鎖シミュレーションの問題を表しているためである。 ギブスサンプラーやハミルトニアンモンテカルロサンプラーのある種の実装は、この例では収束が遅くなることがある。

注目の統計モデルは次のように規定される。

\[ \begin{aligned} y_j &\sim \mathsf{Normal}(\theta_j, \sigma_j), \quad j=1,\ldots,8 \\ \theta_j &\sim \mathsf{Normal}(\mu, \tau), \quad j=1,\ldots,8 \\ p(\mu, \tau) &\propto 1, \end{aligned} \]

ここで、各 \(\sigma_j\) は既知とする。

Stanプログラムを書く

RStanでは、Stanプログラムをテキストファイル（通常、接尾辞は .stan ）またはR文字ベクタ（長さは1）でコード化することができる。私たちは、Eight Schoolsモデルの次のコードをファイル 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);
}

上記の Stan プログラムの最初のセクション、data ブロックは、ベイズルールの条件となるデータを指定する：学校数、 \(J\)、推定値のベクトル、 \((y_1, \ldots, y_J)\)、推定値の標準誤差のベクトル、 \((\sigma_{1}, \ldots, \sigma_{J})\)。データは整数または実数で宣言され、次元が指定されていればベクタ（より一般的には配列）であってもよい。例えば、上記のモデルでは、 \(J\) は少なくとも \(1\) でなければならず、 \(\sigma_y\) の成分はすべて正でなければならないという制約がある。

parameters ブロックは、事後分布が求められるパラメータを宣言している。これらは、学校効果の平均値 \(\mu\) と標準偏差 \(\tau\)、および標準化された学校レベルの効果 \(\eta\) である。このモデルでは、標準化されていない学校レベル効果 \(\theta\) は、パラメータとして \(\theta\) を直接宣言するのではなく、標準化効果を \(\tau\) でスケーリングし、 \(\mu\) でシフトすることで構築される変換パラメータとする。このようにモデルをパラメータ化することで、サンプラーはより効率的に実行される。なぜなら、結果として得られる多変量幾何学は、ハミルトン・モンテカルロ (Neal 2011) に適合しやすいことがある。

最後に、model のブロックは標準的な統計表記と似ている。 (ただ、注意していただきたいのは、Stanの正規 \((\cdot,\cdot)\) 分布の第2引数は標準偏差であり、統計表記で通常使われる分散ではない).私たちはモデルをベクタ表記で書いた。これにより、Stanはより効率的なアルゴリズム微分（AD）を利用できるようになった。また、normal_lpdf(y | theta,sigma) を \(J\) のループに置き換えてモデルを書くことも可能である（ただし、効率は落つ）。

for (j in 1:J) 
  target += normal_lpdf(y[j] | theta[j],sigma[j]);

Stan は、確率分布、行列演算、様々な特殊関数など、統計モデリングに最も有用な R 関数の多く のバージョンを備えている。しかし、Stan 関数の名前は R 関数の名前と異なることがあり、さらに、Stanの確率分布のパラメータ化は、同じ分布についてR関数のものと異なることがある。この問題を軽減するために、lookup 関数に R 関数または R 関数の名前を示す文字列を渡すと、RStan は対応する Stan 関数を検索してその引数を表示し、The Stan Development Team (2016) でその関数が説明されているページ番号を与えようとする。

lookup("dnorm")
          StanFunction
374 normal_id_glm_lpmf
375 normal_id_glm_lpmf
376      normal_id_glm
379        normal_lpdf
380             normal
                                                       Arguments ReturnType
374   (vector y , matrix x, real alpha, vector beta, real sigma)       real
375 (vector y , matrix x, vector alpha, vector beta, real sigma)       real
376                                                            ~       real
379                            (reals y , reals mu, reals sigma)       real
380                                                            ~       real
lookup(dwilcox)   # no corresponding Stan function
[1] "no matching Stan functions"

lookup 関数は、Stan 関数に対応する R 関数を見つけられなかった場合、引数を正規表現として扱い、Stan 関数名との一致を探そうとする。

データの準備

stan 関数は、データを名前付きリスト、オブジェクト名の文字ベクタ、または environment として受け取る。また、data 引数を省略すると、R は Stan プログラムの data ブロックで宣言されたものと同じ名前を持つオブジェクトを検索する。以下は、Eight Schools の例のデータである。

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

また、Rスクリプトに直接数値を入力するのではなく、ファイルからデータを読み込むことも可能（むしろ推奨）であろう。

事後分布からのサンプル

次に、stan 関数を呼び出して、事後標本を描く。

library(rstan)
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 = 1,              # number of cores (could use one per chain)
  refresh = 0             # no progress shown
  )
Trying to compile a simple C file
Warning: There were 1 divergent transitions after warmup. See
https://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
to find out why this is a problem and how to eliminate them.
Warning: Examine the pairs() plot to diagnose sampling problems

stan 関数は、以下の3つのステップを包んでいる。

stan を1回呼び出すだけで、3つのステップをすべて実行するが、1つずつ実行することもできる ( stanc , stan_model , sampling のヘルプページを参照)。 これはデバッグに便利である。さらに、Stan は DSO を保存するので、同じモデルを再び適合させるとき (新しいデータや設定で) 再コンパイルする必要がない。モデルがコンパイルされた後、サンプリング前にエラーが発生した場合（例えば、データや初期値などの入力に問題があった場合）でも、コンパイルされたモデルを再利用することができる。

stan 関数は stanfit オブジェクトを返するが、これはクラス "stanfit" の S4 オブジェクトである。R のクラスと S4 クラスの概念に馴染みのない方は、Chambers (2008) を参照。S4 クラスは、オブジェクトをモデル化するためのいくつかの属性（データ）と、オブジェクトの振る舞いをモデル化するためのいくつかのメソッドから構成されている。ユーザーの立場からすると、スタンフィットオブジェクトが作成された後は、どのようなメソッドが定義されているのかが主な関心事となる。

エラーが発生しない場合、返された stanfit オブジェクトは、モデルパラメータとモデルで定義された他の量に対する事後分布から描かれたサンプルを含んでいる。エラーが発生した場合（Stanプログラムのシンタックスエラーなど）、stan は終了するか、事後描画を含まないstanfitオブジェクトが返されるかのいずれかである。

クラス "stanfit" については、print や plot など、MCMC サンプルで作業するための多くのメソッドが定義されている。例えば、以下は、print メソッドを用いた Eight Schools モデルからのパラメーターの要約である。

print(fit1, pars=c("theta", "mu", "tau", "lp__"), probs=c(.1,.5,.9))
Inference for Stan model: anon_model.
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.38    0.15 8.10   2.63  10.23  21.93  2787    1
theta[2]   7.85    0.09 6.31   0.03   7.76  15.64  4615    1
theta[3]   6.15    0.12 7.49  -3.01   6.67  14.65  3632    1
theta[4]   7.78    0.10 6.44  -0.09   7.90  15.65  4493    1
theta[5]   5.15    0.09 6.22  -3.26   5.61  12.37  4525    1
theta[6]   6.28    0.10 6.65  -2.14   6.68  13.87  4231    1
theta[7]  10.61    0.12 6.65   2.84  10.07  19.29  3261    1
theta[8]   8.36    0.13 7.72  -0.43   8.33  17.62  3387    1
mu         7.86    0.10 4.85   1.85   7.90  13.82  2253    1
tau        6.44    0.14 5.44   0.92   5.15  13.53  1578    1
lp__     -39.51    0.07 2.64 -43.01 -39.26 -36.33  1374    1

Samples were drawn using NUTS(diag_e) at Sat Nov 26 04:37:52 2022.
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).

この出力の最後の行、lp__、は Stan によって計算された（正規化されていない）事後密度の対数である。 この対数密度は、モデルの評価や比較に様々な形で利用することができる（例えば、Vehtari and Ojanen (2012) を参照）。

スタン(stan)関数への引数 {#arguments-to-the-stan-function}

サンプリング（関数 stan および sampling ）の主な引数は、データ、初期値、および chains、iter、warmup などのサンプラーのオプションである。特に warmup は、NUTSサンプラーがサンプリング開始前の適応フェーズに使用する反復回数を指定しま す。 ウォームアップの後、サンプラーはアダプテーションをオフにし、合計 iter のイテレーション（ warmup を含む）が完了するまで継続する。ウォームアップ中に得られたドローが事後分布からのものであるという理論的な保証はないので、ウォームアップのドローは診断にのみ使用し、推論には使用しないでみよう。print 法で示されるパラメータの要約は、ウォームアップ後の描画のみを用いて計算されている。

オプションの引数 init は、マルコフ連鎖の初期値を指定するために使用することができる。初期値を指定する方法はいくつかあり、詳細は stan 関数のドキュメントに記載されている。ほとんどの場合、Stanがそれ自身の初期値をランダムに生成するようにするのが適切である。しかし、Stanプログラムの parameters ブロックで宣言されたオブジェクトの少なくともサブセットの初期値を指定した方が良い場合がある。

Stan は、並列処理をサポートする乱数生成器(RNG)を使用する。RNGの初期化は引数 seed と chain_id で決定される。stan 関数の1回の呼び出しから複数の連鎖を実行する場合でも、指定する必要があるのは1つの種のみで、指定しない場合はRによってランダムに生成される。

データの前処理と受け渡し

stan に渡されたデータは、前処理を経る。この前処理の詳細については、stan 関数のドキュメントに記載されている。ここでは、いくつかの重要なステップを強調する。まず、RStanでは、data ブロックで宣言されている以上のオブジェクトをデータとして渡すことができる（不要なオブジェクトは黙って省かれる）。一般に、RからStanに渡されるデータのリストの要素は数値であるべきで、その次元はモデルの data ブロックの宣言と一致しなければならない。そのため、例えばRの factor 型はRStanのデータ要素としてサポートされておらず、as.integer を介して整数コードに変換する必要がある。Stanモデリング言語では、整数と倍数を区別している（Stanモデリング言語ではそれぞれ int 型と real 型）。stan 関数は、いくつかの R データ（通常は倍精度である）を可能な限り整数に変換する。

Stan 言語には、スカラーや、スカラーの集合であるベクトル、行列、配列などの型があ る。R は真のスカラーを持たないので、RStan は長さ1のベクタをスカラーとして扱いる。しかし、data ブロックが定義されたモデルを考えてみよう。

data {                
  int<lower=1> N;      
  real y[N];
}

ここで、N は特別な場合として \(1\) とすることができる。つまり、N は常に \(1\) よりも大きいことが分かっていれば、y のデータ入力として、Rの長さ N のベクタ（例えば y <- rnorm(N) で作成したベクタ）を使うことができるのである。 \(N`\) が \(1\) のとき、RStan が y の入力データをスカラーとして扱わないようにするには、次のRコードのように明示的に配列にする必要がある。

y <- as.array(y)

Stan はデータの欠損値を自動的に処理できないので、データのどの要素も NA の値を含むことはできない。RStan のデータ前処理の重要なステップは、欠損値をチェックし、欠損値が見つかったらエラーを出すことである。しかし、欠損値を考慮した Stan プログラムを書く方法はいろいろある（The Stan Development Team (2016) 参照）。

"stanfit" クラス用メソッド {#methods-for-the-"stanfit"-class}

rstan パッケージに含まれる他の vignette では、stanfit オブジェクトについてより詳しく説明し、オブジェクトに含まれる最も重要なコンテンツ (例: posterior draws, 診断要約) にアクセスする例を示している。また、利用可能なメソッドの完全なリストは、"stanfit" クラスのドキュメントで見つけることができる。 help("stanfit "," rstan") .ここでは、いくつかの例のみを示す。

stanfit オブジェクトの plot メソッドは、出力の様々なグラフィカルな概観を提供する。デフォルトのプロットは、lp__ （加法定数までの事後密度関数の対数）と同様に、すべてのパラメータの事後不確定性区間（デフォルトでは80% (inner) と95% (outer) ）と事後中央値を表示す。

plot(fit1)
'pars' not specified. Showing first 10 parameters by default.
ci_level: 0.8 (80% intervals)
outer_level: 0.95 (95% intervals)

オプションの引数 plotfun は、利用可能な様々なプロットの中から選択するために使用される。help("plot,stanfit-method") を参照。

traceplot メソッドは、事後ドローの時系列をプロットするために使用される。 inc_warmup=TRUE を設定してウォームアップ描画を含めると、ウォームアップ領域の背景色がウォームアップ後と異なる。

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

マルコフ連鎖の収束を評価するために、トレースプロットを目視で確認することに加えて、split \(\hat{R}\) 統計量を計算することができる。split \(\hat{R}\) は Gelman and Rubin (1992) で提案された \(\hat{R}\) 統計量の更新版で、各鎖を2つの半分に分割することに基づいている。詳細は Stan のマニュアルを参照。各パラメータの推定値 \(\hat{R}\) は summary と print の出力に列の一つとして含まれている。

print(fit1, pars = c("mu", "tau"))
Inference for Stan model: anon_model.
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  2.5%  25%  50% 75% 97.5% n_eff Rhat
mu  7.86    0.10 4.85 -1.76 4.74 7.90  11 17.47  2253    1
tau 6.44    0.14 5.44  0.26 2.43 5.15   9 20.03  1578    1

Samples were drawn using NUTS(diag_e) at Sat Nov 26 04:37:52 2022.
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).

繰り返しになるが、詳細については、stanfit オブジェクトに関する追加の vignette を参照。

サンプリングの難しさ

モデルの出力を可視化する最良の方法は、 shinystan RパッケージからアクセスできるSinyStanインターフェイスを使うことである。 ShinyStan はパラメータ分布の可視化とサンプラーの問題診断の両方を容易にする。 shinystan パッケージのドキュメントには、stanfit オブジェクトでインターフェイスを使用するための指示が記載されている。

ShinyStan の使用に加えて、 rstan パッケージの関数を使用して、いくつかのサンプリング問題を診断することが可能である。get_sampler_params 関数は、サンプラーの性能に関連するパラメータの情報を返す。

# all chains combined
sampler_params <- get_sampler_params(fit1, inc_warmup = TRUE)
summary(do.call(rbind, sampler_params), digits = 2)
 accept_stat__    stepsize__     treedepth__   n_leapfrog__  divergent__    
 Min.   :0.00   Min.   :0.034   Min.   :0.0   Min.   :  1   Min.   :0.0000  
 1st Qu.:0.77   1st Qu.:0.290   1st Qu.:3.0   1st Qu.:  7   1st Qu.:0.0000  
 Median :0.95   Median :0.352   Median :3.0   Median : 15   Median :0.0000  
 Mean   :0.83   Mean   :0.392   Mean   :3.4   Mean   : 12   Mean   :0.0085  
 3rd Qu.:0.99   3rd Qu.:0.423   3rd Qu.:4.0   3rd Qu.: 15   3rd Qu.:0.0000  
 Max.   :1.00   Max.   :9.269   Max.   :7.0   Max.   :159   Max.   :1.0000  
    energy__ 
 Min.   :35  
 1st Qu.:42  
 Median :44  
 Mean   :45  
 3rd Qu.:47  
 Max.   :63  
# each chain separately
lapply(sampler_params, summary, digits = 2)
[[1]]
 accept_stat__    stepsize__     treedepth__   n_leapfrog__  divergent__   
 Min.   :0.00   Min.   :0.048   Min.   :0.0   Min.   : 1    Min.   :0.000  
 1st Qu.:0.78   1st Qu.:0.329   1st Qu.:3.0   1st Qu.: 7    1st Qu.:0.000  
 Median :0.95   Median :0.329   Median :4.0   Median :15    Median :0.000  
 Mean   :0.83   Mean   :0.360   Mean   :3.5   Mean   :13    Mean   :0.004  
 3rd Qu.:0.99   3rd Qu.:0.352   3rd Qu.:4.0   3rd Qu.:15    3rd Qu.:0.000  
 Max.   :1.00   Max.   :3.545   Max.   :6.0   Max.   :63    Max.   :1.000  
    energy__ 
 Min.   :36  
 1st Qu.:42  
 Median :44  
 Mean   :44  
 3rd Qu.:47  
 Max.   :59  

[[2]]
 accept_stat__    stepsize__     treedepth__   n_leapfrog__  divergent__   
 Min.   :0.00   Min.   :0.046   Min.   :0.0   Min.   : 1    Min.   :0.000  
 1st Qu.:0.70   1st Qu.:0.409   1st Qu.:3.0   1st Qu.: 7    1st Qu.:0.000  
 Median :0.92   Median :0.423   Median :3.0   Median : 7    Median :0.000  
 Mean   :0.79   Mean   :0.445   Mean   :3.1   Mean   :10    Mean   :0.011  
 3rd Qu.:0.98   3rd Qu.:0.423   3rd Qu.:3.0   3rd Qu.:15    3rd Qu.:0.000  
 Max.   :1.00   Max.   :9.269   Max.   :6.0   Max.   :63    Max.   :1.000  
    energy__ 
 Min.   :36  
 1st Qu.:42  
 Median :44  
 Mean   :44  
 3rd Qu.:46  
 Max.   :59  

[[3]]
 accept_stat__    stepsize__     treedepth__   n_leapfrog__  divergent__   
 Min.   :0.00   Min.   :0.051   Min.   :0.0   Min.   :  1   Min.   :0.000  
 1st Qu.:0.79   1st Qu.:0.352   1st Qu.:3.0   1st Qu.:  7   1st Qu.:0.000  
 Median :0.95   Median :0.352   Median :3.0   Median : 15   Median :0.000  
 Mean   :0.84   Mean   :0.388   Mean   :3.4   Mean   : 12   Mean   :0.009  
 3rd Qu.:0.99   3rd Qu.:0.376   3rd Qu.:4.0   3rd Qu.: 15   3rd Qu.:0.000  
 Max.   :1.00   Max.   :4.982   Max.   :7.0   Max.   :159   Max.   :1.000  
    energy__ 
 Min.   :36  
 1st Qu.:42  
 Median :44  
 Mean   :45  
 3rd Qu.:47  
 Max.   :57  

[[4]]
 accept_stat__    stepsize__     treedepth__   n_leapfrog__  divergent__   
 Min.   :0.00   Min.   :0.034   Min.   :0.0   Min.   :  1   Min.   :0.000  
 1st Qu.:0.84   1st Qu.:0.290   1st Qu.:3.0   1st Qu.:  7   1st Qu.:0.000  
 Median :0.96   Median :0.290   Median :4.0   Median : 15   Median :0.000  
 Mean   :0.85   Mean   :0.375   Mean   :3.5   Mean   : 13   Mean   :0.011  
 3rd Qu.:0.99   3rd Qu.:0.412   3rd Qu.:4.0   3rd Qu.: 15   3rd Qu.:0.000  
 Max.   :1.00   Max.   :5.780   Max.   :6.0   Max.   :127   Max.   :1.000  
    energy__ 
 Min.   :35  
 1st Qu.:42  
 Median :44  
 Mean   :45  
 3rd Qu.:47  
 Max.   :63

ここでは、少数の発散遷移があることがわかる。これは、divergent__ が \(1\) であることによって識別される。理想的には、ウォームアップ段階以降に発散する遷移がないことが望ましい。 \(0.8\) この場合、accept_stat__ の平均はすべてのチェーンで \(0.8\) に近いが、中央値が \(0.95\) の近くにあるため、非常に歪んだ分布になっている。もう一度 stan を呼び出し、オプションの引数 control=list(adapt_delta=0.9) を指定して、発散する推移性を排除しようとすることもできる。 しかし、目標合格率が高いとき、ステップサイズが非常に小さく、サンプラーが反復あたり取れるリープフロッグステップ数の限界にぶつかることがある。この場合、各鎖の treedepth__ は最大でも \(7\) であり、デフォルトは \(10\) であるため、問題にはならない。しかし、もし treedepth__ の中に \(11\) があれば、control=list(max_treedepth=12) (例えば) を stan に渡して、制限を増やすのが賢明だろう。get_sampler_params が返すオブジェクトの構造については、stanfitオブジェクトに関する vignette を参照。

また、pairs を使って、同じ情報をグラフにすることもできる。pairs プロットは、サンプリングの難しさがテールや最頻値付近で起こっているかどうかを知るために使うことができる。

pairs(fit1, pars = c("mu", "tau", "lp__"), las = 1)
Warning in par(usr): argument 1 does not name a graphical parameter

Warning in par(usr): argument 1 does not name a graphical parameter

Warning in par(usr): argument 1 does not name a graphical parameter

上のプロットでは、選択された各パラメータの周辺分布が対角線に沿ったヒストグラムとして含まれている。デフォルトでは、中央値以下の accept_stat__ (MCMC提案受理率)を持つドローは対角線の下に、中央値以上の accept_stat__ を持つものは対角線の上にプロットされる（これは condition 引数を用いて変更可能）。各対角線外の正方形は、行変数と列変数の交点に対するドローの2変量分布を表している。理想的には、同じ2つの変数の対角線下の交差点と対角線上の交差点は、互いに鏡像の分布を持つべきである。黄色の点は最大値 treedepth__ に当たった推移性を示し、赤色の点は発散した推移性を示すことになる。

その他のトピック

ユーザー定義スタン関数

Stan では、Stan プログラム全体で使用可能な独自の関数を定義することもできる。これらの関数は、functions ブロックで定義される。 functions ブロックはオプションであるが、存在する場合は他のどのブロックよりも前になければならない。このメカニズムにより、ユーザーは統計分布や現在Stanで利用できない他の機能を実装することができる。しかし、ユーザーの関数が既存のStan関数への呼び出しをラップしているだけであっても、1つ（あるいは2つ）のタスクを達成する数行の Stan コードがユーザー定義関数への呼び出しに置き換えられれば、model ブロックのコードははるかに読みやすくなる。

ユーザー定義関数を利用するもう一つの理由は、RStan がそのような関数を R グローバル環境にエクスポートするための expose_stan_functions 関数を提供し、それらが適切に動作していることを確認するために R でテストすることができることである。例えば

model_code <-
'
functions {
  real standard_normal_rng() {
    return normal_rng(0,1);
  }
}
model {}
'
expose_stan_functions(stanc(model_code = model_code))
Trying to compile a simple C file
Running \
  /Volumes/HDD/opt/sw/Library/Frameworks/R.framework/Versions/4.2/Resources/bin/R \
  CMD SHLIB foo.c
clang -I"/Volumes/HDD/opt/sw/Library/Frameworks/R.framework/Versions/4.2/Resources/include" -DNDEBUG   -I"/Users/baba/Library/R/x86_64/4.2/library/Rcpp/include/"  -I"/Users/baba/Library/R/x86_64/4.2/library/RcppEigen/include/"  -I"/Users/baba/Library/R/x86_64/4.2/library/RcppEigen/include/unsupported"  -I"/Users/baba/Library/R/x86_64/4.2/library/BH/include" -I"/Users/baba/Library/R/x86_64/4.2/library/StanHeaders/include/src/"  -I"/Users/baba/Library/R/x86_64/4.2/library/StanHeaders/include/"  -I"/Users/baba/Library/R/x86_64/4.2/library/RcppParallel/include/"  -I"/Users/baba/Library/R/x86_64/4.2/library/rstan/include" -DEIGEN_NO_DEBUG  -DBOOST_DISABLE_ASSERTS  -DBOOST_PENDING_INTEGER_LOG2_HPP  -DSTAN_THREADS  -DUSE_STANC3 -DSTRICT_R_HEADERS  -DBOOST_PHOENIX_NO_VARIADIC_EXPRESSION  -DBOOST_NO_AUTO_PTR  -include '/Users/baba/Library/R/x86_64/4.2/library/StanHeaders/include/stan/math/prim/fun/Eigen.hpp'  -D_REENTRANT -DRCPP_PARALLEL_USE_TBB=1   -I/Volumes/HDD/opt/sw/include   -fPIC  -g -O3  -c foo.c -o foo.o
In file included from <built-in>:1:
In file included from /Users/baba/Library/R/x86_64/4.2/library/StanHeaders/include/stan/math/prim/fun/Eigen.hpp:22:
In file included from /Users/baba/Library/R/x86_64/4.2/library/RcppEigen/include/Eigen/Dense:1:
In file included from /Users/baba/Library/R/x86_64/4.2/library/RcppEigen/include/Eigen/Core:88:
/Users/baba/Library/R/x86_64/4.2/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:628:1: error: unknown type name 'namespace'
namespace Eigen {
^
/Users/baba/Library/R/x86_64/4.2/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:628:16: error: expected ';' after top level declarator
namespace Eigen {
               ^
               ;
In file included from <built-in>:1:
In file included from /Users/baba/Library/R/x86_64/4.2/library/StanHeaders/include/stan/math/prim/fun/Eigen.hpp:22:
In file included from /Users/baba/Library/R/x86_64/4.2/library/RcppEigen/include/Eigen/Dense:1:
/Users/baba/Library/R/x86_64/4.2/library/RcppEigen/include/Eigen/Core:96:10: fatal error: 'complex' file not found
#include <complex>
         ^~~~~~~~~
3 errors generated.
make: *** [foo.o] Error 1
standard_normal_rng()
[1] -0.9529876

Log-Posterior (関数と勾配) {#the-log-posterior-(function-and-gradient)}

Stan は、未知の加法定数までの事後分布の確率密度関数の対数を定義する。各反復におけるこの対数カーネルの実測値を表すために、lp__ を使用する（そして、lp__ は要約と分割 \(\hat{R}\) および有効サンプルサイズの計算において、未知数として扱われる）。

rstan パッケージの素晴らしい特徴は、与えられた stanfit オブジェクトに対して lp__ とその勾配の両方を計算するための関数を公開していることである。この2つの関数はそれぞれ log_prob と grad_log_prob である。どちらも、パラメータのサポートが実線全体でない場合でも、_unconstrained_空間上のパラメータを取る。Stanのマニュアル (The Stan Development Team 2016) には、Stanが実線全体からその部分空間へ（またはその逆）マッピングするために使用する特定の変換についての詳細がある。

制約のないパラメータの数は、パラメータの総数よりも少なくなる可能性がある。例えば、長さが \(K\) の単項の場合、単項のすべての要素は非負で和が 1 でなければならないという制約があるため、実際には制約のないパラメータは \(K-1\) だけである。 get_num_upars メソッドは制約のないパラメータの数を得るために提供されており、unconstrain_pars と constrain_pars メソッドはそれぞれパラメータの制約のない値と制約のある値を計算するために使用することができる。 前者はパラメータのリストを入力として受け取り、それを制約のないベクタに変換し、後者はその逆を行う。これらの関数を用いることで、ベイズモデルの最大事後推定など、他のアルゴリズムを実装することができる。

スタンでの最適化

RStan は Stan の最適化機能へのインターフェースも提供しており、Stan プログラムによって定義された（おそらくペナルティ付きの）尤度関数を最大化することによって、点推定値を得るために使用することができる。このフィーチャは、標準偏差が既知の正規分布から得られたと仮定した標本から平均を推定するという、非常にシンプルな例で説明される。つまり

\[y_n \sim \mathsf{Normal}(\mu,1), \quad n = 1, \ldots, N. \]

事前分布 \(p(\mu) \propto 1\) を指定すると、 \(\mu\) の最大事後推定量は、単に標本平均となる。 \(p(\mu) \propto 1\) は事前分布を指定しない場合のデフォルトなので、 \(\mu\) に対してこの事前分布を明示的にコード化する必要はない。

まず、"stanmodel" クラスのオブジェクトを作成し、optimizing メソッドを使用する。このメソッドには、データなどの引数を与えることができる。

ocode <- "
  data {
    int<lower=1> N;
    real y[N];
  } 
  parameters {
    real mu;
  } 
  model {
    target += normal_lpdf(y | mu, 1);
  } 
"

sm <- stan_model(model_code = ocode)
Trying to compile a simple C file
y2 <- rnorm(20)
mean(y2)
[1] 0.5575213
optimizing(sm, data = list(y = y2, N = length(y2)), hessian = TRUE)
$par
       mu 
0.5575213 

$value
[1] -24.30839

$return_code
[1] 0

$hessian
    mu
mu -20

$theta_tilde
            mu
[1,] 0.5575213

モデルコンパイル

vignette で述べたように、Stan プログラムは Stan モデリング言語で書かれ、C++コードに変換された後、動的共有オブジェクト（DSO）にコンパイルされる。DSOはRに読み込まれ、実行され、事後標本が描かれる。C++のコードをDSOにコンパイルする工程は、時に時間がかかる。モデルが同じであれば、前回実行したDSOを再利用することができる。stan 関数はオプション引数として fit を受け付ける。この引数は、コンパイルしたモデルを再利用するように、既存のフィットしたモデルオブジェクトを渡すために使用する。以前の適合モデルを再利用する場合でも、stan の他の引数に異なる値を指定することができ、data の引数に異なるデータを渡すことができる。

さらに、フィットしたモデルが save や save.image のような関数を使用して保存されている場合、RStan は DSO を保存することができ、R セッションをまたいで使用することができる。DSOを保存しないようにするには、stan 関数を呼び出すときに save_dso=FALSE を指定する。

ユーザが rstan_options(auto_write = TRUE) を実行すると、コンパイルされたモデルのシリアライズ版が自動的にハードディスクの .stan ファイルと同じディレクトリに、Stan プログラムが文字列で表現されている場合は R のテンポラリディレクトリに保存される。このオプションはCRANポリシーによりデフォルトでは有効ではないが、冗長なコンパイルを排除するために、通常はユーザが指定する必要がある。

Stan は、コードが最大レベルの最適化でコンパイルされると、はるかに高速に実行される。これは、ほとんどの C++ コンパイラで -O3 である。しかし、R のデフォルト値は -O2。これはほとんどの R パッケージで適切であるが、Stan では若干の速度低下を伴う。このデフォルトは、 CRAN - Customizing-package-compilationの指示に従ってローカルに変更することができる。 ただし、CXXFLAGS = -O3 を設定すると、他のRパッケージに悪影響が出る可能性があるので、注意が必要である。

Stan プログラムのパースとコンパイルの詳細については、stanc と stan_model 関数のドキュメントを参照。

複数チェインの並列実行

実行するマルコフ連鎖の数は、stan または sampling 関数の chains 引数を使用して指定できる。デフォルトでは、連鎖は親Rプロセスを用いてシリアルに（すなわち、一度に1つずつ）実行される。また、オプションの引数 cores、鎖の数を設定することができる（ハードウェアが十分なプロセッサと RAM を持っている場合）、これはほとんどのラップトップで適切である。cores 引数を手動で指定する必要がなく、利用可能なすべてのコアを使用できるように、通常、最初に R セッションごとに一度 options(mc.cores=parallel::detectCores()) を呼び出すことを勧める。

異なる並列化スキーム(おそらくリモートクラスタで)で作業するユーザーのために、 rstan パッケージは、複数のstanfitオブジェクトのリスト(同じ Stan プログラムから作成し、同じ数のウォームアップとサンプリングの反復を使用)を単一の stanfit オブジェクトに統合するための sflist2stanfit という関数を提供している。すべての連鎖に同じ種を指定することが重要であり、異なる連鎖 ID（引数 chain_id ）を使用することも同様に重要である。これらの組み合わせにより、すべての連鎖に対して Stan で生成される乱数が本質的に独立であることが保証される。これは次のような場合に自動的に(内部的に)処理される。 $ cores > 1$ .

CmdStanを使った作業

rstan__ パッケージは、Stan のコマンドラインインターフェイスである CmdStan のデータを作成し、出力を読み取るためのいくつかの関数を提供する。

まず、Stan がデータや初期値を読み込むとき、R のダンプデータフォーマットの構文のサブセットをサポートしている。そのため、ベースRの dump 関数を使ってデータを用意すると、Stan が中身を読めなくなることがある。 rstan にある stan_rdump 関数は、dump 関数と非常によく似たセマンティクスで、R から Stan がサポートするフォーマットへデータをダンプするように設計されている。

次に、read_stan_csv 関数は、CmdStan が生成した CSV ファイルを読み込んで stanfit オブジェクトを作成する。作成された stanfit オブジェクトは、診断や事後解析の様々な手法に対応している。

See Also

Chambers, John M. 2008. Software for Data Analysis : Programming with r. New York: Springer.
Gelman, Andrew, and Donald B. Rubin. 1992. “Inference from Iterative Simulation Using Multiple Sequences.” Statistical Science 7 (4): 457–72.
Hoffman, Matthew D., and Andrew Gelman. 2012. “The No-U-Turn Sampler: Adaptively Setting Path Lengths in Hamiltonian Monte Carlo.” Journal of Machine Learning Research.
Lunn, D. J., A. Thomas, N. Best, and D. Spiegelhalter. 2000. “WinBUGS — a Bayesian Modelling Framework: Concepts, Structure, and Extensibility.” Statistics and Computing, 325–37.
Neal, Radford. 2011. “MCMC Using Hamiltonian Dynamics.” In Handbook of Markov Chain Monte Carlo, edited by Steve Brooks, Andrew Gelman, Galin L. Jones, and Xiao-Li Meng, 116–62. Chapman; Hall/CRC.
R Core Team. 2014. R: A Language and Environment for Statistical Computing. Vienna, Austria: R Foundation for Statistical Computing. https://www.R-project.org/.
Rubin, Donald B. 1981. “Estimation in Parallel Randomized Experiments.” Journal of Educational and Behavioral Statistics 6 (4): 377–401.
The Stan Development Team. 2014. “RStan Getting Started.” https://mc-stan.org/.
———. 2016. Stan Modeling Language: User’s Guide and Reference Manual. https://mc-stan.org.
Vehtari, A., and J. Ojanen. 2012. “A Survey of Bayesian Predictive Methods for Model Assessment, Selection and Comparison.” Statistics Surveys 6: 142–228.