Stan:「クネクネ」した関係をガウス過程で推定する

Kéry and Royle (2016)の10.14節にある「クネクネ」した(wiggly)関係の推定をStanをつかってガウス過程でやってみました。元ネタでは、在/不在データをさらに検出/不検出するというデータになっていますが、ここではひとまず在/不在の部分だけをあつかっています。

• Kéry M. and Royle J.A. (2016) "Applied Hierarchical Modeling in Ecology" Vol.1, Academic Press

RコードとStanコードです。psi_funは、共変量xに対応した出現率psiを返す関数で、この関係が下の図の赤線のようにクネクネです。このpsiの値に応じて、在/不在Yがベルヌーイ分布により きまります。

Stanコードはマニュアルの18章を参考にしています。X2に予測する共変量の値をいれています。

赤線が、共変量と出現率psiとの関係。赤丸は、乱数で生成された共変量の値に対応するpsiの値。黒丸は、psiから生成された在/不在Yの値。

主なパラメータの事後分布の要約です。sigmaの有効サンプルサイズがちいさめです。

Inference for Stan model: GP_bin.
4 chains, each with iter=6000; warmup=1000; thin=5; 
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
rho   0.82    0.00 0.18 0.50 0.69 0.81 0.93  1.20  2300 1.00
alpha 1.61    0.03 0.48 0.89 1.27 1.56 1.88  2.73   193 1.01
sigma 0.78    0.07 0.49 0.13 0.38 0.69 1.08  1.96    50 1.05

予測値（事後平均）を重ねてプロットしたものです。Kéry and Royle (2016)では、罰則付きスプラインと多項式をそれぞれあてはめているのですが、その両者よりも「クネクネ」に忠実になっているような気がします。

タグ：R STAn