xbinomial
samplingAbstract
Task We test sampling from the custom likelihoods
xbinomial_lpdf
and xbinomial_logit_lpdf
.
xbinomial_lpdf
temp1 <- rstan::stan(file = "xbinom.stan", data = list(m = 10, rho = 0.5), warmup = 500, iter = 1000, refresh = 0)
## Warning: There were 5 divergent transitions after warmup. See
## https://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
## to find out why this is a problem and how to eliminate them.
## Warning: Examine the pairs() plot to diagnose sampling problems
rstan::summary(temp1)$summary
## mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat
## y 5.1715341 0.07317071 1.6431686 2.166292 3.9577373 5.14370274 6.368971 8.393380 504.3009 1.005065
## lp__ -0.3342727 0.02947408 0.7693032 -2.450940 -0.5420835 -0.05227256 0.188336 0.256736 681.2625 1.003222
out1 <- rstan::extract(temp1)
plot(out1$y)
hist(out1$y)
xbinomial_logit_lpdf
temp2 <- rstan::stan(file = "xbinom_logit.stan", data = list(m = 10, eta = 0), warmup = 500, iter = 1000, refresh = 0)
## recompiling to avoid crashing R session
## Warning: There were 11 divergent transitions after warmup. See
## https://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
## to find out why this is a problem and how to eliminate them.
## Warning: Examine the pairs() plot to diagnose sampling problems
## Warning: Tail Effective Samples Size (ESS) is too low, indicating posterior variances and tail quantiles may be unreliable.
## Running the chains for more iterations may help. See
## https://mc-stan.org/misc/warnings.html#tail-ess
rstan::summary(temp2)$summary
## mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat
## y 4.953689 0.07252057 1.5982580 1.832909 3.840170 4.966525 6.016633 8.077665 485.7036 1.005330
## lp__ -1.734385 0.04377748 0.8424262 -4.166279 -1.964789 -1.407077 -1.183856 -1.129246 370.3073 1.006264
out2 <- rstan::extract(temp2)
plot(out2$y)
hist(out2$y)
df1 <- data.frame(x = out1$y)
df2 <- data.frame(x = out2$y)
ks.test(out1$y, out2$y)
## Warning in ks.test.default(out1$y, out2$y): p-value will be approximate in the presence of ties
##
## Asymptotic two-sample Kolmogorov-Smirnov test
##
## data: out1$y and out2$y
## D = 0.0665, p-value = 0.0002883
## alternative hypothesis: two-sided
ggplot() +
stat_ecdf(data = df1, aes(x = x), geom = "step") +
stat_ecdf(data = df2, aes(x = x), geom = "step", col = "red")
xbinom.stan
functions {
real xbinomial_lpdf(real y, real m, real rho) {
return(lchoose(m, y) + y * log(rho) + (m - y) * log(1 - rho));
}
}
data {
real rho;
real m;
}
parameters {
real<lower=0> y;
}
model {
y ~ xbinomial(m, rho);
}
xbinom_logit.stan
functions {
real xbinomial_logit_lpdf(real y, real m, real eta) {
real dens = lchoose(m, y) + y * log(inv_logit(eta)) + (m - y) * log(1 - inv_logit(eta));
real jacobian = log(inv_logit(eta)) + log(inv_logit(-eta)); // Alternative: eta - 2 * log(1 + exp(eta))
return dens + jacobian;
}
}
data {
real eta;
real m;
}
parameters {
real<lower=0> y;
}
model {
y ~ xbinomial_logit(m, eta);
}
sessionInfo()
## R version 4.2.0 (2022-04-22)
## Platform: x86_64-apple-darwin17.0 (64-bit)
## Running under: macOS Monterey 12.6
##
## Matrix products: default
## LAPACK: /Library/Frameworks/R.framework/Versions/4.2/Resources/lib/libRlapack.dylib
##
## locale:
## [1] en_GB.UTF-8/en_GB.UTF-8/en_GB.UTF-8/C/en_GB.UTF-8/en_GB.UTF-8
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] tidyr_1.2.0 magrittr_2.0.3 dplyr_1.0.9 rmarkdown_2.14 bayesplot_1.9.0 rstan_2.21.5 ggplot2_3.3.6 StanHeaders_2.21.0-7
## [9] sf_1.0-7 bsae_0.2.7
##
## loaded via a namespace (and not attached):
## [1] sass_0.4.1 bit64_4.0.5 jsonlite_1.8.0 viridisLite_0.4.0 ids_1.0.1 bslib_0.3.1 RcppParallel_5.1.5 assertthat_0.2.1 askpass_1.1
## [10] highr_0.9 stats4_4.2.0 blob_1.2.3 yaml_2.3.5 pillar_1.7.0 RSQLite_2.2.14 lattice_0.20-45 glue_1.6.2 RcppEigen_0.3.3.9.2
## [19] uuid_1.1-0 digest_0.6.29 colorspace_2.0-3 cowplot_1.1.1 htmltools_0.5.2 Matrix_1.5-1 plyr_1.8.7 pkgconfig_2.0.3 purrr_0.3.4
## [28] scales_1.2.0 processx_3.5.3 tibble_3.1.7 openssl_2.0.1 proxy_0.4-26 generics_0.1.2 farver_2.1.0 ellipsis_0.3.2 cachem_1.0.6
## [37] withr_2.5.0 cli_3.3.0 crayon_1.5.1 memoise_2.0.1 evaluate_0.15 ps_1.7.0 fansi_1.0.3 class_7.3-20 pkgbuild_1.3.1
## [46] tools_4.2.0 loo_2.5.1 prettyunits_1.1.1 lifecycle_1.0.1 matrixStats_0.62.0 stringr_1.4.0 munsell_0.5.0 callr_3.7.0 compiler_4.2.0
## [55] jquerylib_0.1.4 e1071_1.7-9 rlang_1.0.2 classInt_0.4-3 units_0.8-0 grid_4.2.0 ggridges_0.5.3 rstudioapi_0.13 labeling_0.4.2
## [64] orderly_1.4.3 gtable_0.3.0 codetools_0.2-18 inline_0.3.19 DBI_1.1.2 reshape2_1.4.4 R6_2.5.1 gridExtra_2.3 knitr_1.39
## [73] bit_4.0.4 fastmap_1.1.0 utf8_1.2.2 rprojroot_2.0.3 BH_1.78.0-0 KernSmooth_2.23-20 stringi_1.7.6 parallel_4.2.0 Rcpp_1.0.8.3
## [82] vctrs_0.4.1 tidyselect_1.1.2 xfun_0.31