この記事では積分消去する変数が2個ある場合において、階層モデルのWAICを算出します。シンプソンの公式を使う場合とモンテカルロ積分を使う場合の二通り行います。
その前に、積分消去する変数が1個の場合を復習します。以前の記事のモデル2の(2)の予測を考えます。ここではこのモデルをmodel1と呼びます。
シンプソンの公式を使うバージョンのstanコードは以下です。
29～40行目でシンプソンの公式の係数だけ作っておきます。そうすると13～16行目で簡潔に書けます。
WAICの算出には、事前分布に由来する尤度は入れないと理解しているので、2～5行目の関数fで事前分布の項は入れていないです。
モンテカルロ積分を使うバージョンのstanコードは以下です。
45～49行目がモンテカルロ積分をしている部分です。あってると思っていますが、100%の自信はありません。
＊ ＊ ＊
次に標準偏差にもグループ差を考慮したモデルを考えます。このモデルをmodel2と呼びます。この場合、積分消去する変数が2個になります。
シンプソンの公式を使うバージョンのstanコードは以下です。
標準偏差にグループ差が入ったため、2～5行目の対数尤度の関数fが変更されています。
2次元のシンプソンの公式はstackexchangeのこのページを参考にしました。
45行目で2次元のシンプソンの公式の係数を作っています。1次元のシンプソンの公式の係数ベクトルのテンソル積で作るのがラクです。15～20行目で簡潔に書けます。76行目、標準偏差は正の値しか取らないので下限として1e-8を使っています。
モンテカルロ積分を使うバージョンのstanコードは以下です。
51行目で半正規分布からs_rndを生成しています。正規分布から生成した乱数の絶対値を取れば半正規分布から生成した乱数になるので、fabs関数を使っています。
実行するRコードは以下です。model1が真の場合(Data 1)とmodel2が真の場合(Data 2)の両方でWAICを計算しました。
2次元の場合において、シンプソンの公式の分割数Mはおおよそ300あれば十分でした（300 x 300の分割）。モンテカルロ積分のシミュレーション数Mはおおよそ100000あれば十分でした。このあたりは問題によりけりと思います。
個人的には、積分消去する変数が3個以下ならばシンプソンの公式の方をおすすめしたいです。コードをそのまま使う場合や参考にした場合にはできれば引用をお願いしますね。
本記事は伊東さんとの議論がきっかけです。ありがとうございました。
