この記事はStan Advent Calendar 2016およびR Advent Calendar 2016の12月7日の記事です。StanコードとRコードは記事の最後にあります。

背景は以下です。

[1] Aki Vehtari, Andrew Gelman, Jonah Gabry (2015). Practical Bayesian model evaluation using leave-one-out cross-validation and WAIC. arXiv:1507.04544. (url)
[2] 渡辺澄夫. 広く使える情報量規準(WAIC)の続き（注４）【WAICとクロスバリデーションの違いについて】 (url)
[3] Sumio Watanabe. Comparison of PSIS Cross Validation with WAIC. (url)

leave-one-outクロスバリデーション（以下LOOCV）およびWAICは予測のよさをベースにしたモデル選択に用いられる情報量規準であり、ともに汎化誤差（Generalization Error、以下GE）の近似です。それにもかかわらず[1]では、本来性能評価では必須と思われる汎化誤差との比較がありません。実データ（真のモデルが未知の状況）で用いているためと思いますが、これではいけないように思います。この記事では僕が日常的に使用するような5つの基本的なモデルを使い、真のモデルが既知の状況でGE・LOOCV・WAICの性能比較を行いました。

具体的には[2]に以下のようなコメントがあります。なお[1]ではPareto Smoothed Importance SamplingでLOOCVを算出しており、PSISCVとも呼ばれるようです。

(0) まず同じデータに対してマルコフ連鎖モンテカルロ法を何度も行ったときの値の揺らぎを調べてみましょう．WAICの分散はISCVおよびPSISCVの分散よりも小さくなります。つまり、WAICはISCVおよびPSISCVよりもマルコフ連鎖揺らぎに対して強いということができます．

(2) しかしながら，CVもWAICも汎化誤差を推定することが本来の目的です．CVとWAICの両方の厳密値が計算できたとして,（つまりMCMC法で無限にサンプルが取れたとき），どちらの方が汎化誤差の推定として優れているのでしょうか。(中略) 我々の実験では，GEを汎化誤差とするとき，ほぼ，いつでも E[|PSISCV-GE|] > E[|WAIC-GE|] が成り立つのですが・・・。このページをご覧の皆様にはぜひ実験してみていただければと思います。なお、E[ ]は学習用データのでかたについての平均を表しています。

そこでこれら2点について検証してみました。先に「(2)のコメント」について説明します。

GE・LOOCV・WAICの比較

使用した5つのモデルとシミュレーションの手順を説明します。

重回帰

真のモデルは以下です。

　 $Y \sim Normal(1.3 - 3.1 X_1 + 0.7 X_2, 2.5)$

あてはめたモデルは以下です。

　 $Y \sim Normal(b_1 + b_2 X_1 + b_3 X_2, \sigma)$

データ点の数Nについては10,30,100,300を試しました。例としてN = 10の場合を説明します。まず乱数でデータX（すなわち $X_1, X_2$ ）を生成します。次にそのXの値を使って1000通りのYを生成します（学習用データのでかたの平均をとるため）。各YについてStanでiter=11000, warmup=1000, chains=4で実行して合計40000個のMCMCサンプルを得てGE・LOOCV・WAICと「LOOCV - GE」と「WAIC - GE」を求めました。その後、Nごとに「LOOCV - GE」と「WAIC - GE」のboxplotを描きました。

ロジスティック回帰

手順は重回帰の場合と同じです。使用したモデルだけが異なります。真のモデルは以下です。

　 $Y \sim Bernoulli(inv\_logit(0.8 - 1.1 X_1 + 0.1 X_2))$

あてはめたモデルは以下です。

　 $Y \sim Bernoulli(inv\_logit(b_1 + b_2 X_1 + b_3 X_2))$

　 $b_1,b_2,b_3 \sim Student\_t(4,0,1)$

非線形回帰ミカエリス・メンテン型

手順は重回帰の場合と同じです。使用したモデルだけが異なります。真のモデルは以下です。

　 $Y \sim Normal(10.0 X / (2.0 + X), 0.8)$

あてはめたモデルは以下です。

　 $Y \sim Normal(m X / (k + X), \sigma)$

　 $k \sim Uniform(0, 12)$

　 $m \sim Uniform(0, 20)$

真のモデルが含まれない場合

あてはめたモデルが以下の場合も試しました。

　 $Y \sim Normal(a + b X, \sigma)$

ノイズがt分布に従う重回帰

手順は重回帰の場合と同じです。使用したモデルだけが異なります。真のモデルは以下です。

　 $Y \sim Student\_t(4, 1.3 - 3.1 X_1 + 0.7 X_2, 2.5)$

あてはめたモデルは以下です。

　 $Y \sim Student\_t(4, b_1 + b_2 X_1 + b_3 X_2, \sigma)$

階層モデル

手順は重回帰の場合とおおよそ同じですが、データ点の数とモデルが異なります。グループの数を10に固定し、データ点の数Nについては20,50,130,400を試しました（それぞれ各グループで2,5,13,40人）。また真のモデルは以下です。

　 $\mu\left[g\right] \sim Normal(0, 10)$ 　　 $g = 1,2,\dots,10$

　 $Y\left[n\right] \sim Normal(\mu\left[g\left[n\right]\right], 2.0)$ 　　 $n = 1,2,\dots,N$

あてはめたモデルは以下です。

　 $\mu\left[g\right] \sim Normal(\mu_0, \sigma_0)$ 　　 $g = 1,2,\dots,10$

　 $Y\left[n\right] \sim Normal(\mu\left[g\left[n\right]\right], \sigma)$ 　　 $n = 1,2,\dots,N$

　 $\mu_0,\sigma_0 \sim Student\_t(4, 0, 10)$

階層モデルにおいては何を予測したいのか（どういう状況の汎化誤差を知りたいのか）を注意深く考える必要があります。以下の記事を参照。

statmodeling.hatenablog.com

ここでは以下の2つの場合を計算しました。

既存の各グループに、新しく1人ずつ加わる場合
別の新しいクラスができて、新しく1人が加わる場合