特異モデルに対してBridge Samplingで周辺尤度の推定は可能か

暫定的な見解：難しいのでは？（[編集後記 2019/9/28] warp-U bridge samplingを使うとうまくいくかもしれない。推定の分散が他の方法と比べてどれだけ小さいかが気になるところ）

定義：周辺尤度

まず言葉の定義をすると、モデルを $p(X|\theta)$ , 事前分布を $p(\theta)$ としたとき、周辺尤度 $p(X)$ を
$p(X):= \int p(X|\theta) p(\theta) d\theta$
と定義する。

周辺尤度の評価方法

周辺尤度は積分評価なので以下の評価が考えられる：

積分を陽に評価
数値積分(sampling technique)で評価
積分を漸近的に評価

1は古典的なアプローチでは難しいのはみなさんご存じの通り。現代的なアプローチとして代数幾何を使うものもある[1]が一般的に適用するには難しい。
2は次元が高い時には次元の呪いで誤差が累積したり分散が大きくなってしまったりします。
3はBIC, WBIC, sBICなど情報量規準的なアプローチ。

この記事の目的

この記事は上記の2の評価の1つであるbridge samplingが高い評価を得ているようなので簡単な検証をしてみました。
bridge samplingについては岡田先生のこちらのreviewを参考にするとよいかと思います。具体的な計算方法は清水先生のこちらの記事を見てください。どうやら心理学系の方でよく使われているようです。この記事を書こうと思ったモチベーションの1つは清水先生のこちらのスライドの中でWBICとbridge samplingではおそらくbridge samplingの方が精度が良いであろうというところに関する検証になります。

正則モデルのbridge sampling

まずexactな周辺尤度が知られている以下のケースでbridge samplingが正しく推定できるかを見てみます。radiata pine datasetというデータを使って以下の回帰モデルを考えます:
$y_i= \alpha + \beta (x_i - \bar{x}) + \epsilon_i, \quad \epsilon_i \sim N(0, \tau^{-1})$
事前分布の設定についてはこちらの論文に従いました。exactな対数周辺尤度は論文内にも記載があり、また著者のHPにもその計算プログラムが公開されています。
bridge samplingは以下で計算させました。

#radiata.stan
data {
  int<lower=1> N;
  vector[N] Y;
  vector[N] W;
  real meanW;
}

parameters {
  real<lower=0> tau;
  real alpha;
  real beta;
}

model {
  target += gamma_lpdf(tau | 6, 4*300*300);
  target += normal_lpdf(alpha | 3000, 1.0/(0.06*sqrt(tau)));
  target += normal_lpdf(beta | 185, 1.0/(6*sqrt(tau)));
  for(i in 1:N)
    target += normal_lpdf(Y[i] | alpha+beta*(W[i]-meanW), 1.0/sqrt(tau));
}

library(rstan)
library(bridgesampling)

lm.model <- stan_model(file='radiata.stan')
d <- read.table("RadiataPine.csv", sep=",", header=T)
data_bs <- list(N=nrow(d), Y=d$y, W=d$x, meanW=mean(d$x))

fit <- sampling(lm.model, data=data_bs, iter=10000, chains=4)
bs <- bridge_sampler(fit, silent = T)

結果は以下の通りでした：

手法	評価値
exact	-310.128
bridge sampling	-310.679

少しずれていますがこんなもんでしょうか。

特異モデルのbridge sampling

次に特異モデルでbridge samplingをしてみましょう。ここでは混合数2の混合正規分布を考えます:
$\alpha N(\mu_1, 1) + (1-\alpha) N(\mu_2, 1)$
データ生成分布は $N(0, 1)$ とします。
このケースだとexactな周辺尤度は知られていない（はずな）ので、Nealの方法＆プログラム[2]で周辺尤度を評価してbridge samplingの評価値との比較を行います。

# normalmix.stan
data {
  int<lower=1> N;
  vector[N] Y;
}

parameters {
  real<lower=0, upper=1> a;
  ordered[2] mu;
}

model {
  mu[1] ~ normal(0,sqrt(10));
  mu[2] ~ normal(0,sqrt(10));
  for(n in 1:N)
    target += log_sum_exp(
      log(1-a) + normal_lpdf(Y[n] | mu[1], 1),
      log(a) + normal_lpdf(Y[n] | mu[2], 1)
    );
}

library(rstan)
library(bridgesampling)

N <- 50
Y <- rnorm(N, 0, 1)
gmm.model <- stan_model(file='normalmix.stan')

data_1 <- list(N=N, Y=Y)
fit_1 <- sampling(gmm.model, data=data_1, iter=10000, chains=3)
BR_tmp <- bridge_sampler(fit_1, silent = T)
bridge_res <- BR_tmp[[1]]

上記のコードをベースに100回のiteratoinでNealのMonte Carlo evaluation, bridge sampling, WBIC, WBICを補正した値, sBICの平均で比較してみたのが以下の結果です:

手法	mean
Monte Carlo evaluation	-74.11
bridge sampling	-70.70
WBIC	-73.18
WBIC - $\hat{\nu}(t_w)$	-73.72
sBIC	-75.61

このケースであればNealのtechnical reportと同じでパラメータは3次元であるので最初に述べた次元の呪いの影響はあまりないと考えられるためMonte Carlo evaluationが最も精度が高いと考えてよいと思います。それと比べるとbridge samplingは結構overestimationしているように見えます。WBICもoverestimationすると指摘されていますが、それより大きくなっています。

特異モデルのbridge samplingはなぜうまくいかないか

（コードが正しいとするならば）実装で提案分布がnormal distributionになっていることが原因ではないか（要検証）。bridge samplingの提案分布は"事後分布と類似していて重なりが大きいように選ばれる分布"[3]ということだが、特異モデルの事後分布は事前にどうなるか予想が難しいと思うので（というのもデータ生成分布に依存してしまうので）、良い感じの提案分布を一般に用意するのは難しいと思う。となると、特異モデルでbridge samplingをうまく使うのは難しいのではないか。

bridge samplingなにもわからんマンなので、間違っていたらコメントをしていただけるとありがたいです。

[編集後記 2019/9/28] [1609.07690] Warp Bridge Sampling: The Next Generation なら特異モデルでもうまくいくかもしれない。まだpreprintだからか実装コードが公開されていないので公開されたら試してみます。

参考文献：
[1] LIN, S. STURMFELS, B. AND XU, Z. (2009). Marginal likelihood integrals for mixtures of independence models, JMLR (10), 1611–1631.
[2] NEAL, R. (1999). Erroneous results in ‘Marginal likelihood from the Gibbs output’. Tech. rep., Unpublished. University of Toronto, Toronto.
[3] 岡田謙介 (2018). ベイズファクターによる心理学的仮説・モデルの評価. Japanese Psychological Review, Vol. 61, No. 1, 101–115.
[4] 階層ベイズと自由エネルギー. http://norimune.net/3107

あざらしの日常

あざらしが腹を叩く片手間で書いたメモ