Stanで1変数の積分を実行する - StatModeling Memorandum

Stan2.19でモデルの内部で1変数の積分を実行する機能が追加されましたので簡単に紹介します。Stan2.19がCRANにないとお嘆きのアナタ！GitHubからインストールすればよろし。Windowsならば、Benさんがここに書いてくれた方法でバイナリを直接インストールすることもできます。ただし、他のパッケージとの依存関係は壊れてもしーらん。

積分する関数名は現在のところintegrate_1dです。使い方はStanのマニュアルに載っているかなと期待したけどまだ載ってませんでした。reference manualには少しだけ触れられているけど、これだけではさっぱり分からないです。こういう場合は、StanのdiscouseかGitHubのリポジトリで検索して、手あたりしだい見ていくと知見が蓄積されます。結局、GitHubのリポジトリ内で新しい関数のテストをしているディレクトリがあり、そこに使用例がありました。

この記事では階層モデルの場合に、積分を使って対数尤度を算出する例を扱います。階層モデルの場合には複数の予測したい量が考えられて、それに応じた対数尤度（および情報量規準）が考えられるのでした。

statmodeling.hatenablog.com

ここでは上記の記事でいうところの「モデル2の(3)」に対応する予測を扱います。すなわち、別の新しいクラスができて新しく1人が加わる場合です。この場合、予測分布はグループの平均を積分消去した分布になります。

先にデータ生成およびStanを実行するRコードから説明します。

library(rstan)

G      <- 10
N_by_G <- 2
N      <- N_by_G * G
GID    <- rep(1:G, N_by_G)

set.seed(2)
Mu <- rnorm(G, mean=0, sd=10.0)
Y  <- rnorm(N, mean=Mu[GID], sd=2.0)

data <- list(N=N, G=G, N_by_G=N_by_G, GID=GID, Y=Y)
stanmodel <- stan_model(file='model/model.stan')
fit <- sampling(stanmodel, data=data, iter=3500, warmup=1000, seed=123)

library(loo)
ms <- rstan::extract(fit)
waic  <- waic(ms$log_lik)$waic/(2*N)
looic <- loo(ms$log_lik)$looic/(2*N)

3行目：グループの数です。
4行目： 1グループに何人いるかです。
5行目：全部で何人いるかです。ここでは20人です。
9行目：各グループの平均を生成しています。
10行目： 9行目で生成した値を平均として、各グループに含まれる人たちの値を生成しています。
16～19行目：対数尤度のサンプリング結果から(3)の場合のWAICとLOOICを算出しています。

Stanコードは以下です。

functions {
  real f(real x, real xc, real[] par, real[] y_r, int[] y_i) {
    return exp(normal_lpdf(x   | par[1], par[2]) +
               normal_lpdf(y_r | x,      par[3]));
  }

  real log_lik_marginal(real[] Y, real[] par) {
    return log(integrate_1d(f, par[1]-6*par[2], par[1]+6*par[2], par,
                            Y, {0}, 1e-3));
  }
}

data {
  int G;
  int N;
  vector[N] Y;
  int GID[N];
}

parameters {
  real mu0;
  real<lower=0> s0;
  vector[G] mu;
  real<lower=0> sigma;
}

model {
  mu0 ~ student_t(6, 0, 10);
  s0 ~ student_t(6, 0, 10);
  mu ~ normal(mu0, s0);
  Y ~ normal(mu[GID], sigma);
}

generated quantities {
  vector[N] log_lik;
  real mu_pred = normal_rng(mu0, s0);
  real y_pred = normal_rng(mu_pred, sigma);
  for (n in 1:N)
    log_lik[n] = log_lik_marginal({Y[n]}, {mu0, s0, sigma});
}

2～5行目：被積分関数です。現状ではこの引数タイプ（real, real, real[], real[], int[]）に限るようです。1番目の引数は積分で動く変数、2番目は後述、3番目の引数はパラメータの配列、4番目の引数は実数値をとるデータの配列、5番目の引数は整数値をとるデータの配列です。2番目の引数はどんなときに使うべきか不明だけど、たぶん3番目の引数とは切り離して考えたいパラメータを入れるのに使うのかなぁ…。
7～10行目：積分して対数尤度を求める関数の定義です。integrate_1d関数の、1番目の引数は関数名、2番目の引数は積分範囲の下限、3番目の引数は積分範囲の上限、4番目の引数は関数に渡すパラメータの配列、5番目の引数は関数に渡す実数値をとるデータの配列、6番目の引数は関数に渡す整数値をとるデータの配列、7番目の引数は積分の相対的な誤差をどれくらいに抑えるか、です。最近のStanでは{}で囲むと配列となり、[]で囲むとvectorとなります。
39行目：ここで対数尤度を求めています。

実際に実行してみると、Rでデータを作成するときの乱数をset.seed(1)の場合には情報量規準の算出までいきますが、set.seed(2)の場合には以下のエラーが出てサンプリングが止まってしまいます。

...(中略)...
Chain 2: 
[1] "Error in sampler$call_sampler(args_list[[i]]) : "                                                                                                                                                      
[2] "  Exception: Exception: integrate: error estimate of integral above zero 1.80655e-006 exceeds the given relative tolerance times norm of integral above zero  (in 'model2f8834cb4d93_model' at line 8)"
[3] "  (in 'model2f8834cb4d93_model' at line 39)"                                                                                                                                                           
[1] "error occurred during calling the sampler; sampling not done"

Stanのintegrate_1dは内部的にはBoostのintegrate_1dを使っているようですが、どうもやや不安定な気がします。GSLのCQUADの方が安定しています（{RcppGSL}で使う例はこちら）。また、尤度を算出するような場合、どこかに対数をいれて積分しないと値が小さくなりすぎて数値的に不安定になります。以前のlog_sum_expとシンプソンの公式を使う方法は一例ですが、特化した効率の良い方法がありそうです。今後に期待します。いや、お前がやるんだよ。