diff --git a/DESCRIPTION b/DESCRIPTION
index 5315e06..d4ead43 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -8,8 +8,8 @@ Description: Tool-set to support Bayesian evidence synthesis.  This
     for details on applying this package while Neuenschwander et al. (2010)
     <doi:10.1177/1740774509356002> and Schmidli et al. (2014)
     <doi:10.1111/biom.12242> explain details on the methodology.
-Version: 1.7-3
-Date: 2024-01-02
+Version: 1.7-4
+Date: 2024-11-21
 Authors@R: c(person("Novartis", "Pharma AG", role = "cph")
             ,person("Sebastian", "Weber", email="sebastian.weber@novartis.com", role=c("aut", "cre"))
             ,person("Beat", "Neuenschwander", email="beat.neuenschwander@novartis.com", role="ctb")
@@ -25,8 +25,8 @@ Imports:
     methods,
     Rcpp (>= 0.12.0),
     RcppParallel (>= 5.0.1),
-    rstan (>= 2.26.0),
-    rstantools (>= 2.3.1),
+    rstan (>= 2.32.0),
+    rstantools (>= 2.4.0),
     assertthat,
     mvtnorm,
     Formula,
@@ -44,8 +44,8 @@ LinkingTo:
     Rcpp (>= 0.12.0),
     RcppEigen (>= 0.3.3.3.0),
     RcppParallel (>= 5.0.1),
-    rstan (>= 2.26.0),
-    StanHeaders (>= 2.26.0)
+    rstan (>= 2.32.0),
+    StanHeaders (>= 2.32.0)
 License: GPL (>=3)
 LazyData: true
 Biarch: true
@@ -71,4 +71,5 @@ Suggests:
 VignetteBuilder: knitr
 SystemRequirements: GNU make, pandoc (>= 1.12.3), pngquant, C++17
 Encoding: UTF-8
-RoxygenNote: 7.2.3
+RoxygenNote: 7.3.1
+Config/testthat/edition: 3
diff --git a/Makefile b/Makefile
index 2c955c4..ce32c5c 100644
--- a/Makefile
+++ b/Makefile
@@ -13,8 +13,11 @@ PROJROOT_ABS=$(abspath .)
 
 RPKG=$(patsubst ‘%’, %, $(word 2, $(shell grep ^Package: DESCRIPTION)))
 INCS = 
-R_PKG_SRCS = $(wildcard R/*.R)
+R_PKG_SRCS = $(wildcard R/*.R inst/examples/*R)
 R_SRCS = $(wildcard *.R $(foreach fd, $(SRCDIR), $(fd)/*.R))
+R_TEST_SRCS = $(wildcard tests/testthat/test*.R)
+R_TEST_OBJS = $(R_TEST_SRCS:.R=.Rtest)
+R_TESTFAST_OBJS = $(R_TEST_SRCS:.R=.Rtestfast)
 RMD_SRCS = $(wildcard *.Rmd $(foreach fd, $(SRCDIR), $(fd)/x*.Rmd))
 STAN_SRCS = $(wildcard *.stan $(foreach fd, $(SRCDIR), $(fd)/*.stan))
 SRCS = $(R_PKG_SRCS) $(R_SRCS) $(RMD_SRCS) $(STAN_SRCS)
@@ -28,6 +31,7 @@ PKG_VERSION ?= $(patsubst ‘%’, %, $(word 2, $(shell grep ^Version DESCRIPTIO
 GIT_TAG ?= v$(PKG_VERSION)
 
 MD5 ?= md5sum
+TMPDIR := $(realpath $(shell mktemp -d))
 
 all : $(TARGET)
 
@@ -45,31 +49,47 @@ all : $(TARGET)
 	cd $(@D); echo running $(RCMD) -e "rmarkdown::render('$(<F)', output_format=rmarkdown::html_document(self_contained=TRUE))"
 	cd $(@D); $(RCMD) -e "rmarkdown::render('$(<F)', output_format=rmarkdown::html_document(self_contained=TRUE))"
 
-R/stanmodels.R:
+tests/%.Rtest : tests/%.R
+	NOT_CRAN=true $(RCMD) -e "devtools::load_all()" -e "test_file('$<')" > $@ 2>&1
+	@printf "Test summary for $(<F): "
+	@grep '^\[' $@ | tail -n 1
+
+tests/%.Rtestfast : tests/%.R
+	NOT_CRAN=false $(RCMD) -e "devtools::load_all()" -e "test_file('$<')" > $@ 2>&1
+	@printf "Test summary for $(<F): "
+	@grep '^\[' $@ | tail -n 1
+
+
+R/stanmodels.R: $(STAN_SRCS)
 	## ensure that NAMESPACE contains load directive
 	echo "# Generated by roxygen2: do not edit by hand" > NAMESPACE
 	echo "import(Rcpp)" >> NAMESPACE
 	echo "import(methods)" >> NAMESPACE
 	echo "importFrom(rstan, sampling)" >> NAMESPACE
 	echo "useDynLib($(RPKG), .registration = TRUE)" >> NAMESPACE
+	install -d src
 	"${R_HOME}/bin/Rscript" -e "rstantools::rstan_config()"
+	touch R/stanmodels.R
 
-src/package-binary: $(STAN_SRCS) R/stanmodels.R
+src/package-binary: R/stanmodels.R
 	## ensure that NAMESPACE contains load directive
 	echo "# Generated by roxygen2: do not edit by hand" > NAMESPACE
 	echo "import(Rcpp)" >> NAMESPACE
 	echo "import(methods)" >> NAMESPACE
 	echo "importFrom(rstan, sampling)" >> NAMESPACE
 	echo "useDynLib($(RPKG), .registration = TRUE)" >> NAMESPACE
-	"${R_HOME}/bin/R" --slave -e 'library(pkgbuild); pkgbuild::compile_dll()'
+	install -d src
+	"${R_HOME}/bin/Rscript" -e 'pkgbuild::compile_dll(debug=FALSE)'
 	touch src/package-binary
 
 man/package-doc: $(R_PKG_SRCS) $(BIN_OBJS)
-	"${R_HOME}/bin/R" --slave -e 'library(roxygen2); roxygen2::roxygenize()'
+	"${R_HOME}/bin/Rscript" -e 'roxygen2::roxygenize()'
 	touch man/package-doc
 
+inst/sbc/sbc_report.html : inst/sbc/calibration.rds
+
 inst/sbc/calibration.rds :
-	@echo Please run inst/sbc/make_reference_rankhist.R
+	echo "Please run inst/sbc/make_reference_rankhist.R"
 	exit 1
 
 R/sysdata.rda: inst/sbc/calibration.rds
@@ -85,12 +105,21 @@ NAMESPACE: man/package-doc
 
 
 PHONY := $(TARGET)
-$(TARGET): build/r-source
+$(TARGET): build/r-source-fast
 
-build/r-source : $(BIN_OBJS) $(DOC_OBJS) $(SRCS)
+build/r-source-fast : $(BIN_OBJS) $(DOC_OBJS) $(SRCS)
 	install -d build
-	cd build; $(RCMD) CMD build .. --no-build-vignettes --no-manual
-	cd build; mv $(RPKG)_$(PKG_VERSION).tar.gz $(RPKG)-source.tar.gz
+	git archive --format=tar.gz --prefix $(RPKG)-$(GIT_TAG)/ HEAD > build/$(RPKG)-$(GIT_TAG).tar.gz
+	rm -rf build/$(RPKG)-$(GIT_TAG)
+	cd build; tar x -C $(TMPDIR) -f $(RPKG)-$(GIT_TAG).tar.gz
+	rm -f build/$(RPKG)-$(GIT_TAG).tar.gz
+	cp -v NAMESPACE $(TMPDIR)/$(RPKG)-$(GIT_TAG)
+	install -d $(TMPDIR)/$(RPKG)-$(GIT_TAG)/man
+	cp -v man/*.Rd $(TMPDIR)/$(RPKG)-$(GIT_TAG)/man
+	cd $(TMPDIR)/$(RPKG)-$(GIT_TAG); "${R_HOME}/bin/R" --slave --file=tools/make-ds.R
+	cd $(TMPDIR); NOT_CRAN=false "${R_HOME}/bin/R" CMD build $(RPKG)-$(GIT_TAG) --no-build-vignettes --no-manual
+	rm -rf $(TMPDIR)/$(RPKG)-$(GIT_TAG)
+	mv $(TMPDIR)/$(RPKG)_$(PKG_VERSION).tar.gz build/$(RPKG)-source.tar.gz
 	touch build/r-source-fast
 
 build/r-source-release : $(BIN_OBJS) $(DOC_OBJS) $(SRCS) inst/sbc/sbc_report.html
@@ -124,8 +153,33 @@ derived : NAMESPACE $(BIN_OBJS) $(DOC_OBJS)
 
 PHONY += r-source-check
 r-source-check : r-source
-	cd build; tar xvzf RBesT-source.tar.gz
-	cd build; NOT_CRAN=true $(RCMD) CMD check RBesT
+	cd build; tar xvzf $(RPKG)-source.tar.gz
+	cd build; NOT_CRAN=true $(RCMD) CMD check $(RPKG)
+
+PHONY += r-source-release-check
+r-source-release-check : r-source-release
+	cd build; tar xvzf $(RPKG)_$(PKG_VERSION).tar.gz
+	cd build; NOT_CRAN=true $(RCMD) CMD check $(RPKG)
+
+build/installed/$(RPKG)/DESCRIPTION : build/r-source-fast
+	rm -rf build/installed
+	install -d build/installed
+	cd build; $(RCMD) CMD INSTALL --library=./installed --no-docs --no-multiarch --no-test-load --no-clean-on-error $(RPKG)-source.tar.gz
+
+PHONY += dev-install
+dev-install: build/installed/$(RPKG)/DESCRIPTION
+
+PHONY += test-all
+test-all : $(R_TEST_OBJS)
+
+PHONY += testfast-all
+testfast-all : $(R_TESTFAST_OBJS)
+
+PHONY += retestfast-all
+retestfast-all : clean-test $(R_TESTFAST_OBJS)
+
+PHONY += retest-all
+retest-all : clean-test $(R_TEST_OBJS)
 
 #$(DIR_OBJ)/%.o: %.c $(INCS)
 #    mkdir -p $(@D)
@@ -146,6 +200,14 @@ clean:
 	rm -f vignettes/*.html
 	rm -f vignettes/*.docx
 	rm -rf .Rd2pdf*
+	rm -f $(R_TEST_OBJS)
+	rm -f $(R_TESTFAST_OBJS)
+	rm -rf src
+	rm -f R/stanmodels.R
+
+clean-test:
+	rm -f $(R_TEST_OBJS)
+	rm -f $(R_TESTFAST_OBJS)
 
 PHONY += doc
 doc: $(DOC_OBJS)
diff --git a/NAMESPACE b/NAMESPACE
index 377f92c..103bfa0 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -76,6 +76,7 @@ S3method(print,betaBinomialMix)
 S3method(print,betaMix)
 S3method(print,decision1S)
 S3method(print,decision2S)
+S3method(print,dlink)
 S3method(print,gMAP)
 S3method(print,gMAPpred)
 S3method(print,gMAPsummary)
diff --git a/NEWS.md b/NEWS.md
index f08c827..bcc1ea5 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -1,3 +1,26 @@
+# RBesT 1.7-4 - November 21st, 2024
+
+## Enhancements
+
+* Added for `mixstanvar` automatic generation of distribution
+  functions allowing truncated priors in `brms` with mixtures
+* Slight speed increase for Stan model by more efficient construction
+  of likelihood. Also reducing object size of `gMAP` objects by
+  avoiding to store redundant variables.
+* Upgraded `testthat` edition to version 3.
+
+## Bugfixes
+
+* Fix issue #18 of ESS ELIR aborting whenever one mixture component
+  has zero weight.
+* Fixed a rare issue when estimating ESS ELIR for beta mixtures. The
+  calculation aborted due to boundary issues which are now avoided.
+* Avoid over and underflow of beta-binomial distribution.
+* Replace calls to deprecated `ggplot2::qplot` with respective calls
+  to `ggplot2::ggplot`.
+* Use new `array` notation for `mixstanvar` generated Stan code to
+  make generated Stan programs compatible with Stan >= 2.33
+
 # RBesT 1.7-3 - January 2nd, 2024
 
 ## Enhancements
diff --git a/R/RBesT-package.R b/R/RBesT-package.R
index d128d85..14cc304 100644
--- a/R/RBesT-package.R
+++ b/R/RBesT-package.R
@@ -38,6 +38,7 @@
 #' \tab \code{rel.tol=.Machine$double.eps^0.25,} \tab \cr
 #' \tab \code{abs.tol=.Machine$double.eps^0.25,} \tab \cr
 #' \tab \code{subdivisions=1E3)} \tab \cr
+#' \code{RBesT.integrate_prob_eps} \tab \code{1E-6} \tab probability mass left out from tails if integration needs to be restricted in range \cr
 #' }
 #'
 #' @section Version History:
diff --git a/R/dlink.R b/R/dlink.R
index 19d6653..7f1aac4 100644
--- a/R/dlink.R
+++ b/R/dlink.R
@@ -97,5 +97,6 @@ is.dlink <- function(x)
 is.dlink_identity <- function(x)
     is.dlink(x) & x$name == "identity"
 
+#' @export
 print.dlink <- function(x, ...)
     cat("Link:", x$name, "\n")
diff --git a/R/integrate_logit_log.R b/R/integrate_logit_log.R
index 651be98..653ff4f 100644
--- a/R/integrate_logit_log.R
+++ b/R/integrate_logit_log.R
@@ -16,8 +16,8 @@
 #' @param Lpupper logit of upper cumulative density
 #'
 #' @keywords internal
-integrate_density_log <- function(log_integrand, mix, Lplower=-Inf, Lpupper=Inf, eps=1E-9) {
-    .integrand_comp <- function(mix_comp) {
+integrate_density_log <- function(log_integrand, mix, Lplower=-Inf, Lpupper=Inf, eps=getOption("RBesT.integrate_prob_eps", 1E-6)) {
+    .integrand_comp_logit <- function(mix_comp) {
         function(l) {
             u   <- inv_logit(l)
             lp  <- log_inv_logit(l)
@@ -37,10 +37,9 @@ integrate_density_log <- function(log_integrand, mix, Lplower=-Inf, Lpupper=Inf,
     upper  <- inv_logit(Lpupper)
     return(sum(vapply(1:Nc, function(comp) {
         mix_comp  <- mix[[comp, rescale=TRUE]]
-        mix_comp_identity <- mix_comp
-        dlink(mix_comp_identity) <- identity_dlink
-        if (all(dmix(mix_comp_identity, support(mix_comp_identity)) == 0 )) {
-            return(.integrate(.integrand_comp(mix_comp), Lplower, Lpupper))
+        fn_integrand_comp_logit <- .integrand_comp_logit(mix_comp)
+        if (all(!is.na(fn_integrand_comp_logit(c(Lplower, Lpupper))))) {
+            return(.integrate(fn_integrand_comp_logit, Lplower, Lpupper))
         }
         lower_comp  <- ifelse(Lplower==-Inf, qmix(mix_comp, eps), qmix(mix_comp, lower))
         upper_comp  <- ifelse(Lpupper==Inf, qmix(mix_comp, 1-eps), qmix(mix_comp, upper))
@@ -48,8 +47,8 @@ integrate_density_log <- function(log_integrand, mix, Lplower=-Inf, Lpupper=Inf,
     }, c(0.1)) * mix[1,]))
 }
 
-integrate_density <- function(integrand, mix, Lplower=-Inf, Lpupper=Inf, eps=1E-9) {
-    .integrand_comp <- function(mix_comp) {
+integrate_density <- function(integrand, mix, Lplower=-Inf, Lpupper=Inf, eps=getOption("RBesT.integrate_prob_eps", 1E-6)) {
+    .integrand_comp_logit <- function(mix_comp) {
         function(l) {
             u   <- inv_logit(l)
             lp  <- log_inv_logit(l)
@@ -57,18 +56,19 @@ integrate_density <- function(integrand, mix, Lplower=-Inf, Lpupper=Inf, eps=1E-
             exp(lp + lnp) * integrand(qmix(mix_comp, u))
         }
     }
-
     Nc <- ncol(mix)
 
     lower  <- inv_logit(Lplower)
     upper  <- inv_logit(Lpupper)
+    
     return(sum(vapply(1:Nc, function(comp) {
         mix_comp  <- mix[[comp, rescale=TRUE]]
-        mix_comp_identity <- mix_comp
-        dlink(mix_comp_identity) <- identity_dlink
-        if (all(dmix(mix_comp_identity, support(mix_comp_identity)) == 0 )) {
-            return(.integrate(.integrand_comp(mix_comp), Lplower, Lpupper))
+        ## ensure that the integrand is defined at the boundaries...
+        fn_integrand_comp_logit <- .integrand_comp_logit(mix_comp)
+        if (all(!is.na(fn_integrand_comp_logit(c(Lplower, Lpupper))))) {
+            return(.integrate(fn_integrand_comp_logit, Lplower, Lpupper))
         }
+        ## ... otherwise we avoid the boundaries by eps prob density:
         lower_comp  <- ifelse(Lplower==-Inf, qmix(mix_comp, eps), qmix(mix_comp, lower))
         upper_comp  <- ifelse(Lpupper==Inf, qmix(mix_comp, 1-eps), qmix(mix_comp, upper))
         return(.integrate(function(x) integrand(x) * dmix(mix_comp, x), lower_comp, upper_comp))
@@ -80,7 +80,8 @@ integrate_density <- function(integrand, mix, Lplower=-Inf, Lpupper=Inf, eps=1E-
     args <- modifyList(list(lower=lower, upper=upper,
                             rel.tol=.Machine$double.eps^0.25,
                             abs.tol=.Machine$double.eps^0.25,
-                            subdivisions=1000),
+                            subdivisions=1000,
+                            stop.on.error=TRUE),
                        integrate_args_user)
 
     integrate(integrand,
@@ -88,5 +89,6 @@ integrate_density <- function(integrand, mix, Lplower=-Inf, Lpupper=Inf, eps=1E-
               upper=args$upper,
               rel.tol=args$rel.tol,
               abs.tol=args$abs.tol,
-              subdivisions=args$subdivisions)$value
+              subdivisions=args$subdivisions,
+              stop.on.error=args$stop.on.error)$value
 }
diff --git a/R/mix.R b/R/mix.R
index 1ac1da1..5388492 100644
--- a/R/mix.R
+++ b/R/mix.R
@@ -222,20 +222,32 @@ pmix.normMix  <- function(mix, q, lower.tail = TRUE, log.p=FALSE) pmix_impl(pnor
 ## pmix.betaBinomialMix <- function(mix, q, lower.tail = TRUE, log.p=FALSE) pmix_impl(Curry(pBetaBinomial, n=attr(mix, "n")), mix, q, lower.tail, log.p)
 pmix.betaBinomialMix <- function(mix, q, lower.tail = TRUE, log.p=FALSE) {
     assert_that(is.dlink_identity(attr(mix, "link")))
-##     ## the analytic solution needs the generalized hypergeometric
-##     ## function which is only available in the hyper-geo package which
-##     ## is why a numeric integration is performed here
-##     ## treat out-of-bounds quantiles special
-     out_low  <- q<0
-     out_high <- q>attr(mix, "n")
-     q[out_low | out_high] <- 0
-     dist <- cumsum(dmix.betaBinomialMix(mix, seq(0,max(q))))
-     p <- dist[q+1]
-     p[out_low] <- 0
-     p[out_high] <- 1
-     if(!lower.tail) p <- 1-p
-     if(log.p) p <- log(p)
-     return(p)
+    ##     ## the analytic solution needs the generalized hypergeometric
+    ##     ## function which is only available in the hyper-geo package which
+    ##     ## is why a numeric integration is performed here
+    ##     ## treat out-of-bounds quantiles special
+##    if (requireNamespace("extraDistr", quietly = TRUE)) {
+##        Nq <- length(q)
+##        size <- attr(mix, "n")
+##        if(!log.p)
+##            return(matrixStats::colSums2(matrix(mix[1,] * extraDistr::pbbinom(rep(q, each=Nc), size, rep(mix[2,], times=Nq), rep(mix[3,], times=Nq), lower.tail=lower.tail, log.p=FALSE), nrow=Nc)))
+        ## log version is slower, but numerically more stable
+##        res <- matrixStats::colLogSumExps(matrix(log(mix[1,]) + dist(rep(q, each=Nc), rep(mix[2,], times=Nq), rep(mix[3,], times=Nq), lower.tail=lower.tail, log.p=TRUE), nrow=Nc))
+##        return(res)
+##    }
+    
+    ## default implementation uses brute-force integration
+    out_low  <- q<0
+    out_high <- q>attr(mix, "n")
+    q[out_low | out_high] <- 0
+    dist <- cumsum(dmix.betaBinomialMix(mix, seq(0,max(q))))
+    dist <- pmin(pmax(dist, 0), 1) ## avoid over and underflows
+    p <- dist[q+1]
+    p[out_low] <- 0
+    p[out_high] <- 1
+    if(!lower.tail) p <- 1-p
+    if(log.p) p <- log(p)
+    return(p)
 }
 
 ## internal redefinition of negative binomial
diff --git a/R/mixcombine.R b/R/mixcombine.R
index 43033c9..c1e5310 100644
--- a/R/mixcombine.R
+++ b/R/mixcombine.R
@@ -19,14 +19,8 @@
 #' @family mixdist
 #' @seealso \code{\link{robustify}}
 #'
-#' @examples
-#' # beta with two informative components
-#' bm <- mixbeta(inf=c(0.5, 10, 100), inf2=c(0.5, 30, 80))
-#'
-#' # robustified with mixcombine, i.e. a 10% uninformative part added
-#' unif <- mixbeta(rob=c(1,1,1))
-#' mixcombine(bm, unif, weight=c(9, 1))
-#'
+#' @example inst/examples/mixcombine.R
+#' 
 #' @export
 mixcombine <- function(..., weight, rescale=TRUE) {
     comp <- list(...)
diff --git a/R/mixess.R b/R/mixess.R
index 829d4b1..03ae16a 100644
--- a/R/mixess.R
+++ b/R/mixess.R
@@ -57,62 +57,8 @@
 #' Predictively Consistent Prior Effective Sample Sizes.
 #' \emph{pre-print} 2019; arXiv:1907.04185
 #'
-#' @examples
-#' # Conjugate Beta example
-#' a <- 5
-#' b <- 15
-#' prior <- mixbeta(c(1, a, b))
-#'
-#' ess(prior)
-#' (a+b)
-#'
-#' # Beta mixture example
-#' bmix <- mixbeta(rob=c(0.2, 1, 1), inf=c(0.8, 10, 2))
-#'
-#' ess(bmix, "elir")
-#'
-#' ess(bmix, "moment")
-#' # moments method is equivalent to
-#' # first calculate moments
-#' bmix_sum <- summary(bmix)
-#' # then calculate a and b of a matching beta
-#' ab_matched <- ms2beta(bmix_sum["mean"], bmix_sum["sd"])
-#' # finally take the sum of a and b which are equivalent
-#' # to number of responders/non-responders respectivley
-#' round(sum(ab_matched))
-#'
-#' ess(bmix, method="morita")
-#'
-#' # Predictive consistency of elir
-#' \donttest{
-#' n_forward <- 1E2
-#' bmixPred <- preddist(bmix, n=n_forward)
-#' pred_samp <- rmix(bmixPred, 1E3)
-#' pred_ess <- sapply(pred_samp, function(r) ess(postmix(bmix, r=r, n=n_forward), "elir") )
-#' ess(bmix, "elir")
-#' mean(pred_ess) - n_forward
-#' }
-#'
-#' # Normal mixture example
-#' nmix <- mixnorm(rob=c(0.5, 0, 2), inf=c(0.5, 3, 4), sigma=10)
-#'
-#' ess(nmix, "elir")
-#'
-#' ess(nmix, "moment")
-#'
-#' ## the reference scale determines the ESS
-#' sigma(nmix) <- 20
-#' ess(nmix)
-#'
-#' # Gamma mixture example
-#' gmix <- mixgamma(rob=c(0.3, 20, 4), inf=c(0.7, 50, 10))
-#'
-#' ess(gmix) ## interpreted as appropriate for a Poisson likelihood (default)
-#'
-#' likelihood(gmix) <- "exp"
-#' ess(gmix) ## interpreted as appropriate for an exponential likelihood
-#'
-#'
+#' @example inst/examples/ess.R
+#' 
 #' @export
 ess <- function(mix, method=c("elir", "moment", "morita"), ...) UseMethod("ess")
 #' @export
@@ -156,12 +102,12 @@ mixInfo <- function(mix, x, dens, gradl, hessl) {
     dhl <- (hessl(x,a,b) + dgl^2)
     ## attempt numerically more stable log calculations if possible,
     ## i.e. if all sings are the same
-    if(all(dgl < 0) || all(dgl > 0)) {
+    if(all(!is.na(dgl)) && ( all(dgl < 0) || all(dgl > 0))) {
         gsum <- exp(2*matrixStats::logSumExp(lwdensComp + log(abs(dgl))))
     } else {
         gsum <- (sum(exp(lwdensComp)*dgl))^2
     }
-    if(all(dhl < 0) || all(dhl > 0)) {
+    if(all(!is.na(dhl)) && (all(dhl < 0) || all(dhl > 0))) {
         hsum <- sign(dhl[1]) * exp(matrixStats::logSumExp(lwdensComp + log(abs(dhl))))
     } else {
         hsum <- (sum(exp(lwdensComp)*dhl))
diff --git a/R/mixplot.R b/R/mixplot.R
index 974098e..972ef1a 100644
--- a/R/mixplot.R
+++ b/R/mixplot.R
@@ -12,7 +12,7 @@
 #' which will display colour-coded each mixture component of the
 #' density in addition to the density.
 #' @param size controls the linesize in plots.
-#' @param ... extra arguments passed on to the \code{\link[ggplot2]{qplot}} call.
+#' @param ... extra arguments passed on to the plotted function.
 #'
 #' @details Plot function for mixture distribution objects. It shows
 #' the density/quantile/cumulative distribution (corresponds to
@@ -85,7 +85,7 @@ plot.mix <- function(x, prob=0.99, fun=dmix, log=FALSE, comp=TRUE, size=1.25, ..
         plot_fun <- function(x, ...) fun(floor(x), ...)
         plot_geom <- "step"
     } else {
-        plot_fun <- fun
+        plot_fun <- function(x, ...) fun(x, ...)
         plot_geom <- "line"
     }
     n_fun <- 501
diff --git a/R/mixstanvar.R b/R/mixstanvar.R
index db1756d..6da68e8 100644
--- a/R/mixstanvar.R
+++ b/R/mixstanvar.R
@@ -136,7 +136,7 @@ mixstanvar <- function(..., verbose=FALSE) {
 
     if(includes_density("mvnormMix")) {
         sv <- sv + brms::stanvar(name="mixmvnorm_lpdf", scode="
-real mixmvnorm_lpdf(vector y, vector w, vector[] m, matrix[] L) {
+real mixmvnorm_lpdf(vector y, vector w, array[] vector m, array[] matrix L) {
   int Nc = rows(w);
   vector[Nc] lp_mix;
   for(i in 1:Nc) {
@@ -154,7 +154,32 @@ real {{mixdens}}_lpdf(real y, vector w, vector {{arg1}}, vector {{arg2}}) {
      lp_mix[i] = {{standens}}_lpdf(y | {{arg1}}[i], {{arg2}}[i]);
   }
   return log_sum_exp(log(w) + lp_mix);
-}", .open="{{", .close="}}"), block="functions")
+}
+real {{mixdens}}_lcdf(real y, vector w, vector {{arg1}}, vector {{arg2}}) {
+  int Nc = rows(w);
+  vector[Nc] lp_mix;
+  for(i in 1:Nc) {
+     lp_mix[i] = {{standens}}_lcdf(y | {{arg1}}[i], {{arg2}}[i]);
+  }
+  return log_sum_exp(log(w) + lp_mix);
+}
+real {{mixdens}}_lccdf(real y, vector w, vector {{arg1}}, vector {{arg2}}) {
+  int Nc = rows(w);
+  vector[Nc] lp_mix;
+  for(i in 1:Nc) {
+     lp_mix[i] = {{standens}}_lccdf(y | {{arg1}}[i], {{arg2}}[i]);
+  }
+  return log_sum_exp(log(w) + lp_mix);
+}
+real {{mixdens}}_cdf(real y, vector w, vector {{arg1}}, vector {{arg2}}) {
+  int Nc = rows(w);
+  vector[Nc] p_mix;
+  for(i in 1:Nc) {
+     p_mix[i] = {{standens}}_cdf(y | {{arg1}}[i], {{arg2}}[i]);
+  }
+  return sum(w + p_mix);
+}
+", .open="{{", .close="}}"), block="functions")
     }
     if(includes_density("normMix")) {
         sv <- sv + mixdensity_template("mixnorm", "normal", "m", "s")
@@ -193,10 +218,10 @@ mix2brms.mvnormMix <- function(mix, name, verbose=FALSE) {
     mvprior <-  brms::stanvar(Nc, glue::glue("{prefix}Nc")) +
         brms::stanvar(p, glue::glue("{prefix}p")) +
         brms::stanvar(array(mix[1,,drop=TRUE], dim=Nc), glue::glue("{prefix}w"), scode=glue::glue("vector[{prefix}Nc] {prefix}w;")) +
-        brms::stanvar(t(mix[2:(p+1),,drop=FALSE]), glue::glue("{prefix}m"), scode=glue::glue("vector[{prefix}p] {prefix}m[{prefix}Nc];")) +
-        brms::stanvar(Sigma, glue::glue("{prefix}sigma"), scode=glue::glue("matrix[{prefix}p, {prefix}p] {prefix}sigma[{prefix}Nc];")) +
+        brms::stanvar(t(mix[2:(p+1),,drop=FALSE]), glue::glue("{prefix}m"), scode=glue::glue("array[{prefix}Nc] vector[{prefix}p] {prefix}m;")) +
+        brms::stanvar(Sigma, glue::glue("{prefix}sigma"), scode=glue::glue("array[{prefix}Nc] matrix[{prefix}p, {prefix}p] {prefix}sigma;")) +
         brms::stanvar(scode=glue::glue("
-matrix[{{prefix}}p, {{prefix}}p] {{prefix}}sigma_L[{{prefix}}Nc];
+array[{{prefix}}Nc] matrix[{{prefix}}p, {{prefix}}p] {{prefix}}sigma_L;
 for (i in 1:{{prefix}}Nc) {
     {{prefix}}sigma_L[i] = cholesky_decompose({{prefix}}sigma[i]);
 }", .open="{{", .close="}}"), block="tdata")
diff --git a/R/plot_gMAP.R b/R/plot_gMAP.R
index 73d4c5a..a9a9d83 100644
--- a/R/plot_gMAP.R
+++ b/R/plot_gMAP.R
@@ -11,7 +11,7 @@
 #'
 #' @template plot-help
 #' 
-#' @return The function returns a list of \code{\link{ggplot}}
+#' @return The function returns a list of \code{\link[ggplot2:ggplot]{ggplot}}
 #' objects.
 #'
 #' @method plot gMAP
diff --git a/R/predict_gMAP.R b/R/predict_gMAP.R
index c7cd732..b9f19c4 100644
--- a/R/predict_gMAP.R
+++ b/R/predict_gMAP.R
@@ -23,41 +23,7 @@
 #' @seealso \code{\link{gMAP}}, \code{\link{predict.glm}}
 #'
 #' @template example-start
-#' @examples
-#' # create a fake data set with a covariate
-#' trans_cov <- transform(transplant, country=cut(1:11, c(0,5,8,Inf), c("CH", "US", "DE")))
-#' set.seed(34246)
-#' map <- gMAP(cbind(r, n-r) ~ 1 + country | study,
-#'             data=trans_cov,
-#'             tau.dist="HalfNormal",
-#'             tau.prior=1,
-#'             # Note on priors: we make the overall intercept weakly-informative
-#'             # and the regression coefficients must have tighter sd as these are
-#'             # deviations in the default contrast parametrization
-#'             beta.prior=rbind(c(0,2), c(0,1), c(0,1)),
-#'             family=binomial,
-#'             ## ensure fast example runtime
-#'             thin=1, chains=1)
-#'
-#' # posterior predictive distribution for each input data item (shrinkage estimates)
-#' pred_cov <- predict(map)
-#' pred_cov
-#'
-#' # extract sample as matrix
-#' samp <- as.matrix(pred_cov)
-#'
-#' # predictive distribution for each input data item (if the input studies were new ones)
-#' pred_cov_pred <- predict(map, trans_cov)
-#' pred_cov_pred
-#'
-#'
-#' # a summary function returns the results as matrix
-#' summary(pred_cov)
-#'
-#' # obtain a prediction for new data with specific covariates
-#' pred_new <- predict(map, data.frame(country="CH", study=12))
-#' pred_new
-#'
+#' @example inst/examples/predict_gMAP.R
 #' @template example-stop
 #' @rdname predict.gMAP
 #' @method predict gMAP
diff --git a/R/sysdata.rda b/R/sysdata.rda
index b219798..f58c2da 100644
Binary files a/R/sysdata.rda and b/R/sysdata.rda differ
diff --git a/R/zzz.R b/R/zzz.R
index 1a480f6..083a9df 100644
--- a/R/zzz.R
+++ b/R/zzz.R
@@ -8,3 +8,4 @@
     ver <- utils::packageVersion("RBesT")
     packageStartupMessage("This is RBesT version ", ver, " (released ", format(pkg_create_date, "%F"), ", git-sha ", pkg_sha, ")")
 }
+
diff --git a/configure b/configure
index 7e23937..0304fc5 100755
--- a/configure
+++ b/configure
@@ -3,4 +3,3 @@
 # Generated by rstantools.  Do not edit by hand.
 
 "${R_HOME}/bin/Rscript" -e "rstantools::rstan_config()"
-
diff --git a/data/AS.rda b/data/AS.rda
index 290f944..97965a4 100644
Binary files a/data/AS.rda and b/data/AS.rda differ
diff --git a/data/colitis.rda b/data/colitis.rda
index 7af25b0..937189f 100644
Binary files a/data/colitis.rda and b/data/colitis.rda differ
diff --git a/data/crohn.rda b/data/crohn.rda
index 746c481..de50e5d 100644
Binary files a/data/crohn.rda and b/data/crohn.rda differ
diff --git a/data/transplant.rda b/data/transplant.rda
index 3ca48e6..0c5c615 100644
Binary files a/data/transplant.rda and b/data/transplant.rda differ
diff --git a/demo/robustMAP.R b/demo/robustMAP.R
index 2f7e040..5b4e0da 100644
--- a/demo/robustMAP.R
+++ b/demo/robustMAP.R
@@ -170,15 +170,14 @@ kable(ocAdaptSamp, digits=1, caption="Sample size")
 #' ## Additional power Figure under varying pc
 #'
 
-qplot(pt-pc, power, data=power, geom=c("line"), colour=prior) +
+ggplot(power, aes(pt-pc, power, colour=prior)) + geom_line()  +
     facet_wrap(~pc) +
     scale_y_continuous(breaks=seq(0,1,by=0.2)) +
-        scale_x_continuous(breaks=seq(-0.8,0.8,by=0.2)) +
-            coord_cartesian(xlim=c(-0.15,0.5)) +
-            geom_hline(yintercept=0.025, linetype=2) +
-                geom_hline(yintercept=0.8, linetype=2) +
-                    ggtitle("Prob. for alternative for different pc")
-
+    scale_x_continuous(breaks=seq(-0.8,0.8,by=0.2)) +
+    coord_cartesian(xlim=c(-0.15,0.5)) +
+    geom_hline(yintercept=0.025, linetype=2) +
+    geom_hline(yintercept=0.8, linetype=2) +
+    ggtitle("Prob. for alternative for different pc")
 
 #'
 #' ## Bias and rMSE, Figure 1
@@ -265,10 +264,8 @@ est <- foreach(case=names(map), .combine=rbind) %do% {
 }
 
 
-
-qplot(p, 100*bias, data=est, geom="line", colour=prior, main="Bias")
-qplot(p, 100*rMSE, data=est, geom="line", colour=prior, main="rMSE")
-
+ggplot(est, aes(p, 100*bias, colour=prior)) + geom_line() + ggtitle("Bias")
+ggplot(est, aes(p, 100*rMSE, colour=prior)) + geom_line() + ggtitle("rMSE")
 
 
 #'
@@ -330,9 +327,8 @@ post <- foreach(prior=names(mapCol), .combine=rbind) %do% {
     res
 }
 
-qplot(r, mean, data=post, colour=prior, shape=prior) + geom_abline(slope=1/20)
-qplot(r, sd, data=post, colour=prior, shape=prior) + coord_cartesian(ylim=c(0,0.17))
-
+ggplot(post, aes(r, mean, colour=prior, shape=prior)) + geom_point() + geom_abline(slope=1/20)
+ggplot(post, aes(r, sd, colour=prior, shape=prior)) + geom_point() + coord_cartesian(ylim=c(0,0.17))
 
 sessionInfo()
 
diff --git a/inst/examples/ess.R b/inst/examples/ess.R
new file mode 100644
index 0000000..dfe433e
--- /dev/null
+++ b/inst/examples/ess.R
@@ -0,0 +1,55 @@
+# Conjugate Beta example
+a <- 5
+b <- 15
+prior <- mixbeta(c(1, a, b))
+
+ess(prior)
+(a+b)
+
+# Beta mixture example
+bmix <- mixbeta(rob=c(0.2, 1, 1), inf=c(0.8, 10, 2))
+
+ess(bmix, "elir")
+
+ess(bmix, "moment")
+# moments method is equivalent to
+# first calculate moments
+bmix_sum <- summary(bmix)
+# then calculate a and b of a matching beta
+ab_matched <- ms2beta(bmix_sum["mean"], bmix_sum["sd"])
+# finally take the sum of a and b which are equivalent
+# to number of responders/non-responders respectivley
+round(sum(ab_matched))
+
+ess(bmix, method="morita")
+
+# Predictive consistency of elir
+n_forward <- 1E1
+bmixPred <- preddist(bmix, n=n_forward)
+pred_samp <- rmix(bmixPred, 1E2)
+# use more samples here for greater accuracy
+# pred_samp <- rmix(bmixPred, 1E3) 
+pred_ess <- sapply(pred_samp, function(r) ess(postmix(bmix, r=r, n=n_forward), "elir") )
+ess(bmix, "elir")
+mean(pred_ess) - n_forward
+
+# Normal mixture example
+nmix <- mixnorm(rob=c(0.5, 0, 2), inf=c(0.5, 3, 4), sigma=10)
+
+ess(nmix, "elir")
+
+ess(nmix, "moment")
+
+## the reference scale determines the ESS
+sigma(nmix) <- 20
+ess(nmix)
+
+# Gamma mixture example
+gmix <- mixgamma(rob=c(0.3, 20, 4), inf=c(0.7, 50, 10))
+
+ess(gmix) ## interpreted as appropriate for a Poisson likelihood (default)
+
+likelihood(gmix) <- "exp"
+ess(gmix) ## interpreted as appropriate for an exponential likelihood
+
+
diff --git a/inst/examples/mixcombine.R b/inst/examples/mixcombine.R
new file mode 100644
index 0000000..f2ca36f
--- /dev/null
+++ b/inst/examples/mixcombine.R
@@ -0,0 +1,6 @@
+# beta with two informative components
+bm <- mixbeta(inf=c(0.5, 10, 100), inf2=c(0.5, 30, 80))
+
+# robustified with mixcombine, i.e. a 10% uninformative part added
+unif <- mixbeta(rob=c(1,1,1))
+mixcombine(bm, unif, weight=c(9, 1))
diff --git a/inst/examples/predict_gMAP.R b/inst/examples/predict_gMAP.R
new file mode 100644
index 0000000..56dc59c
--- /dev/null
+++ b/inst/examples/predict_gMAP.R
@@ -0,0 +1,33 @@
+# create a fake data set with a covariate
+trans_cov <- transform(transplant, country=cut(1:11, c(0,5,8,Inf), c("CH", "US", "DE")))
+set.seed(34246)
+map <- gMAP(cbind(r, n-r) ~ 1 + country | study,
+            data=trans_cov,
+            tau.dist="HalfNormal",
+            tau.prior=1,
+            # Note on priors: we make the overall intercept weakly-informative
+            # and the regression coefficients must have tighter sd as these are
+            # deviations in the default contrast parametrization
+            beta.prior=rbind(c(0,2), c(0,1), c(0,1)),
+            family=binomial,
+            ## ensure fast example runtime
+            thin=1, chains=1)
+
+# posterior predictive distribution for each input data item (shrinkage estimates)
+pred_cov <- predict(map)
+pred_cov
+
+# extract sample as matrix
+samp <- as.matrix(pred_cov)
+
+# predictive distribution for each input data item (if the input studies were new ones)
+pred_cov_pred <- predict(map, trans_cov)
+pred_cov_pred
+
+
+# a summary function returns the results as matrix
+summary(pred_cov)
+
+# obtain a prediction for new data with specific covariates
+pred_new <- predict(map, data.frame(country="CH", study=12))
+pred_new
diff --git a/inst/examples/robustify.R b/inst/examples/robustify.R
new file mode 100644
index 0000000..3b4c40e
--- /dev/null
+++ b/inst/examples/robustify.R
@@ -0,0 +1,22 @@
+bmix <- mixbeta(inf1=c(0.2, 8, 3), inf2=c(0.8, 10, 2))
+plot(bmix)
+
+rbmix <- robustify(bmix, weight=0.1, mean=0.5)
+rbmix
+plot(rbmix)
+
+
+gmnMix <- mixgamma(inf1=c(0.2, 2, 3), inf2=c(0.8, 2, 5), param="mn")
+plot(gmnMix)
+
+rgmnMix <- robustify(gmnMix, weight=0.1, mean=2)
+rgmnMix
+plot(rgmnMix)
+
+
+nm <- mixnorm(inf1=c(0.2, 0.5, 0.7), inf2=c(0.8, 2, 1), sigma=2)
+plot(nm)
+
+rnMix <- robustify(nm, weight=0.1, mean=0, sigma=2)
+rnMix
+plot(rnMix)
diff --git a/inst/sbc/calibration.md5 b/inst/sbc/calibration.md5
index fd88f2f..90b6094 100644
--- a/inst/sbc/calibration.md5
+++ b/inst/sbc/calibration.md5
@@ -1,3 +1,3 @@
-Created:  2023-08-03 13:56:36 UTC
-git hash: 47e501da2dc04c25171cd0a1b14da6cf5115b0ed
-MD5:      067f0aaa153b8c6ee2d100b34fc7dd99
+Created:  2024-11-21 10:30:15 UTC
+git hash: 5f5afb101f98533e4229f9dc267b3540b172c573
+MD5:      c8f566f591c59887653bac53ab03bb45
diff --git a/inst/sbc/calibration.rds b/inst/sbc/calibration.rds
index 8447e65..2fc287a 100644
Binary files a/inst/sbc/calibration.rds and b/inst/sbc/calibration.rds differ
diff --git a/inst/sbc/make_reference_rankhist.R b/inst/sbc/make_reference_rankhist.R
index 91c46d1..20fc4e3 100755
--- a/inst/sbc/make_reference_rankhist.R
+++ b/inst/sbc/make_reference_rankhist.R
@@ -5,20 +5,29 @@ start_time  <- Sys.time()
 here::i_am("inst/sbc/make_reference_rankhist.R")
 library(here)
 
+rbest_lib_dir <- here("build", "installed")
+rbest_source_dir <- here()
+
+## ensure that the current dev version of OncoBayes2 is build and
+## loaded.
+## code to use if one wants the more conservative load_OB2_dev_install routine to be used
 setwd(here())
-system("make binary")
+system("make dev-install", ignore.stdout=TRUE, ignore.stderr=TRUE)
 setwd(here("inst", "sbc"))
 
 pkg <- c("assertthat", "rstan", "mvtnorm", "checkmate", "Formula", "abind", "dplyr", "tidyr", "here", "bayesplot")
 sapply(pkg, require, character.only=TRUE)
 
-
 library(clustermq)
 library(data.table)
 library(knitr)
 sbc_tools <- new.env()
 source(here("inst", "sbc", "sbc_tools.R"), local=sbc_tools)
-set.seed(453453)
+set.seed(45348346)
+
+sbc_tools$rbest_lib_dir <- rbest_lib_dir
+sbc_tools$rbest_source_dir <- rbest_source_dir
+sbc_tools$load_rbest_dev()
 
 scheduler <- getOption("clustermq.scheduler")
 
@@ -71,7 +80,7 @@ num_simulations <- nrow(scenarios)
 cat("Total number of jobs to dispatch:", num_simulations, "\n")
 
 RNGkind("L'Ecuyer-CMRG")
-set.seed(56969)
+set.seed(269698974)
 rng_seeds <- sbc_tools$setup_lecuyer_seeds(.Random.seed, num_simulations)
 
 sim_result <- Q_rows(scenarios, sbc_tools$run_sbc_case, const=list(base_scenarios=base_scenarios, seeds=rng_seeds), export=as.list(sbc_tools), n_jobs=n_jobs, pkgs=pkg)
diff --git a/inst/sbc/sbc_tools.R b/inst/sbc/sbc_tools.R
index 453f6e6..db65eb4 100644
--- a/inst/sbc/sbc_tools.R
+++ b/inst/sbc/sbc_tools.R
@@ -2,13 +2,13 @@
 #' Utilities for SBC validation
 #'
 load_rbest_dev <- function() {
-    here::i_am("inst/sbc/sbc_tools.R")
+    if(rbest_lib_dir != .libPaths()[1]) {
+        cat("Unloading RBesT and setting libPaths to dev RBesT install.\n")
+        unloadNamespace("RBesT")
+        .libPaths(c(rbest_lib_dir, .libPaths()))
+    }
     if(!("RBesT" %in% .packages())) {
-        cat("RBesT not yet loaded, bringing up local dev version.\n")
-        library(devtools)
-        devtools::load_all(here())
-    } else {
-        cat("RBesT is already loaded.\n")
+        require(RBesT, lib.loc=rbest_lib_dir)
     }
 }
 
@@ -87,12 +87,16 @@ fit_rbest <- function(fake, draw, draw_theta, family, prior_mean_mu, prior_sd_mu
                     gaussian=cbind(y, y_se) ~ 1 | group)
 
     options(RBesT.MC.warmup=2000, RBesT.MC.iter=4000, RBesT.MC.thin=1, RBesT.MC.init=0.1,
-            RBesT.MC.control=list(adapt_delta=0.999, stepsize=0.01, max_treedepth=10,
+            RBesT.MC.control=list(##adapt_delta=0.999, ## 2024-11-21 lowered to 0.95 for the sake of performance
+                                  adapt_delta=0.95,
+                                  stepsize=0.01, max_treedepth=10,
                                   adapt_init_buffer=100, adapt_term_buffer=300))
 
     if(Ng > 5) {
         ## for the dense case we can be a bit less aggressive with the sampling tuning parameters
-        options(RBesT.MC.control=list(adapt_delta=0.95, stepsize=0.01, max_treedepth=10,
+        options(RBesT.MC.control=list(##adapt_delta=0.95,
+                                      adapt_delta=0.90,
+                                      stepsize=0.01, max_treedepth=10,
                                       adapt_init_buffer=100, adapt_term_buffer=300))
     }
     fit <- gMAP(model, data=fake, family=family,
@@ -159,58 +163,3 @@ run_sbc_case <- function(job.id, repl, data_scenario, family, mean_mu, sd_mu, sd
     c(list(job.id=job.id, time.running=runtime["elapsed"]), fit)
 }
 
-## Submits to batchtools cluster with fault tolerance, i.e.
-## resubmitting failed jobs max_num_tries times
-auto_submit <- function(jobs, registry, resources=list(), max_num_tries = 10) {
-  all_unfinished_jobs <- jobs
-
-  num_unfinished_jobs <- nrow(all_unfinished_jobs)
-  num_all_jobs <- num_unfinished_jobs
-  remaining_tries <- max_num_tries
-  all_jobs_finished <- FALSE
-  while (remaining_tries > 0 && !all_jobs_finished) {
-    remaining_tries <- remaining_tries - 1
-
-    message("Submitting jobs at ", Sys.time())
-    # Once things run fine let's submit this work to the cluster.
-    submitJobs(all_unfinished_jobs, resources=resources)
-    # Wait for results.
-    waitForJobs()
-    message("Finished waiting for jobs at ", Sys.time())
-
-    # Check status:
-    print(getStatus())
-
-    # Ensure that all jobs are done
-    if (nrow(findNotDone()) != 0) {
-      not_done_jobs <- findNotDone()
-      print(getErrorMessages(not_done_jobs))
-      ##browser()
-      ##invisible(readline(prompt="Press [enter] to continue"))
-
-      message("Some jobs did not complete. Please check the batchtools registry ", registry$file.dir)
-      all_unfinished_jobs <- inner_join(not_done_jobs, all_unfinished_jobs)
-
-      if (num_unfinished_jobs == nrow(all_unfinished_jobs) &&  nrow(all_unfinished_jobs) > 0.25 * num_all_jobs)
-      {
-        # Unfinished job count did not change -> retrying will probably not help. Abort!
-        warning("Error: unfinished job count is not decreasing. Aborting job retries.")
-        remaining_tries <- 0
-      }
-
-      if (num_unfinished_jobs == nrow(jobs))
-      {
-        # All jobs errored -> retrying will probably not help. Abort!
-        warning("Error: all jobs errored. Aborting job retries.")
-        remaining_tries <- 0
-      }
-
-      num_unfinished_jobs <- nrow(all_unfinished_jobs)
-      message("Trying to resubmit jobs. Remaining tries: ", remaining_tries, " / ", max_num_tries)
-    } else {
-      all_jobs_finished <- TRUE
-    }
-  }
-
-  invisible(all_jobs_finished)
-}
diff --git a/inst/stan/gMAP.stan b/inst/stan/gMAP.stan
index 597642f..b9bcc4d 100644
--- a/inst/stan/gMAP.stan
+++ b/inst/stan/gMAP.stan
@@ -134,10 +134,8 @@ parameters {
 }
 transformed parameters {
   vector[H] theta;
-  vector[n_groups] eta;
   vector[mX] beta;
   vector[n_tau_strata] tau;
-  vector[n_groups] tau_group;
 
   beta = beta_raw_guess[1] + beta_raw_guess[2] .* beta_raw;
 
@@ -147,14 +145,24 @@ transformed parameters {
   else
     tau = exp(tau_raw_guess[1] + tau_raw_guess[2] * tau_raw);
 
-  tau_group = tau[tau_strata_gindex];
-  
-  if (ncp)  // NCP
-    eta = xi_eta .* tau_group;
-  else // CP places overall intercept into random effect
-    eta = beta_raw_guess[1,1] + beta_raw_guess[2,1] * xi_eta;
-
-  theta = X_param * beta + eta[group_index];
+  // expand random effect to groups in loop for performance reasons
+  if(ncp) {
+    if(n_tau_strata == 1) {
+      // most common case of just one stratum which simplifies things
+      // and in ncp mode
+      for(h in 1:H) {
+        theta[h] = X_param[h] * beta + xi_eta[ group_index[h] ] * tau[1];
+      }
+    } else {
+      for(h in 1:H) {
+        theta[h] = X_param[h] * beta + xi_eta[ group_index[h] ] * tau[ tau_strata_gindex[ group_index[h] ] ];
+      }
+    }
+  } else {
+    for(h in 1:H) {
+      theta[h] = X_param[h] * beta + beta_raw_guess[1,1] + beta_raw_guess[2,1] * xi_eta[ group_index[h] ];
+    }
+  }
 }
 model {
   if (ncp) {
@@ -163,8 +171,8 @@ model {
     if(re_dist == 1) xi_eta ~ student_t(re_dist_t_df, 0, 1);
   } else {
     // random effect distribution
-    if(re_dist == 0) xi_eta ~ normal( (beta[1] - beta_raw_guess[1,1])/beta_raw_guess[2,1], tau_group / beta_raw_guess[2,1]);
-    if(re_dist == 1) xi_eta ~ student_t(re_dist_t_df, (beta[1] - beta_raw_guess[1,1])/beta_raw_guess[2,1], tau_group / beta_raw_guess[2,1]);
+    if(re_dist == 0) xi_eta ~ normal( (beta[1] - beta_raw_guess[1,1])/beta_raw_guess[2,1], tau[tau_strata_gindex] / beta_raw_guess[2,1]);
+    if(re_dist == 1) xi_eta ~ student_t(re_dist_t_df, (beta[1] - beta_raw_guess[1,1])/beta_raw_guess[2,1], tau[tau_strata_gindex] / beta_raw_guess[2,1]);
   }
 
   // assign priors to coefficients
diff --git a/man/RBesT-package.Rd b/man/RBesT-package.Rd
index e321ebc..e2802a1 100644
--- a/man/RBesT-package.Rd
+++ b/man/RBesT-package.Rd
@@ -45,6 +45,7 @@ Option \tab Default \tab Description \cr
 \tab \code{rel.tol=.Machine$double.eps^0.25,} \tab \cr
 \tab \code{abs.tol=.Machine$double.eps^0.25,} \tab \cr
 \tab \code{subdivisions=1E3)} \tab \cr
+\code{RBesT.integrate_prob_eps} \tab \code{1E-6} \tab probability mass left out from tails if integration needs to be restricted in range \cr
 }
 }
 
diff --git a/man/ess.Rd b/man/ess.Rd
index e277ae9..464388a 100644
--- a/man/ess.Rd
+++ b/man/ess.Rd
@@ -120,14 +120,14 @@ round(sum(ab_matched))
 ess(bmix, method="morita")
 
 # Predictive consistency of elir
-\donttest{
-n_forward <- 1E2
+n_forward <- 1E1
 bmixPred <- preddist(bmix, n=n_forward)
-pred_samp <- rmix(bmixPred, 1E3)
+pred_samp <- rmix(bmixPred, 1E2)
+# use more samples here for greater accuracy
+# pred_samp <- rmix(bmixPred, 1E3) 
 pred_ess <- sapply(pred_samp, function(r) ess(postmix(bmix, r=r, n=n_forward), "elir") )
 ess(bmix, "elir")
 mean(pred_ess) - n_forward
-}
 
 # Normal mixture example
 nmix <- mixnorm(rob=c(0.5, 0, 2), inf=c(0.5, 3, 4), sigma=10)
diff --git a/man/integrate_density_log.Rd b/man/integrate_density_log.Rd
index 1989106..a1682ac 100644
--- a/man/integrate_density_log.Rd
+++ b/man/integrate_density_log.Rd
@@ -14,7 +14,7 @@ integrate_density_log(
   mix,
   Lplower = -Inf,
   Lpupper = Inf,
-  eps = 1e-09
+  eps = getOption("RBesT.integrate_prob_eps", 1e-06)
 )
 }
 \arguments{
diff --git a/man/mixbeta.Rd b/man/mixbeta.Rd
index 8f2f6a5..e930df4 100644
--- a/man/mixbeta.Rd
+++ b/man/mixbeta.Rd
@@ -86,11 +86,11 @@ print(bm4)
 }
 \seealso{
 Other mixdist: 
+\code{\link{mix}},
 \code{\link{mixcombine}()},
 \code{\link{mixgamma}()},
 \code{\link{mixmvnorm}()},
 \code{\link{mixnorm}()},
-\code{\link{mixplot}},
-\code{\link{mix}}
+\code{\link{mixplot}}
 }
 \concept{mixdist}
diff --git a/man/mixcombine.Rd b/man/mixcombine.Rd
index 4925797..eb43676 100644
--- a/man/mixcombine.Rd
+++ b/man/mixcombine.Rd
@@ -35,17 +35,16 @@ bm <- mixbeta(inf=c(0.5, 10, 100), inf2=c(0.5, 30, 80))
 # robustified with mixcombine, i.e. a 10\% uninformative part added
 unif <- mixbeta(rob=c(1,1,1))
 mixcombine(bm, unif, weight=c(9, 1))
-
 }
 \seealso{
 \code{\link{robustify}}
 
 Other mixdist: 
+\code{\link{mix}},
 \code{\link{mixbeta}()},
 \code{\link{mixgamma}()},
 \code{\link{mixmvnorm}()},
 \code{\link{mixnorm}()},
-\code{\link{mixplot}},
-\code{\link{mix}}
+\code{\link{mixplot}}
 }
 \concept{mixdist}
diff --git a/man/mixgamma.Rd b/man/mixgamma.Rd
index e1d97cb..8d0950b 100644
--- a/man/mixgamma.Rd
+++ b/man/mixgamma.Rd
@@ -104,11 +104,11 @@ gfmix <- mixgamma(rob1=c(0.15, mn2gamma(2, 1)), rob2=c(0.15, ms2gamma(2, 5)), in
 }
 \seealso{
 Other mixdist: 
+\code{\link{mix}},
 \code{\link{mixbeta}()},
 \code{\link{mixcombine}()},
 \code{\link{mixmvnorm}()},
 \code{\link{mixnorm}()},
-\code{\link{mixplot}},
-\code{\link{mix}}
+\code{\link{mixplot}}
 }
 \concept{mixdist}
diff --git a/man/mixmvnorm.Rd b/man/mixmvnorm.Rd
index 3c352d1..d97dc2a 100644
--- a/man/mixmvnorm.Rd
+++ b/man/mixmvnorm.Rd
@@ -81,11 +81,11 @@ colMeans(mixSamp1)
 }
 \seealso{
 Other mixdist: 
+\code{\link{mix}},
 \code{\link{mixbeta}()},
 \code{\link{mixcombine}()},
 \code{\link{mixgamma}()},
 \code{\link{mixnorm}()},
-\code{\link{mixplot}},
-\code{\link{mix}}
+\code{\link{mixplot}}
 }
 \concept{mixdist}
diff --git a/man/mixnorm.Rd b/man/mixnorm.Rd
index cac7eb3..e28033e 100644
--- a/man/mixnorm.Rd
+++ b/man/mixnorm.Rd
@@ -107,11 +107,11 @@ ess(nm)
 }
 \seealso{
 Other mixdist: 
+\code{\link{mix}},
 \code{\link{mixbeta}()},
 \code{\link{mixcombine}()},
 \code{\link{mixgamma}()},
 \code{\link{mixmvnorm}()},
-\code{\link{mixplot}},
-\code{\link{mix}}
+\code{\link{mixplot}}
 }
 \concept{mixdist}
diff --git a/man/mixplot.Rd b/man/mixplot.Rd
index 7fe703c..8bd5706 100644
--- a/man/mixplot.Rd
+++ b/man/mixplot.Rd
@@ -25,7 +25,7 @@ density in addition to the density.}
 
 \item{size}{controls the linesize in plots.}
 
-\item{...}{extra arguments passed on to the \code{\link[ggplot2]{qplot}} call.}
+\item{...}{extra arguments passed on to the plotted function.}
 }
 \value{
 A \code{\link[ggplot2]{ggplot}} object is returned.
@@ -92,11 +92,11 @@ pl + ggtitle("Normal 2-Component Mixture")
 }
 \seealso{
 Other mixdist: 
+\code{\link{mix}},
 \code{\link{mixbeta}()},
 \code{\link{mixcombine}()},
 \code{\link{mixgamma}()},
 \code{\link{mixmvnorm}()},
-\code{\link{mixnorm}()},
-\code{\link{mix}}
+\code{\link{mixnorm}()}
 }
 \concept{mixdist}
diff --git a/man/oc1S.Rd b/man/oc1S.Rd
index c1797da..9884172 100644
--- a/man/oc1S.Rd
+++ b/man/oc1S.Rd
@@ -143,8 +143,8 @@ designC2_nL(theta_eval)
 }
 \seealso{
 Other design1S: 
-\code{\link{decision1S_boundary}()},
 \code{\link{decision1S}()},
+\code{\link{decision1S_boundary}()},
 \code{\link{pos1S}()}
 }
 \concept{design1S}
diff --git a/man/oc2S.Rd b/man/oc2S.Rd
index a817d1f..7ff428b 100644
--- a/man/oc2S.Rd
+++ b/man/oc2S.Rd
@@ -152,8 +152,8 @@ Robust meta-analytic-predictive priors in clinical trials with historical contro
 }
 \seealso{
 Other design2S: 
-\code{\link{decision2S_boundary}()},
 \code{\link{decision2S}()},
+\code{\link{decision2S_boundary}()},
 \code{\link{pos2S}()}
 }
 \concept{design2S}
diff --git a/man/plot.gMAP.Rd b/man/plot.gMAP.Rd
index 538f1c4..ffa83f0 100644
--- a/man/plot.gMAP.Rd
+++ b/man/plot.gMAP.Rd
@@ -14,7 +14,7 @@
 \item{...}{Ignored.}
 }
 \value{
-The function returns a list of \code{\link{ggplot}}
+The function returns a list of \code{\link[ggplot2:ggplot]{ggplot}}
 objects.
 }
 \description{
diff --git a/man/pos1S.Rd b/man/pos1S.Rd
index fb18882..63d9aa9 100644
--- a/man/pos1S.Rd
+++ b/man/pos1S.Rd
@@ -121,8 +121,8 @@ pos_ia(post_ia)
 }
 \seealso{
 Other design1S: 
-\code{\link{decision1S_boundary}()},
 \code{\link{decision1S}()},
+\code{\link{decision1S_boundary}()},
 \code{\link{oc1S}()}
 }
 \concept{design1S}
diff --git a/man/pos2S.Rd b/man/pos2S.Rd
index 4e3f026..4486039 100644
--- a/man/pos2S.Rd
+++ b/man/pos2S.Rd
@@ -155,8 +155,8 @@ pos_final(postP_interim, postT_interim)
 }
 \seealso{
 Other design2S: 
-\code{\link{decision2S_boundary}()},
 \code{\link{decision2S}()},
+\code{\link{decision2S_boundary}()},
 \code{\link{oc2S}()}
 }
 \concept{design2S}
diff --git a/man/predict.gMAP.Rd b/man/predict.gMAP.Rd
index 2c46d1a..c87e32d 100644
--- a/man/predict.gMAP.Rd
+++ b/man/predict.gMAP.Rd
@@ -89,7 +89,6 @@ summary(pred_cov)
 # obtain a prediction for new data with specific covariates
 pred_new <- predict(map, data.frame(country="CH", study=12))
 pred_new
-
 ## Recover user set sampling defaults
 options(.user_mc_options)
 
diff --git a/tests/testthat/helper-utils.R b/tests/testthat/helper-utils.R
new file mode 100644
index 0000000..95db7f2
--- /dev/null
+++ b/tests/testthat/helper-utils.R
@@ -0,0 +1,10 @@
+
+source_example <- function(example, env=parent.frame(), disable_plots=TRUE) {
+    ex_source <- readLines(system.file("examples", example, package="RBesT", mustWork=TRUE))
+    if(disable_plots) {
+        ex_source <- grep("plot\\(", ex_source, value=TRUE, invert=TRUE)
+    }
+    suppressMessages(ex <- source(textConnection(ex_source),
+                                  local=env, echo=FALSE, verbose=FALSE))
+    invisible(ex)
+}
diff --git a/tests/testthat/test-EM.R b/tests/testthat/test-EM.R
index 3fd0705..c1eae92 100644
--- a/tests/testthat/test-EM.R
+++ b/tests/testthat/test-EM.R
@@ -1,4 +1,4 @@
-context("EM: Expectation-Maximization")
+#context("EM: Expectation-Maximization")
 
 ## test that the EM algorithms recover reliably test distributions;
 ## test criterium is a "sufficiently" small KL divergence
@@ -145,23 +145,23 @@ EM_mvn_test <- function(mixTest, seed, Nsim=1e4, verbose=FALSE, ...) {
     expect_true(all(EMmix1 == EMmix2), info="Result of EM is independent of random seed.")
 }
 
-test_that("Normal EM fits single component",     EM_test(ref$norm_single, 3453563, Nsim, verbose))
-test_that("Normal EM fits heavy-tailed mixture", EM_test(ref$norm_heavy,  9275624, Nsim, verbose))
-test_that("Normal EM fits bi-modal mixture",     EM_test(ref$norm_bi,     9345726, Nsim, verbose))
+test_that("Normal EM fits single component",     { EM_test(ref$norm_single, 3453563, Nsim, verbose) })
+test_that("Normal EM fits heavy-tailed mixture", { EM_test(ref$norm_heavy,  9275624, Nsim, verbose) })
+test_that("Normal EM fits bi-modal mixture",     { EM_test(ref$norm_bi,     9345726, Nsim, verbose) })
 
-test_that("Multivariate Normal EM fits single component",     EM_mvn_test(ref$mvnorm_single, 3453563, Nsim, verbose))
-test_that("Multivariate Normal EM fits heavy-tailed mixture", EM_mvn_test(ref$mvnorm_heavy,  9275624, Nsim, verbose))
-test_that("Multivariate Normal EM fits bi-modal mixture",     EM_mvn_test(ref$mvnorm_bi,     9345726, Nsim, verbose))
-test_that("Multivariate Normal EM fits bi-modal mixture 1D",     EM_mvn_test(ref$mvnorm_bi_1D,     9345726, Nsim, verbose))
+test_that("Multivariate Normal EM fits single component",     { EM_mvn_test(ref$mvnorm_single, 3453563, max(1E4, Nsim), verbose) })
+test_that("Multivariate Normal EM fits heavy-tailed mixture", { EM_mvn_test(ref$mvnorm_heavy,  9275624, max(1E4, Nsim), verbose) })
+test_that("Multivariate Normal EM fits bi-modal mixture",     { EM_mvn_test(ref$mvnorm_bi,     9345726, max(1E4, Nsim), verbose) })
+test_that("Multivariate Normal EM fits bi-modal mixture 1D",     { EM_mvn_test(ref$mvnorm_bi_1D,     9345726, max(1E4, Nsim), verbose) })
 
-test_that("Gamma EM fits single component",      EM_test(ref$gamma_single, 9345835, Nsim, verbose))
-test_that("Gamma EM fits heavy-tailed mixture",  EM_test(ref$gamma_heavy,  5629389, Nsim, verbose))
-test_that("Gamma EM fits bi-modal mixture",      EM_test(ref$gamma_bi,     9373515, Nsim, verbose))
+test_that("Gamma EM fits single component",      { EM_test(ref$gamma_single, 9345835, Nsim, verbose) })
+test_that("Gamma EM fits heavy-tailed mixture",  { EM_test(ref$gamma_heavy,  5629389, Nsim, verbose) })
+test_that("Gamma EM fits bi-modal mixture",      { EM_test(ref$gamma_bi,     9373515, Nsim, verbose) })
 
-test_that("Beta EM fits single component",       EM_test(ref$beta_single, 7265355, Nsim, verbose))
-test_that("Beta EM fits single component with mass at boundary", EM_test(ref$beta_single_alt, 7265355, Nsim, verbose, constrain_gt1=FALSE))
-test_that("Beta EM fits heavy-tailed mixture",   EM_test(ref$beta_heavy,  2946562, Nsim, verbose))
-test_that("Beta EM fits bi-modal mixture",       EM_test(ref$beta_bi,     9460370, Nsim, verbose))
+test_that("Beta EM fits single component",       { EM_test(ref$beta_single, 7265355, Nsim, verbose) })
+test_that("Beta EM fits single component with mass at boundary", { EM_test(ref$beta_single_alt, 7265355, Nsim, verbose, constrain_gt1=FALSE) })
+test_that("Beta EM fits heavy-tailed mixture",   { EM_test(ref$beta_heavy,  2946562, Nsim, verbose) })
+test_that("Beta EM fits bi-modal mixture",       { EM_test(ref$beta_bi,     9460370, Nsim, verbose) })
 
 test_that("Constrained Beta EM respects a>1 & b>1", {
     unconstrained  <- mixbeta(c(0.6, 2.8, 64), c(0.25, 0.5, 0.92), c(0.15, 3, 15))
diff --git a/tests/testthat/test-gMAP.R b/tests/testthat/test-gMAP.R
index 512ee92..be4b222 100644
--- a/tests/testthat/test-gMAP.R
+++ b/tests/testthat/test-gMAP.R
@@ -1,4 +1,3 @@
-context("gMAP: Generalized Meta Analytic Predictive")
 
 ## test gMAP results using SBC and with matching rstanarm models
 
@@ -106,14 +105,14 @@ rate <- round(-log(0.05)/2, 1)
 test_that("gMAP matches RStanArm binomial family", {
     skip("RStanArm has issues loading since 2024-01-02 on CI/CD systems.")
               skip_on_cran()
-              best_run <- gMAP(cbind(r, n-r) ~ 1 | study,
+              suppressWarnings( best_run <- gMAP(cbind(r, n-r) ~ 1 | study,
                                data=AS,
                                family=binomial,
                                tau.dist="Exp",
                                tau.prior=c(rate),
                                beta.prior=cbind(0, 2)
-                               )
-              out <- capture.output(rstanarm_run <- make_rstanarm_ref(
+                               ) )
+              suppressWarnings( out <- capture.output(rstanarm_run <- make_rstanarm_ref(
                   stan_glmer(cbind(r, n-r) ~ 1 + (1|study),
                              data=AS,
                              family=binomial,
@@ -127,6 +126,7 @@ test_that("gMAP matches RStanArm binomial family", {
                              prior_intercept=normal(0,2,autoscale=FALSE),
                              prior_covariance=decov(1, 1, 1, 1/rate)
                              )))
+                  )
               cmp_reference(best_gmap=best_run, OB_ref=rstanarm_run)
           })
 
@@ -148,9 +148,10 @@ test_that("gMAP processes single trial case", {
           )
 
 test_that("gMAP processes not continuously labeled studies", {
-              out <- capture.output(map1 <- gMAP(cbind(r, n-r) ~ 1 | study, data=AS[-1,],
+              suppressWarnings( out <- capture.output(map1 <- gMAP(cbind(r, n-r) ~ 1 | study, data=AS[-1,],
                                                  family=binomial, tau.dist="HalfNormal", tau.prior=0.5,
                                                  iter=100, warmup=50, chains=1, thin=1))
+                               )
               expect_true(nrow(fitted(map1)) == nrow(AS) - 1)
           })
 
@@ -163,7 +164,7 @@ test_that("gMAP reports divergences", {
                                                                  tau.dist="Uniform", tau.prior=cbind(0, 1000),
                                                                  beta.prior=cbind(0,1E5),
                                                                  iter=1000, warmup=0, chains=1, thin=1, init=10)))
-              sp <- rstan::get_sampler_params(mcmc_div$fit)[[1]]
+              sp <- rstan::get_sampler_params(mcmc_div$fit, inc_warmup=FALSE)[[1]]
               expect_true(sum(sp[,"divergent__"]) > 0)
           })
 
@@ -173,31 +174,31 @@ do.call(options, std_sampling)
 test_that("gMAP handles extreme response rates", {
               n <- 5
               data1 <- data.frame(n=c(n,n,n,n),r=c(5,5,5,5), study=1)
-              map1 <- gMAP(cbind(r, n-r) ~ 1 | study, family=binomial,
-                           data=data1, tau.dist="HalfNormal",
-                           tau.prior=2.0, beta.prior=2,
-                           warmup=100, iter=200, chains=1, thin=1)
+              suppressWarnings( map1 <- gMAP(cbind(r, n-r) ~ 1 | study, family=binomial,
+                                             data=data1, tau.dist="HalfNormal",
+                                             tau.prior=2.0, beta.prior=2,
+                                             warmup=100, iter=200, chains=1, thin=1) )
               expect_true(nrow(fitted(map1)) == 4)
               data2 <- data.frame(n=c(n,n,n,n),r=c(0,0,0,0), study=1)
-              map2 <- gMAP(cbind(r, n-r) ~ 1 | study, family=binomial,
-                           data=data2, tau.dist="HalfNormal",
-                           tau.prior=2.0, beta.prior=2,
-                           warmup=100, iter=200, chains=1, thin=1)
+              suppressWarnings( map2 <- gMAP(cbind(r, n-r) ~ 1 | study, family=binomial,
+                                             data=data2, tau.dist="HalfNormal",
+                                             tau.prior=2.0, beta.prior=2,
+                                             warmup=100, iter=200, chains=1, thin=1) )
               expect_true(nrow(fitted(map2)) == 4)
               data3 <- data.frame(n=c(n,n,n,n),r=c(5,5,5,5), study=c(1,1,2,2))
-              map3 <- gMAP(cbind(r, n-r) ~ 1 | study, family=binomial,
-                           data=data3, tau.dist="HalfNormal",
-                           tau.prior=2.0, beta.prior=2,
-                           warmup=100, iter=200, chains=1, thin=1)
+              suppressWarnings( map3 <- gMAP(cbind(r, n-r) ~ 1 | study, family=binomial,
+                                             data=data3, tau.dist="HalfNormal",
+                                             tau.prior=2.0, beta.prior=2,
+                                             warmup=100, iter=200, chains=1, thin=1) )
               expect_true(nrow(fitted(map3)) == 4)
           })
 
 test_that("gMAP handles fixed tau case", {
-              map1 <- gMAP(cbind(r, n-r) ~ 1 | study, family=binomial,
-                           data=AS, tau.dist="Fixed",
-                           tau.prior=0.5, beta.prior=2,
-                           warmup=100, iter=200, chains=1, thin=1)
-              expect_true(map1$Rhat.max >= 1)
+    suppressWarnings( map1 <- gMAP(cbind(r, n-r) ~ 1 | study, family=binomial,
+                                   data=AS, tau.dist="Fixed",
+                                   tau.prior=0.5, beta.prior=2,
+                                   warmup=100, iter=200, chains=1, thin=1) )
+    expect_true(map1$Rhat.max >= 1)
           })
 
 test_that("gMAP labels data rows correctly when using covariates", {
diff --git a/tests/testthat/test-mixdiff.R b/tests/testthat/test-mixdiff.R
index fbe7fae..154d05c 100644
--- a/tests/testthat/test-mixdiff.R
+++ b/tests/testthat/test-mixdiff.R
@@ -1,4 +1,3 @@
-context("mixdiff: Mixture Difference Distribution")
 
 ## test that calculations for the cumulative distribution function of
 ## differences in mixture distributions are correct by comparison with
@@ -64,18 +63,18 @@ mixdiff_cmp_norm <- function(case, mixref) {
     expect_true(all(res_quants < eps))
 }
 
-test_that("Difference in beta variates evaluates correctly", mixdiff_cmp(beta))
-test_that("Difference in beta mixture variates evaluates correctly", mixdiff_cmp(betaMix))
+test_that("Difference in beta variates evaluates correctly", { mixdiff_cmp(beta) })
+test_that("Difference in beta mixture variates evaluates correctly", { mixdiff_cmp(betaMix) })
 
-test_that("Difference in gamma variates evaluates correctly", mixdiff_cmp(gamma))
-test_that("Difference in gamma mixture variates evaluates correctly", mixdiff_cmp(gammaMix))
+test_that("Difference in gamma variates evaluates correctly", { mixdiff_cmp(gamma) })
+test_that("Difference in gamma mixture variates evaluates correctly", { mixdiff_cmp(gammaMix) })
 
-test_that("Difference in normal variates evaluates correctly", mixdiff_cmp(norm))
-test_that("Difference in normal mixture variates evaluates correctly", mixdiff_cmp(normMix))
+test_that("Difference in normal variates evaluates correctly", { mixdiff_cmp(norm) })
+test_that("Difference in normal mixture variates evaluates correctly", { mixdiff_cmp(normMix) })
 
 ## for the normal we can use exact analytical results
-test_that("Difference in normal variates evaluates (analytically) correctly", mixdiff_cmp_norm(norm,norm_ref))
-test_that("Difference in normal mixture variates evaluates (analytically) correctly", mixdiff_cmp_norm(normMix,normMix_ref))
+test_that("Difference in normal variates evaluates (analytically) correctly", { mixdiff_cmp_norm(norm,norm_ref) })
+test_that("Difference in normal mixture variates evaluates (analytically) correctly", { mixdiff_cmp_norm(normMix,normMix_ref) })
 
 ## now test difference distributions on the link-transformed scales
 ## (the cannonical cases)
@@ -86,14 +85,14 @@ apply_link <- function(dists, link) {
 }
  
 ## log-odds
-test_that("Difference in beta variates with logit link evaluates correctly", mixdiff_cmp(apply_link(beta, RBesT:::logit_dlink)))
-test_that("Difference in beta mixture variates with logit link evaluates correctly", mixdiff_cmp(apply_link(betaMix, RBesT:::logit_dlink)))
+test_that("Difference in beta variates with logit link evaluates correctly", { mixdiff_cmp(apply_link(beta, RBesT:::logit_dlink)) })
+test_that("Difference in beta mixture variates with logit link evaluates correctly", { mixdiff_cmp(apply_link(betaMix, RBesT:::logit_dlink)) })
 
 ## relative risk
-test_that("Difference in beta variates with log link evaluates correctly", mixdiff_cmp(apply_link(beta, RBesT:::log_dlink)))
-test_that("Difference in beta mixture variates with log link evaluates correctly", mixdiff_cmp(apply_link(betaMix, RBesT:::log_dlink)))
+test_that("Difference in beta variates with log link evaluates correctly", { mixdiff_cmp(apply_link(beta, RBesT:::log_dlink)) })
+test_that("Difference in beta mixture variates with log link evaluates correctly", { mixdiff_cmp(apply_link(betaMix, RBesT:::log_dlink)) })
 
 ## relative counts
-test_that("Difference in gamma variates with log link evaluates correctly", mixdiff_cmp(apply_link(gamma, RBesT:::log_dlink)))
-test_that("Difference in gamma mixture variates with log link evaluates correctly", mixdiff_cmp(apply_link(gammaMix, RBesT:::log_dlink)))
+test_that("Difference in gamma variates with log link evaluates correctly", { mixdiff_cmp(apply_link(gamma, RBesT:::log_dlink)) })
+test_that("Difference in gamma mixture variates with log link evaluates correctly", { mixdiff_cmp(apply_link(gammaMix, RBesT:::log_dlink)) })
 
diff --git a/tests/testthat/test-mixdist.R b/tests/testthat/test-mixdist.R
index 202e618..84a9cdf 100644
--- a/tests/testthat/test-mixdist.R
+++ b/tests/testthat/test-mixdist.R
@@ -1,4 +1,3 @@
-context("mixdist: Mixture Distribution")
 
 ## various tests around mixture distributions
 set.seed(234534)
@@ -41,14 +40,14 @@ pmix_lower_tail_test <- function(mix, N=Nsamp_quant) {
     do_test(preddist(mix, n=100))
 }
 
-test_that("Cumulative beta distribution function evaluates lower.tail correctly", pmix_lower_tail_test(beta))
-test_that("Cumulative beta mixture distribution function evaluates lower.tail correctly", pmix_lower_tail_test(betaMix))
+test_that("Cumulative beta distribution function evaluates lower.tail correctly", { pmix_lower_tail_test(beta) })
+test_that("Cumulative beta mixture distribution function evaluates lower.tail correctly", { pmix_lower_tail_test(betaMix) })
 
-test_that("Cumulative normal distribution function evaluates lower.tail correctly", pmix_lower_tail_test(norm))
-test_that("Cumulative normal mixture distribution function evaluates lower.tail correctly", pmix_lower_tail_test(normMix))
+test_that("Cumulative normal distribution function evaluates lower.tail correctly", { pmix_lower_tail_test(norm) })
+test_that("Cumulative normal mixture distribution function evaluates lower.tail correctly", { pmix_lower_tail_test(normMix) })
 
-test_that("Cumulative gamma distribution function evaluates lower.tail correctly", pmix_lower_tail_test(gamma))
-test_that("Cumulative gamma mixture distribution function evaluates lower.tail correctly", pmix_lower_tail_test(gammaMix))
+test_that("Cumulative gamma distribution function evaluates lower.tail correctly", { pmix_lower_tail_test(gamma) })
+test_that("Cumulative gamma mixture distribution function evaluates lower.tail correctly", { pmix_lower_tail_test(gammaMix) })
 
 ## tests the quantile and distribution function against simulated samples
 mix_simul_test <- function(mix, eps, qtest, ptest = seq(0.1, 0.9, by=0.1), S=Nsamp_equant) {
@@ -63,40 +62,52 @@ mix_simul_test <- function(mix, eps, qtest, ptest = seq(0.1, 0.9, by=0.1), S=Nsa
     expect_true(all(res_probs < eps))
 }
 
-test_that("Beta quantile function is correct", mix_simul_test(beta, eps, c(0.1, 0.9)))
-test_that("Beta mixture quantile function is correct", mix_simul_test(betaMix, eps, c(0.1, 0.9)))
+test_that("Beta quantile function is correct", { mix_simul_test(beta, eps, c(0.1, 0.9)) })
+test_that("Beta mixture quantile function is correct", { mix_simul_test(betaMix, eps, c(0.1, 0.9)) })
 
-test_that("Normal quantile function is correct", mix_simul_test(norm, eps, c(-1, 0)))
-test_that("Normal mixture quantile function is correct", mix_simul_test(normMix, eps, c(4, 1)))
-test_that("Normal mixture with very weak component quantile function is correct", mix_simul_test(normMixWeak, eps, c(4, 1)))
+test_that("Normal quantile function is correct", { mix_simul_test(norm, eps, c(-1, 0)) })
+test_that("Normal mixture quantile function is correct", { mix_simul_test(normMix, eps, c(4, 1)) })
+test_that("Normal mixture with very weak component quantile function is correct", { mix_simul_test(normMixWeak, eps, c(4, 1)) })
 
-test_that("Gamma quantile function is correct", mix_simul_test(gamma, eps, c(2, 7)))
-test_that("Gamma mixture quantile function is correct", mix_simul_test(gammaMix, eps, c(2, 7), ptest = seq(0.2, 0.8, by=0.1)))
+test_that("Gamma quantile function is correct", { mix_simul_test(gamma, eps, c(2, 7)) })
+test_that("Gamma mixture quantile function is correct", { mix_simul_test(gammaMix, eps, c(2, 7), ptest = seq(0.2, 0.8, by=0.1)) })
 
 ## problematic gamma (triggers internally a fallback to root finding)
 gammaMix2 <- mixgamma(c(8.949227e-01, 7.051570e-01, 6.125121e-02),
                       c(1.049106e-01, 3.009986e-01, 5.169626e-04),
                       c(1.666667e-04, 1.836051e+04, 1.044005e-02))
 
-test_that("Singular gamma mixture quantile function is correct", mix_simul_test(gammaMix2, 10*eps, c(1, 1E3), ptest = seq(0.2, 0.8, by=0.1)))
+test_that("Singular gamma mixture quantile function is correct", { mix_simul_test(gammaMix2, 10*eps, c(1, 1E3), ptest = seq(0.2, 0.8, by=0.1)) })
 
 
 consistent_cdf <- function(mix, values) {
     dens <- dmix(mix, values)
     cdf <- pmix(mix, values)
+    lcdf <- pmix(mix, values, log.p=TRUE)
     expect_true(all(diff(cdf) >= 0))
-    expect_numeric(dens, any.missing=FALSE)
-    expect_numeric(cdf, any.missing=FALSE)
+    expect_numeric(dens, lower=0, finite=TRUE, any.missing=FALSE)
+    expect_numeric(cdf, lower=0, upper=1, finite=TRUE, any.missing=FALSE)
+    expect_numeric(lcdf, upper=0, finite=FALSE, any.missing=FALSE)
 }
 
-test_that("Beta CDF function is consistent", consistent_cdf(beta, seq(0.1, 0.9, by=0.1)))
-test_that("Beta mixture CDF function is consistent", consistent_cdf(betaMix, seq(0.1, 0.9, by=0.1)))
+consistent_ccdf <- function(mix, values) {
+    dens <- dmix(mix, values)
+    ccdf <- pmix(mix, values, FALSE)
+    lccdf <- pmix(mix, values, FALSE, TRUE)
+    expect_true(all(diff(ccdf) <= 0))
+    expect_numeric(dens, lower=0, finite=TRUE, any.missing=FALSE)
+    expect_numeric(ccdf, lower=0, upper=1, finite=TRUE, any.missing=FALSE)
+    expect_numeric(lccdf, upper=0, finite=FALSE, any.missing=FALSE)
+}
+
+test_that("Beta CDF function is consistent", { consistent_cdf(beta, seq(0.1, 0.9, by=0.1)) })
+test_that("Beta mixture CDF function is consistent", { consistent_cdf(betaMix, seq(0.1, 0.9, by=0.1)) })
 
-test_that("Normal CDF is consistent", consistent_cdf(norm, seq(-2, 2, by=0.1)))
-test_that("Normal mixture CDF is consistent", consistent_cdf(norm, seq(-2, 2, by=0.1)))
+test_that("Normal CDF is consistent", { consistent_cdf(norm, seq(-2, 2, by=0.1)) })
+test_that("Normal mixture CDF is consistent", { consistent_cdf(norm, seq(-2, 2, by=0.1)) })
 
-test_that("Gamma CDF function is consistent", consistent_cdf(gamma, seq(2, 7, by=0.1)))
-test_that("Gamma mixture CDF function is consistent", consistent_cdf(gammaMix, seq(2, 7, by=0.1)))
+test_that("Gamma CDF function is consistent", { consistent_cdf(gamma, seq(2, 7, by=0.1)) })
+test_that("Gamma mixture CDF function is consistent", { consistent_cdf(gammaMix, seq(2, 7, by=0.1)) })
 
 
 ## problematic beta which triggers that the cumulative of the
@@ -104,8 +115,12 @@ test_that("Gamma mixture CDF function is consistent", consistent_cdf(gammaMix, s
 ## again once 2.18 is out)
 
 ## problematic Beta density
-bm <- mixbeta(c(1.0, 298.30333970, 146.75306521))
-test_that("Problematic (1) BetaBinomial CDF function is consistent", consistent_cdf(preddist(bm, n=50), 0:50))
+bm1 <- mixbeta(c(1.0, 298.30333970, 146.75306521))
+test_that("Problematic (1) BetaBinomial CDF function is consistent", { consistent_cdf(preddist(bm1, n=50), 0:50) })
+test_that("Problematic (1) BetaBinomial CCDF function is consistent", { consistent_ccdf(preddist(bm1, n=50), 0:50) })
+bm2 <- mixbeta(c(1.0, 3 + 1/3, 47 + 1/3))
+test_that("Problematic (2) BetaBinomial CDF function is consistent", { consistent_cdf(preddist(bm2, n=50), 0:50) })
+test_that("Problematic (2) BetaBinomial CCDF function is consistent", { consistent_ccdf(preddist(bm2, n=50), 0:50) })
 
 ## tests for the multivariate normal mixture density
 p <- 4
diff --git a/tests/testthat/test-mixstanvar.R b/tests/testthat/test-mixstanvar.R
index b5b463c..813e814 100644
--- a/tests/testthat/test-mixstanvar.R
+++ b/tests/testthat/test-mixstanvar.R
@@ -1,4 +1,4 @@
-context("mixstanvar: Mixture Distribution brms Adapter")
+
 
 skip_if_not_installed("brms")
 
@@ -56,9 +56,17 @@ mixstanvar_simul_test <- function(mix,
 mixstanvar_test <- function(mix,
                             brms_args) {
     brms_prior_empty <- do.call(brms::brm, c(brms_args, list(seed=1423545, refresh=0, sample_prior="only", stanvars=mixstanvar(prior=mix), empty=TRUE)))
-    stan_dist <- paste0("mix", gsub("Mix$", "", class(mix)[1]), "_lpdf")
+    mix_class <- gsub("Mix$", "", class(mix)[1])
+    stan_dist_lpdf <- paste0("mix", mix_class, "_lpdf")
+    stan_dist_lcdf <- paste0("mix", mix_class, "_lcdf")
+    stan_dist_lccdf <- paste0("mix", mix_class, "_lccdf")
+    stan_dist_cdf <- paste0("mix", mix_class, "_cdf")
     ## look for the declared density in Stan
-    expect_true(grep(stan_dist, brms::stancode(brms_prior_empty)) == 1, info="Looking for declared Stan mixture density in generated brms Stan code.")
+    stan_code <- brms::stancode(brms_prior_empty)
+    expect_true(grep(stan_dist_lpdf, stan_code) == 1, info="Looking for declared Stan mixture density pdf in generated brms Stan code.")
+    expect_true(grep(stan_dist_lcdf, stan_code) == 1, info="Looking for declared Stan mixture density cdf in generated brms Stan code.")
+    expect_true(grep(stan_dist_lccdf, stan_code) == 1, info="Looking for declared Stan mixture density ccdf in generated brms Stan code.")
+    expect_true(grep(stan_dist_cdf, stan_code) == 1, info="Looking for declared Stan mixture density natural scale cdf in generated brms Stan code.")
     ## now check for the mixture being passed to Stan as data
     stan_data <- brms::standata(brms_prior_empty)
     for(i in 1:3) {    
@@ -71,27 +79,46 @@ brms_beta_args <- list(formula=brms::bf(r | trials(n) ~ 1, family=brms::brmsfami
                        data=data.frame(r=0, n=0),
                        prior=brms::prior(mixbeta(prior_w, prior_a, prior_b), coef=Intercept))
 
-test_that("Beta quantiles are correct for brms sampled prior", mixstanvar_simul_test(beta, brms_beta_args, eps, c(0.1, 0.9)))
-test_that("Beta prior is declared correctly in brms generated model and data", mixstanvar_test(beta, brms_beta_args))
-test_that("Beta mixture quantiles are correct for brms sampled prior", mixstanvar_simul_test(betaMix, brms_beta_args, eps, c(0.1, 0.9)))
-test_that("Beta mixture prior is declared correctly in brms generated model and data", mixstanvar_test(betaMix, brms_beta_args))
+test_that("Beta quantiles are correct for brms sampled prior", { mixstanvar_simul_test(beta, brms_beta_args, eps, c(0.1, 0.9)) })
+test_that("Beta prior is declared correctly in brms generated model and data", { mixstanvar_test(beta, brms_beta_args) })
+test_that("Beta mixture quantiles are correct for brms sampled prior", { mixstanvar_simul_test(betaMix, brms_beta_args, eps, c(0.1, 0.9)) })
+test_that("Beta mixture prior is declared correctly in brms generated model and data", { mixstanvar_test(betaMix, brms_beta_args) })
+
+brms_beta_trunc_args <- list(formula=brms::bf(r | trials(n) ~ 1, family=brms::brmsfamily("binomial", link="identity"), center=FALSE),
+                             data=data.frame(r=0, n=0),
+                             prior=brms::prior(mixbeta(prior_w, prior_a, prior_b), class=b, lb=0.1, ub=0.9))
+test_that("Beta (truncated) prior is declared correctly in brms generated model and data", { mixstanvar_test(beta, brms_beta_trunc_args) })
+test_that("Beta mixture (truncated) prior is declared correctly in brms generated model and data", { mixstanvar_test(betaMix, brms_beta_trunc_args) })
 
 brms_normal_args <- list(formula=brms::bf(y ~ 1, family=brms::brmsfamily("gaussian", link="identity"), center=FALSE),
                          data=data.frame(y=0),
                          prior=brms::prior(mixnorm(prior_w, prior_m, prior_s), coef=Intercept) + brms::prior(constant(1), class=sigma))
-test_that("Normal quantiles are correct for brms sampled prior", mixstanvar_simul_test(norm, brms_normal_args, eps, c(-1, 0)))
-test_that("Normal prior is declared correctly in brms generated model and data", mixstanvar_test(norm, brms_normal_args))
-test_that("Normal mixture quantiles are correct for brms sampled prior", mixstanvar_simul_test(normMix, brms_normal_args, eps, c(2, 1), ptest=c(0.3, 0.5, 0.7)))
-test_that("Normal mixture prior is declared correctly in brms generated model and data", mixstanvar_test(normMix, brms_normal_args))
+test_that("Normal quantiles are correct for brms sampled prior", { mixstanvar_simul_test(norm, brms_normal_args, eps, c(-1, 0)) })
+test_that("Normal prior is declared correctly in brms generated model and data", { mixstanvar_test(norm, brms_normal_args) })
+test_that("Normal mixture quantiles are correct for brms sampled prior", { mixstanvar_simul_test(normMix, brms_normal_args, eps, c(2, 1), ptest=c(0.3, 0.5, 0.7)) })
+test_that("Normal mixture prior is declared correctly in brms generated model and data", { mixstanvar_test(normMix, brms_normal_args) })
+
+brms_normal_trunc_args <- list(formula=brms::bf(y ~ 1, family=brms::brmsfamily("gaussian", link="identity"), center=FALSE),
+                         data=data.frame(y=0),
+                         prior=brms::prior(mixnorm(prior_w, prior_m, prior_s), class=b, lb=-5, ub=5) + brms::prior(constant(1), class=sigma))
+test_that("Normal (truncated) prior is declared correctly in brms generated model and data", { mixstanvar_test(norm, brms_normal_trunc_args) })
+test_that("Normal mixture (truncated) prior is declared correctly in brms generated model and data", { mixstanvar_test(normMix, brms_normal_trunc_args) })
 
 brms_gamma_args <- list(formula=brms::bf(y ~ 1, family=brms::brmsfamily("gaussian", link="identity"), center=FALSE),
                         data=data.frame(y=1),
                         prior=brms::prior(mixgamma(prior_w, prior_a, prior_b), coef=Intercept) + brms::prior(constant(1), class=sigma))
 
-test_that("Gamma quantiles are correct for brms sampled prior", mixstanvar_simul_test(gamma, brms_gamma_args, eps, c(2, 7)))
-test_that("Gamma prior is declared correctly in brms generated model and data", mixstanvar_test(gamma, brms_gamma_args))
-test_that("Gamma mixture quantile function is correct for brms sampled prior", mixstanvar_simul_test(gammaMix, brms_gamma_args, eps, c(2, 7), ptest = seq(0.2, 0.8, by=0.2)))
-test_that("Gamma mixture prior is declared correctly in brms generated model and data", mixstanvar_test(gammaMix, brms_gamma_args))
+test_that("Gamma quantiles are correct for brms sampled prior", { mixstanvar_simul_test(gamma, brms_gamma_args, eps, c(2, 7)) })
+test_that("Gamma prior is declared correctly in brms generated model and data", { mixstanvar_test(gamma, brms_gamma_args) })
+test_that("Gamma mixture quantile function is correct for brms sampled prior", { mixstanvar_simul_test(gammaMix, brms_gamma_args, eps, c(2, 7), ptest = seq(0.2, 0.8, by=0.2)) })
+test_that("Gamma mixture prior is declared correctly in brms generated model and data", { mixstanvar_test(gammaMix, brms_gamma_args) })
+
+brms_gamma_trunc_args <- list(formula=brms::bf(y ~ 1, family=brms::brmsfamily("gaussian", link="identity"), center=FALSE),
+                              data=data.frame(y=1),
+                              prior=brms::prior(mixgamma(prior_w, prior_a, prior_b), class=b, lb=0.1, ub=10) + brms::prior(constant(1), class=sigma))
+
+test_that("Gamma (truncated) prior is declared correctly in brms generated model and data", { mixstanvar_test(gamma, brms_gamma_trunc_args) })
+test_that("Gamma mixture (truncated) prior is declared correctly in brms generated model and data", { mixstanvar_test(gammaMix, brms_gamma_trunc_args) })
 
 # Here we approximate the samples using a multi-variante normal via a
 # moment based approxmation and compare this to the respective
@@ -170,10 +197,10 @@ brms_mvn_4_args <- list(formula=brms::bf(y ~ 1 + l1 + l2 + l3, family=brms::brms
                         data=data.frame(y=1, l1=0, l2=0, l3=0),
                         prior=brms::prior(mixmvnorm(prior_w, prior_m, prior_sigma_L), class=b) + brms::prior(constant(1), class=sigma))
 
-test_that("Multivariate normal (4D) is correct for brms sampled prior", mixstanvar_simul_mv_test(mvnorm_single_4, brms_mvn_4_args, eps))
-test_that("Multivariate normal (4D) prior is declared correctly in brms generated model and data", mixstanvar_test_mvnormMix(mvnorm_single_4, brms_mvn_4_args))
-test_that("Multivariate normal with heavy (4D) tails is correct for brms sampled prior", mixstanvar_simul_mv_test(mvnorm_heavy_4, brms_mvn_4_args, eps))
-test_that("Multivariate normal with heavy (4D) tails is declared correctly in brms generated model and data", mixstanvar_test_mvnormMix(mvnorm_heavy_4, brms_mvn_4_args))
+test_that("Multivariate normal (4D) is correct for brms sampled prior", { mixstanvar_simul_mv_test(mvnorm_single_4, brms_mvn_4_args, eps) })
+test_that("Multivariate normal (4D) prior is declared correctly in brms generated model and data", { mixstanvar_test_mvnormMix(mvnorm_single_4, brms_mvn_4_args) })
+test_that("Multivariate normal with heavy (4D) tails is correct for brms sampled prior", { mixstanvar_simul_mv_test(mvnorm_heavy_4, brms_mvn_4_args, eps) })
+test_that("Multivariate normal with heavy (4D) tails is declared correctly in brms generated model and data", { mixstanvar_test_mvnormMix(mvnorm_heavy_4, brms_mvn_4_args) })
 
 mvnorm_single_2 <- mixmvnorm(c(1, m1[1:2], 5),
                              param="mn", sigma=S[1:2,1:2])
@@ -187,7 +214,7 @@ brms_mvn_2_args <- list(formula=brms::bf(y ~ 1 + l1, family=brms::brmsfamily("ga
                         data=data.frame(y=1, l1=0),
                         prior=brms::prior(mixmvnorm(prior_w, prior_m, prior_sigma_L), class=b) + brms::prior(constant(1), class=sigma))
 
-test_that("Multivariate normal (2D) is correct for brms sampled prior", mixstanvar_simul_mv_test(mvnorm_single_2, brms_mvn_2_args, eps))
-test_that("Multivariate normal (2D) is declared correctly in brms generated model and data", mixstanvar_test_mvnormMix(mvnorm_single_2, brms_mvn_2_args))
-test_that("Multivariate normal with heavy (2D) tails is correct for brms sampled prior", mixstanvar_simul_mv_test(mvnorm_heavy_2, brms_mvn_2_args, eps))
-test_that("Multivariate normal with heavy (2D) is declared correctly in brms generated model and data", mixstanvar_test_mvnormMix(mvnorm_heavy_2, brms_mvn_2_args))
+test_that("Multivariate normal (2D) is correct for brms sampled prior", { mixstanvar_simul_mv_test(mvnorm_single_2, brms_mvn_2_args, eps) })
+test_that("Multivariate normal (2D) is declared correctly in brms generated model and data", { mixstanvar_test_mvnormMix(mvnorm_single_2, brms_mvn_2_args) })
+test_that("Multivariate normal with heavy (2D) tails is correct for brms sampled prior", { mixstanvar_simul_mv_test(mvnorm_heavy_2, brms_mvn_2_args, eps) })
+test_that("Multivariate normal with heavy (2D) is declared correctly in brms generated model and data", { mixstanvar_test_mvnormMix(mvnorm_heavy_2, brms_mvn_2_args) })
diff --git a/tests/testthat/test-oc1S.R b/tests/testthat/test-oc1S.R
index 0d16cc5..910470c 100644
--- a/tests/testthat/test-oc1S.R
+++ b/tests/testthat/test-oc1S.R
@@ -1,4 +1,3 @@
-context("oc1S: 1-Sample Operating Characteristics")
 
 ## test the analytical OC function via brute force simulation
 set.seed(12354)
@@ -34,13 +33,13 @@ test_scenario <- function(oc_res, ref) {
     expect_true(all(abs(resA) < eps))
 }
 
-test_that("Classical NI design critical value", expect_true( abs(decision1S_boundary(prior, 155, decA) - c1) < eps))
+test_that("Classical NI design critical value", { expect_true( abs(decision1S_boundary(prior, 155, decA) - c1) < eps) })
 
 ## n set to give power 80% to detect 0 and type I error 5% for no
 ## better than theta_ni
-test_that("Classical NI design at target    sample size for OCs", test_scenario(oc1S(prior, 155, decA)(thetaA), c(1-beta, alpha)) )
-test_that("Classical NI design at increased sample size for OCs", test_scenario(oc1S(prior, 233, decA)(thetaA), c(1-0.08, alpha)) )
-test_that("Classical NI design at decreased sample size for OCs", test_scenario(oc1S(prior,  77, decA)(thetaA), c(1-0.45, alpha)) )
+test_that("Classical NI design at target    sample size for OCs", { test_scenario(oc1S(prior, 155, decA)(thetaA), c(1-beta, alpha)) })
+test_that("Classical NI design at increased sample size for OCs", { test_scenario(oc1S(prior, 233, decA)(thetaA), c(1-0.08, alpha)) })
+test_that("Classical NI design at decreased sample size for OCs", { test_scenario(oc1S(prior,  77, decA)(thetaA), c(1-0.45, alpha)) })
 
 ## now double criterion with indecision point (mean estimate must be
 ## lower than this)
@@ -57,17 +56,17 @@ thetaD <- c(theta_c, theta_ni)
 
 ## since theta_c == c1, both decision criteria are the same for n =
 ## 155
-test_that("Double criterion NI design at target    sample size for OCs, combined      ", test_scenario(oc1S(prior, 155, decComb)(thetaD), c(0.50, alpha)) )
-test_that("Double criterion NI design at target    sample size for OCs, stat criterion", test_scenario(oc1S(prior, 155, dec1)(thetaD),    c(0.50, alpha)) )
-test_that("Double criterion NI design at target    sample size for OCs, mean criterion", test_scenario(oc1S(prior, 155, dec2)(thetaD),    c(0.50, alpha)) )
+test_that("Double criterion NI design at target    sample size for OCs, combined      ", { test_scenario(oc1S(prior, 155, decComb)(thetaD), c(0.50, alpha)) })
+test_that("Double criterion NI design at target    sample size for OCs, stat criterion", { test_scenario(oc1S(prior, 155, dec1)(thetaD),    c(0.50, alpha)) })
+test_that("Double criterion NI design at target    sample size for OCs, mean criterion", { test_scenario(oc1S(prior, 155, dec2)(thetaD),    c(0.50, alpha)) })
 
 ## at an increased sample size only the mean criterion is active
-test_that("Double criterion NI design at increased sample size for OCs, combined      ", test_scenario(oc1S(prior, 233, decComb)(thetaD), c(0.50, 0.02)) )
-test_that("Double criterion NI design at increased sample size for OCs, mean criterion", test_scenario(oc1S(prior, 233, dec2)(thetaD),    c(0.50, 0.02)) )
+test_that("Double criterion NI design at increased sample size for OCs, combined      ", { test_scenario(oc1S(prior, 233, decComb)(thetaD), c(0.50, 0.02)) })
+test_that("Double criterion NI design at increased sample size for OCs, mean criterion", { test_scenario(oc1S(prior, 233, dec2)(thetaD),    c(0.50, 0.02)) })
 
 ## at a  decreased sample size only the stat criterion is active
-test_that("Double criterion NI design at decreased sample size for OCs, combined      ", test_scenario(oc1S(prior,  78, decComb)(thetaD), c(1-0.68, alpha)) )
-test_that("Double criterion NI design at decreased sample size for OCs, stat criterion", test_scenario(oc1S(prior,  78, dec1)(thetaD),    c(1-0.68, alpha)) )
+test_that("Double criterion NI design at decreased sample size for OCs, combined      ", { test_scenario(oc1S(prior,  78, decComb)(thetaD), c(1-0.68, alpha)) })
+test_that("Double criterion NI design at decreased sample size for OCs, stat criterion", { test_scenario(oc1S(prior,  78, dec1)(thetaD),    c(1-0.68, alpha)) })
 
 ## test type 1 error and correctness of critical values wrt to
 ## lower.tail=TRUE/FALSE
@@ -95,12 +94,12 @@ design_binaryB <- oc1S(beta_prior, 1000, dec1b)
 crit1 <- decision1S_boundary(beta_prior, 1000, dec1)
 crit1B <- decision1S_boundary(beta_prior, 1000, dec1b)
 posterior_binary <- function(r) postmix(beta_prior, r=r, n=1000)
-test_that("Binary type I error rate", test_scenario(design_binary(theta_ni), alpha))
-test_that("Binary crticial value, lower.tail=TRUE",  test_critical_discrete(crit1,  dec1,  posterior_binary))
-test_that("Binary crticial value, lower.tail=FALSE", test_critical_discrete(crit1B, dec1b, posterior_binary))
+test_that("Binary type I error rate", { test_scenario(design_binary(theta_ni), alpha) })
+test_that("Binary crticial value, lower.tail=TRUE",  { test_critical_discrete(crit1,  dec1,  posterior_binary) })
+test_that("Binary crticial value, lower.tail=FALSE", { test_critical_discrete(crit1B, dec1b, posterior_binary) })
 
-test_that("Binary boundary case, lower.tail=TRUE",  expect_numeric(design_binary( 1), lower=0, upper=1, finite=TRUE, any.missing=FALSE))
-test_that("Binary boundary case, lower.tail=FALSE", expect_numeric(design_binaryB(0), lower=0, upper=1, finite=TRUE, any.missing=FALSE))
+test_that("Binary boundary case, lower.tail=TRUE",  { expect_numeric(design_binary( 1), lower=0, upper=1, finite=TRUE, any.missing=FALSE) })
+test_that("Binary boundary case, lower.tail=FALSE", { expect_numeric(design_binaryB(0), lower=0, upper=1, finite=TRUE, any.missing=FALSE) })
 
 ## poisson case
 gamma_prior <- mixgamma(c(1, 2, 2))
@@ -111,6 +110,6 @@ design_poissonB <- oc1S(gamma_prior, 1000, dec_countB)
 pcrit1 <- decision1S_boundary(gamma_prior, 1000, dec_count)
 pcrit1B <- decision1S_boundary(gamma_prior, 1000, dec_countB)
 posterior_poisson <- function(m) postmix(gamma_prior, m=m/1000, n=1000)
-test_that("Poisson type I error rate", test_scenario(design_poisson(1), alpha) )
-test_that("Poisson critical value, lower.tail=TRUE",  test_critical_discrete(pcrit1,  dec_count,  posterior_poisson))
-test_that("Poisson critical value, lower.tail=FALSE", test_critical_discrete(pcrit1B, dec_countB, posterior_poisson))
+test_that("Poisson type I error rate", { test_scenario(design_poisson(1), alpha) } )
+test_that("Poisson critical value, lower.tail=TRUE",  { test_critical_discrete(pcrit1,  dec_count,  posterior_poisson) })
+test_that("Poisson critical value, lower.tail=FALSE", { test_critical_discrete(pcrit1B, dec_countB, posterior_poisson) })
diff --git a/tests/testthat/test-oc2S.R b/tests/testthat/test-oc2S.R
index 4255f1c..1978bef 100644
--- a/tests/testthat/test-oc2S.R
+++ b/tests/testthat/test-oc2S.R
@@ -1,4 +1,3 @@
-context("oc2S: 2-Sample Operating Characteristics")
 
 ## test the analytical OC function via brute force simulation
 set.seed(12354)
@@ -20,8 +19,7 @@ theta2 <- 0.5
 
 Nsim <- 1e4
 
-run_on_cran <- function()
-{
+run_on_cran <- function() {
     if (identical(Sys.getenv("NOT_CRAN"), "true")) {
         return(FALSE)
     }
@@ -87,7 +85,7 @@ test_that("Type I error is matching between MC and analytical computations in th
 test_that("Power is matching between MC and analytical computations in the normal mixture case", {
               skip_on_cran()
 
-              power   <- oc2S(prior1, prior2, N1, N2, decision2S(pcrit, qcrit))(theta1, theta2)
+              power   <- oc2S(prior1, prior2, N1, N2, decision2S(pcrit, qcrit), sigma1=sigma(prior1), sigma2=sigma(prior2))(theta1, theta2)
               powerMC <- oc2S_normal_MC(prior1, prior2, N1, N2, theta1, theta2, pcrit, qcrit)
               res <- 100 * abs(power - powerMC)
               expect_equal(sum(res > 2) , 0)
@@ -125,10 +123,10 @@ test_that("Gsponer et al. results match (normal end-point)", {
               ## spline function
 
               ## Table 1, probability for interim for success
-              oc$success <- oc2S(priorP, priorT, nP1, nT1, successCrit, Ngrid=1)(-49, -49-oc$delta)
+              oc$success <- oc2S(priorP, priorT, nP1, nT1, successCrit, Ngrid=1, sigma1=sigmaFixed, sigma2=sigmaFixed)(-49, -49-oc$delta)
 
               ## Table 1, probability for interim for futility
-              oc$futile <- oc2S(priorP, priorT, nP1, nT1, futilityCrit, Ngrid=1)(-49, -49-oc$delta)
+              oc$futile <- oc2S(priorP, priorT, nP1, nT1, futilityCrit, Ngrid=1, sigma1=sigmaFixed, sigma2=sigmaFixed)(-49, -49-oc$delta)
 
               ## Table 1, first three columns, page 74
               oc[-1] <- lapply(100*oc[-1], round, 1)
@@ -236,9 +234,9 @@ expect_equal_each <- function(test, expected) {
 ## design object, decision function, posterior function must return
 ## posterior after updatding the prior with the given value; we assume
 ## that the priors are the same for sample 1 and 2
-test_critical_discrete <- function(design, decision, posterior, y2) {
+test_critical_discrete <- function(boundary_design, decision, posterior, y2) {
     lower.tail <- attr(decision, "lower.tail")
-    crit <- design(y2=y2)
+    crit <- boundary_design(y2)
     post2 <- posterior(y2)
     if(lower.tail) {
         expect_equal(decision(posterior(crit-1), post2), 1)
@@ -263,13 +261,20 @@ decB <- decision2S(1-alpha, 0, lower.tail=FALSE)
 beta_prior <- mixbeta(c(1, 1, 1))
 if(!run_on_cran()) {
     design_binary  <- oc2S(beta_prior, beta_prior, 100, 100, dec)
+    boundary_design_binary  <- decision2S_boundary(beta_prior, beta_prior, 100, 100, dec)
     design_binaryB <- oc2S(beta_prior, beta_prior, 100, 100, decB)
+    boundary_design_binaryB <- decision2S_boundary(beta_prior, beta_prior, 100, 100, decB)
+} else {
+    design_binary <- function(...) { return(0.1) }
+    design_binaryB <- function(...) { return(0.1) }
+    boundary_design_binary <- function(...) { return(0.1) }
+    boundary_design_binaryB <- function(...) { return(0.1) }
 }
 posterior_binary <- function(r) postmix(beta_prior, r=r, n=100)
 p_test <- 1:9 / 10
 test_that("Binary type I error rate", { skip_on_cran(); test_scenario(design_binary(p_test, p_test), alpha) })
-test_that("Binary crticial value, lower.tail=TRUE", { skip_on_cran(); test_critical_discrete(design_binary, dec, posterior_binary, 30) })
-test_that("Binary crticial value, lower.tail=FALSE", { skip_on_cran(); test_critical_discrete(design_binaryB, decB, posterior_binary, 30) })
+test_that("Binary crticial value, lower.tail=TRUE", { skip_on_cran(); test_critical_discrete(boundary_design_binary, dec, posterior_binary, 30) })
+test_that("Binary crticial value, lower.tail=FALSE", { skip_on_cran(); test_critical_discrete(boundary_design_binaryB, decB, posterior_binary, 30) })
 test_that("Binary boundary case, lower.tail=TRUE",  { skip_on_cran(); expect_numeric(design_binary( 1, 1), lower=0, upper=1, finite=TRUE, any.missing=FALSE) })
 test_that("Binary boundary case, lower.tail=FALSE", { skip_on_cran(); expect_numeric(design_binaryB(0, 0), lower=0, upper=1, finite=TRUE, any.missing=FALSE) })
 
@@ -280,9 +285,11 @@ beta_prior1 <- mixbeta(c(1, 0.9, 1000), param="mn")
 beta_prior2 <- mixbeta(c(1, 0.1, 1000), param="mn")
 design_lower <- oc2S(beta_prior1, beta_prior2, 20, 20, dec)  ## always 0
 design_upper <- oc2S(beta_prior1, beta_prior2, 20, 20, decB) ## always 1
+boundary_design_lower <- decision2S_boundary(beta_prior1, beta_prior2, 20, 20, dec)  ## always 0
+boundary_design_upper <- decision2S_boundary(beta_prior1, beta_prior2, 20, 20, decB) ## always 1
 
-test_that("Binary case, no decision change, lower.tail=TRUE, critical value", { skip_on_cran(); expect_equal_each(design_lower(y2=0:20), -1) })
-test_that("Binary case, no decision change, lower.tail=FALSE, critical value", { skip_on_cran(); expect_equal_each(design_upper(y2=0:20), 21) })
+test_that("Binary case, no decision change, lower.tail=TRUE, critical value", { skip_on_cran(); expect_equal_each(boundary_design_lower(0:20), -1) })
+test_that("Binary case, no decision change, lower.tail=FALSE, critical value", { skip_on_cran(); expect_equal_each(boundary_design_upper(0:20), 21) })
 test_that("Binary case, no decision change, lower.tail=TRUE, frequency=0", { skip_on_cran(); expect_equal_each(design_lower(p_test, p_test), 0.0) })
 test_that("Binary case, no decision change, lower.tail=FALSE, frequency=1",  { skip_on_cran(); expect_equal_each(design_upper(p_test, p_test), 1.0) })
 
@@ -290,10 +297,17 @@ test_that("Binary case, no decision change, lower.tail=FALSE, frequency=1",  { s
 if(!run_on_cran()) {
     design_lower_rev <- oc2S(beta_prior2, beta_prior1, 20, 20, dec)  ## always 1
     design_upper_rev <- oc2S(beta_prior2, beta_prior1, 20, 20, decB) ## always 0
+    boundary_design_lower_rev <- decision2S_boundary(beta_prior2, beta_prior1, 20, 20, dec)  ## always 1
+    boundary_design_upper_rev <- decision2S_boundary(beta_prior2, beta_prior1, 20, 20, decB) ## always 0
+} else {
+    design_lower_rev <- function(...) return(1)
+    design_upper_rev <- function(...) return(0)
+    boundary_design_lower_rev <- function(...) return(1)
+    boundary_design_upper_rev <- function(...) return(0)
 }
 
-test_that("Binary case, no decision change (reversed), lower.tail=TRUE, critical value",  { skip_on_cran(); expect_equal_each(design_lower_rev(y2=0:20), 20) })
-test_that("Binary case, no decision change (reversed), lower.tail=FALSE, critical value", { skip_on_cran(); expect_equal_each(design_upper_rev(y2=0:20), -1) })
+test_that("Binary case, no decision change (reversed), lower.tail=TRUE, critical value",  { skip_on_cran(); expect_equal_each(boundary_design_lower_rev(0:20), 20) })
+test_that("Binary case, no decision change (reversed), lower.tail=FALSE, critical value", { skip_on_cran(); expect_equal_each(boundary_design_upper_rev(0:20), -1) })
 test_that("Binary case, no decision change (reversed), lower.tail=TRUE, frequency=0",  { skip_on_cran(); expect_equal_each(design_lower_rev(p_test, p_test), 1.0) })
 test_that("Binary case, no decision change (reversed), lower.tail=FALSE, frequency=1",  { skip_on_cran(); expect_equal_each(design_upper_rev(p_test, p_test), 0.0) })
 test_that("Binary case, log-link", {
@@ -346,11 +360,13 @@ gamma_prior <- mixgamma(c(1, 2, 2))
 
 design_poisson  <- oc2S(gamma_prior, gamma_prior, 100, 100, dec)
 design_poissonB <- oc2S(gamma_prior, gamma_prior, 100, 100, decB)
+boundary_design_poisson  <- decision2S_boundary(gamma_prior, gamma_prior, 100, 100, dec)
+boundary_design_poissonB <- decision2S_boundary(gamma_prior, gamma_prior, 100, 100, decB)
 posterior_poisson <- function(m) postmix(gamma_prior, m=m/100, n=100)
 lambda_test <- seq(0.5, 1.3, by=0.1)
 test_that("Poisson type I error rate", { skip_on_cran(); test_scenario(design_poisson(lambda_test, lambda_test), alpha) })
-test_that("Poisson crticial value, lower.tail=TRUE", { skip_on_cran(); test_critical_discrete(design_poisson, dec, posterior_poisson, 90) })
-test_that("Poisson crticial value, lower.tail=FALSE", { skip_on_cran(); test_critical_discrete(design_poissonB, decB, posterior_poisson, 90) })
+test_that("Poisson crticial value, lower.tail=TRUE", { skip_on_cran(); test_critical_discrete(boundary_design_poisson, dec, posterior_poisson, 90) })
+test_that("Poisson crticial value, lower.tail=FALSE", { skip_on_cran(); test_critical_discrete(boundary_design_poissonB, decB, posterior_poisson, 90) })
 ## 22 Nov 2017: disabled test as we trigger always calculation of the
 ## boundaries as of now.
 ##test_that("Poisson results cache expands", {
@@ -413,8 +429,8 @@ test_that("Normal OC 2-sample avoids undefined behavior, example 1", {
     dec  <- decision2S(1-alpha, 0, lower.tail=FALSE)
     n <- 58
     k <- 2
-    design_map  <- oc2S(prior_flat, map_ref, n, n/k, dec)
-    design_map_2  <- oc2S(prior_flat, map_ref, n, n/k, dec)
+    design_map  <- oc2S(prior_flat, map_ref, n, n/k, dec, sigma1=sigma_ref, sigma2=sigma_ref)
+    design_map_2  <- oc2S(prior_flat, map_ref, n, n/k, dec, sigma1=sigma_ref, sigma2=sigma_ref)
 
     x <- seq(-2.6, -1.6, by=0.1)
     expect_numeric(design_map(x, x), lower=0, upper=1, any.missing=FALSE)
diff --git a/tests/testthat/test-pos1S.R b/tests/testthat/test-pos1S.R
index 45c620b..151c469 100644
--- a/tests/testthat/test-pos1S.R
+++ b/tests/testthat/test-pos1S.R
@@ -1,5 +1,3 @@
-context("pos1S: 1-Sample Probability of Success")
-
 ## test the analytical OC function via brute force simulation
 set.seed(12354)
 
@@ -50,19 +48,18 @@ test_pos1S <- function(prior, ia_dist, n, dec, decU) {
 }
 
 
-test_that("Normal PoS 1 sample function matches MC integration of CPO", test_pos1S(prior, ia_dist, n1, decA, decAU))
+test_that("Normal PoS 1 sample function matches MC integration of CPO", { test_pos1S(prior, ia_dist, n1, decA, decAU) })
 
 beta_prior <- mixbeta(c(1, 1, 1))
 beta_ia <- postmix(beta_prior, r=20, n=50)
-test_that("Binomial PoS 1 sample function matches MC integration of CPO", test_pos1S(beta_prior, beta_ia, n1, decA, decAU))
+test_that("Binomial PoS 1 sample function matches MC integration of CPO", { test_pos1S(beta_prior, beta_ia, n1, decA, decAU) })
 
 gamma_prior <- mixgamma(c(1, 1, 1), param="mn")
 dec_count  <- decision1S(1-alpha, 1, lower.tail=TRUE)
 dec_countU  <- decision1S(1-alpha, 1, lower.tail=FALSE)
 gamma_ia <- postmix(gamma_prior, m=0.9, n=40)
-test_that("Poisson PoS 1 sample function matches MC integration of CPO", test_pos1S(gamma_prior, gamma_ia, n1, dec_count, dec_countU))
+test_that("Poisson PoS 1 sample function matches MC integration of CPO", { test_pos1S(gamma_prior, gamma_ia, n1, dec_count, dec_countU) })
 
 prior_unit_inf <- mixnorm(c(1, 0, 1), sigma=s, param="mn")
 post_ia_unit_inf <- postmix(prior_unit_inf, m=-1, n=162)
-test_pos1S(prior_unit_inf, post_ia_unit_inf, 459-162, decA, decAU)
-test_that("Normal PoS 1 sample function matches MC integration of CPO (more extreme case)", test_pos1S(prior_unit_inf, post_ia_unit_inf, 459-162, decA, decAU))
+test_that("Normal PoS 1 sample function matches MC integration of CPO (more extreme case)", { test_pos1S(prior_unit_inf, post_ia_unit_inf, 459-162, decA, decAU) })
diff --git a/tests/testthat/test-pos2S.R b/tests/testthat/test-pos2S.R
index 2f52a97..b6eef52 100644
--- a/tests/testthat/test-pos2S.R
+++ b/tests/testthat/test-pos2S.R
@@ -1,4 +1,3 @@
-context("pos2S: 2-Sample Probability of Success")
 
 ## test the analytical OC function via brute force simulation
 set.seed(12354)
@@ -51,11 +50,12 @@ test_pos2S <- function(prior1, prior2, ia_dist1, ia_dist2, n1, n2, dec, decU) {
 ia1 <- postmix(prior1, m=0.2, se=s/sqrt(15))
 ia2 <- postmix(prior2, m=0, se=s/sqrt(15))
 
-test_that("Normal PoS 2 sample function matches MC integration of CPO",
-          test_pos2S(prior1, prior2,
-                     ia1, ia2,           
-                     N1, N2,
-                     dec, decU))
+test_that("Normal PoS 2 sample function matches MC integration of CPO",{
+    test_pos2S(prior1, prior2,
+               ia1, ia2,           
+               N1, N2,
+               dec, decU)
+})
 
 ## also run a MC comparison
 pos2S_normal_MC <- function(prior1, prior2, N1, N2, dtheta1, dtheta2, pcrit=0.975, qcrit=0) {
@@ -98,11 +98,12 @@ test_that("Normal PoS 2 sample function matches MC integration",
 beta_prior <- mixbeta(c(1, 1, 1))
 beta_ia1 <- postmix(beta_prior, r=20, n=50)
 beta_ia2 <- postmix(beta_prior, r=30, n=50)
-test_that("Binomial PoS 2 sample function matches MC integration of CPO",
+test_that("Binomial PoS 2 sample function matches MC integration of CPO", {
           test_pos2S(beta_prior, beta_prior,
                      beta_ia1, beta_ia2,           
                      N1, N2,
-                     dec, decU))
+                     dec, decU)
+})
 
 
 gamma_prior <- mixgamma(c(1, 1, 1), param="mn")
@@ -111,9 +112,26 @@ dec_count  <- decision2S(1-alpha, 0, lower.tail=TRUE)
 dec_countU  <- decision2S(1-alpha, 0, lower.tail=FALSE)
 gamma_ia1 <- postmix(gamma_prior, m=0.7, n=60)
 gamma_ia2 <- postmix(gamma_prior, m=1.2, n=60)
-test_that("Poisson PoS 2 sample function matches MC integration of CPO",
+test_that("Poisson PoS 2 sample function matches MC integration of CPO", {
           test_pos2S(gamma_prior, gamma_prior,
                      gamma_ia1, gamma_ia2,           
                      N1, N2,
-                     dec_count, dec_countU))
-
+                     dec_count, dec_countU)
+})
+
+test_that("Binomial PoS 2 with IA returns results", {
+    ## reported by user
+    successCrit  <- decision2S(c(0.9), c(0), lower.tail=FALSE)
+    n0 <- 50
+    n <- 100
+    n_alt <- 140
+    neutr_prior <- mixbeta(c(1,1/3,1/3))
+    post_placeboIA <- postmix(neutr_prior, r=13, n=n0)
+    post_treatIA   <- postmix(neutr_prior, r=3, n=n0)
+    # Criterion for PPoS at IA  
+    pos_final <- pos2S(post_treatIA, post_placeboIA, n-n0, n-n0, successCrit) 
+    pos_final_alt <- pos2S(post_treatIA, post_placeboIA, n_alt-n0, n_alt-n0, successCrit) 
+    #Predictive proba of success at the end
+    expect_number(pos_final_alt(post_treatIA, post_placeboIA), na.ok=FALSE, lower=0, upper=1, finite=TRUE, null.ok=FALSE)
+    expect_number(pos_final(post_treatIA, post_placeboIA), na.ok=FALSE, lower=0, upper=1, finite=TRUE, null.ok=FALSE)
+})
diff --git a/tests/testthat/test-postmix.R b/tests/testthat/test-postmix.R
index c2a8ad8..2492979 100644
--- a/tests/testthat/test-postmix.R
+++ b/tests/testthat/test-postmix.R
@@ -1,4 +1,3 @@
-context("postmix: Posterior Mixture Distribution")
 
 norm <- mixnorm(c(1, 0, 0.5), sigma=1)
 
diff --git a/tests/testthat/test-preddist.R b/tests/testthat/test-preddist.R
index 635eebc..0da4b37 100644
--- a/tests/testthat/test-preddist.R
+++ b/tests/testthat/test-preddist.R
@@ -1,4 +1,3 @@
-context("preddist: Mixture Predictive Distribution")
 
 ## check that predictive distributions hold what they promise,
 ## i.e. that they describe the sum of n new data points.
@@ -71,11 +70,11 @@ preddist_cmp <- function(mix,  n, n_rng, N=Nsamp, qntls=p_quants, stat=c("sum",
 }
 
 
-test_that("Predictive for a beta evaluates correctly (binary)", preddist_cmp(beta, n, Curry(rbinom, n=n, size=1)))
-test_that("Predictive for a beta mixture evaluates correctly (binary)", preddist_cmp(betaMix, n, Curry(rbinom, n=n, size=1)))
+test_that("Predictive for a beta evaluates correctly (binary)", { preddist_cmp(beta, n, Curry(rbinom, n=n, size=1)) })
+test_that("Predictive for a beta mixture evaluates correctly (binary)", { preddist_cmp(betaMix, n, Curry(rbinom, n=n, size=1)) })
 
-test_that("Predictive for a gamma evaluates correctly (poisson)", preddist_cmp(gamma, n, Curry(rpois, n=n)))
-test_that("Predictive for a gamma mixture evaluates correctly (poisson)", preddist_cmp(gammaMix, n, Curry(rpois, n=n), Teps=1E-1))
+test_that("Predictive for a gamma evaluates correctly (poisson)", { preddist_cmp(gamma, n, Curry(rpois, n=n)) })
+test_that("Predictive for a gamma mixture evaluates correctly (poisson)", { preddist_cmp(gammaMix, n, Curry(rpois, n=n), Teps=1E-1) })
 
-test_that("Predictive for a normal evaluates correctly (normal)", preddist_cmp(norm, n, Curry(rnorm, n=n, sd=sigma(norm)), stat="mean"))
-test_that("Predictive for a normal mixture evaluates correctly (normal)", preddist_cmp(normMix, n, Curry(rnorm, n=n, sd=sigma(normMix)), stat="mean", Teps=1E-1))
+test_that("Predictive for a normal evaluates correctly (normal)", { preddist_cmp(norm, n, Curry(rnorm, n=n, sd=sigma(norm)), stat="mean") })
+test_that("Predictive for a normal mixture evaluates correctly (normal)", { preddist_cmp(normMix, n, Curry(rnorm, n=n, sd=sigma(normMix)), stat="mean", Teps=1E-1) })
diff --git a/tests/testthat/test-utils.R b/tests/testthat/test-utils.R
index 2d0e391..81ea59e 100644
--- a/tests/testthat/test-utils.R
+++ b/tests/testthat/test-utils.R
@@ -1,10 +1,9 @@
-context("Utilities: predict")
 
 ## run the example from predict.gMAP
-suppressMessages(example("predict.gMAP", package="RBesT", echo=FALSE, ask=FALSE, verbose=FALSE))
+source_example("predict_gMAP.R")
 
 ## check that we got for each input data item a prediction
-test_that("correct # of predictions are generated", expect_equal(nrow(map$data), ncol(samp)))
+test_that("correct # of predictions are generated", { expect_equal(nrow(map$data), ncol(samp)) })
 
 ## check that the predictive distribution has a variance which is
 ## larger in accordance to the betwee-trial heterogeniety (needs to be
@@ -28,7 +27,7 @@ test_that("variances have correct ordering", {
 
 
 ## new prediction was done for a single data item
-test_that("correct # of new predictions are generated", expect_equal(ncol(pred_new), 1))
+test_that("correct # of new predictions are generated", { expect_equal(ncol(pred_new), 1) } )
 
 ## must have larger sd than between-trial alone (on link scale)
 test_that("predictive variances have correct ordering",{
@@ -46,8 +45,6 @@ test_that("predictive distributions for the same study & covariate must match ex
     expect_equal(post_trans[,1], post_trans[,2])
 })
 
-context("Utilities: (auto)mixfit")
-
 test_that("automixfit attempts K=4 different models and returns best fitting", {
               auto_map <- automixfit(map, Nc=1:4, k=6)
               models <- attr(auto_map, "models")
@@ -75,16 +72,14 @@ test_that("mixfit for prediction handles response and link scale", {
           })
 
 
-context("Utilities: mixcombine")
-
-example("mixcombine", package="RBesT", echo=FALSE, ask=FALSE, verbose=FALSE)
+source_example("mixcombine.R")
 
 test_that("combination of mixtures", {
               m1 <- mixcombine(bm, unif, weight=c(9, 1))
               m2 <- mixcombine(bm, unif, unif, weight=c(8, 1, 1))
-              expect_equivalent(m1[1,], c(bm[1,] - 0.1/2, 0.1))
-              expect_equivalent(m1[2:3,1:2], bm[2:3,1:2])
-              expect_equivalent(m2[2:3,1:2], bm[2:3,1:2])
+              expect_equal(m1[1,], c(bm[1,] - 0.1/2, 0.1), ignore_attr=TRUE)
+              expect_equal(m1[2:3,1:2], bm[2:3,1:2], ignore_attr=TRUE)
+              expect_equal(m2[2:3,1:2], bm[2:3,1:2], ignore_attr=TRUE)
           })
 
 test_that("throws an error if more weights than mixtures given", {
@@ -98,50 +93,47 @@ test_that("combination of normal mixtures without default sigma works", {
               expect_true(ncol(norm_ui_mix) == 2)
           })
 
-context("Utilities: robustify")
-
-example("robustify", package="RBesT", echo=FALSE, ask=FALSE, verbose=FALSE)
+source_example("robustify.R")
 
 test_that("beta mixture is robustified with Beta(1,1)", {
               expect_equal(ncol(bmix)+1, ncol(rbmix))
-              expect_equivalent(rbmix[,ncol(rbmix)], c(0.1, 1, 1))
+              expect_equal(rbmix[,ncol(rbmix)], c(0.1, 1, 1), ignore_attr=TRUE)
           })
 
 test_that("beta mixture is robustified with Beta(0.5,0.5)", {
-              rbmix2 <- robustify(bmix, w=0.1, n=0)
+              rbmix2 <- robustify(bmix, w=0.1, n=0, mean=0.5)
               expect_equal(ncol(bmix)+1, ncol(rbmix2))
-              expect_equivalent(rbmix2[,ncol(rbmix2)], c(0.1, 0.5, 0.5))
+              expect_equal(rbmix2[,ncol(rbmix2)], c(0.1, 0.5, 0.5), ignore_attr=TRUE)
           })
 
 test_that("gamma mixture is robustified with n=1 equivalent prior", {
               m <- summary(gmnMix)["mean"]
               nr <- ncol(rgmnMix)
-              expect_equivalent(rgmnMix[[nr, rescale=TRUE]], mixgamma(c(1, m, 1), param="mn"))
+              expect_equal(rgmnMix[[nr, rescale=TRUE]], mixgamma(c(1, m, 1), param="mn"), ignore_attr=TRUE)
               expect_equal(rgmnMix[1,nr], 0.1)
           })
 
 test_that("gamma mixture is robustified with n=5 equivalent prior", {
               m <- summary(gmnMix)["mean"]
-              rgmnMix2 <- robustify(gmnMix, w=0.1, n=5)
+              rgmnMix2 <- robustify(gmnMix, w=0.1, n=5, mean=2)
               nr <- ncol(rgmnMix2)
-              expect_equivalent(rgmnMix2[[nr, rescale=TRUE]], mixgamma(c(1, m, 5), param="mn"))
+              expect_equal(rgmnMix2[[nr, rescale=TRUE]], mixgamma(c(1, m, 5), param="mn"), ignore_attr=TRUE)
               expect_equal(rgmnMix2[1,nr], 0.1)
           })
 
 test_that("normal mixture is robustified with n=1 equivalent prior", {
               nr <- ncol(rnMix)
-              expect_equivalent(rnMix[[nr, rescale=TRUE]], mixnorm(c(1, 0, 1), param="mn", sigma=sigma(nm)))
+              expect_equal(rnMix[[nr, rescale=TRUE]], mixnorm(c(1, 0, 1), param="mn", sigma=sigma(nm)), ignore_attr=TRUE)
               expect_equal(rnMix[1,nr], 0.1)
           })
 
 test_that("normal mixture is robustified with n=5 equivalent prior", {
               rnMix2 <- robustify(nm, w=0.1, mean=0, n=5, sigma=sigma(nm))
               nr <- ncol(rnMix2)
-              expect_equivalent(rnMix2[[nr, rescale=TRUE]], mixnorm(c(1, 0, 5), param="mn", sigma=sigma(nm)))
+              expect_equal(rnMix2[[nr, rescale=TRUE]], mixnorm(c(1, 0, 5), param="mn", sigma=sigma(nm)), ignore_attr=TRUE)
               expect_equal(rnMix2[1,nr], 0.1)
           })
 
-context("Utilities: Plotting of Mixtures")
 test_that("plotting of normal mixtures without default sigma works", {
               norm_ui <- mixnorm(c(1, 0, 2))
               norm_mix_ui <- mixcombine(norm_ui, norm_ui, weight=c(0.5,0.5))
@@ -149,9 +141,7 @@ test_that("plotting of normal mixtures without default sigma works", {
               expect_true(inherits(pl, "ggplot"))
           })
 
-context("Utilities: Mixture Effective Sample Size")
-
-example("ess", package="RBesT", echo=FALSE, ask=FALSE, verbose=FALSE)
+source_example("ess.R")
 
 test_that("conjugate beta case matches canonical formula", {
               expect_equal(a+b, ess(prior, "moment"))
@@ -173,7 +163,7 @@ test_that("ess elir for beta mixtures gives a warning for a<1 & b<1 densities",
 
 test_that("ess elir for normal mixtures returns correct values", {
     mix <- mixnorm( inf1=c(0.5026,-191.1869,127.4207),inf2=c(0.2647,-187.5895,31.6130),inf3=c(0.2326,-184.7445,345.3849), sigma=270.4877)
-    expect_gt(ess(mix), 0)
+    expect_gt(ess(mix, sigma=270.4877), 0)
 })
 
 test_that("moment matching for beta mixtures is correct", {
@@ -297,3 +287,31 @@ test_that("ESS elir is predictively consistent for gamma mixtures (Poisson likel
     gmixP <- mixgamma(rob=c(0.3, 20, 4), inf=c(0.7, 50, 10), likelihood="poisson")
     elir_predictive_consistent(gmixP, m=1E2, Nsim=1E3, seed=355435, stat="m", n=1E2)
 })
+
+test_that("ess elir for problematic beta mixtures gives correct result 1", {
+    ## by user reported beta mixture density which triggers this erros
+    ## with RBesT 1.7.2 & 1.7.3 (others not tested):
+    ## Error in if (all(dgl < 0) || all(dgl > 0)) { : 
+    ##   missing value where TRUE/FALSE needed
+
+    mixmat <- matrix(c(0.06429517, 0.03301215, 0.00269268, 0.90000000,
+                       437.32302999, 64.04211307, 5.92543558, 1.00000000,
+                       10.71709277, 2.14157953, 1.00000001, 1.00000000), byrow=TRUE, ncol=4)
+
+    mixb <- do.call(mixbeta, apply(mixmat,2,c,simplify=FALSE))
+
+    expect_double(ess(mixb), lower=0, finite=TRUE, any.missing=FALSE, len=1)
+})
+ 
+test_that("ess elir for problematic beta mixtures gives correct result 2", {
+    mixmat <- matrix(c(0.7237396, 0.1665037,  0.1097567,
+                       53.3721902, 44.3894573,  9.8097062,
+                       1.4301638,  4.3842200,  1.8492197
+                       ), byrow=TRUE, ncol=3)
+
+    mixb <- do.call(mixbeta, apply(mixmat,2,c,simplify=FALSE))
+    
+    expect_double(ess(robustify(mixb, 0.05, 0.5)), lower=0, finite=TRUE, any.missing=FALSE, len=1)
+    expect_double(ess(robustify(mixb, 0.95, 0.5)), lower=0, finite=TRUE, any.missing=FALSE, len=1)
+})
+ 
diff --git a/tools/make-ds.R b/tools/make-ds.R
index a7274c9..9fd1a33 100644
--- a/tools/make-ds.R
+++ b/tools/make-ds.R
@@ -43,7 +43,7 @@ make_internal_ds <- function() {
     calibration_meta["MD5"]  <- vals["MD5"]
 
     pkg_create_date  <- Sys.time()
-    pkg_sha <- "511a0f1"
+    pkg_sha <- "896a402"
 
     if (gsub("\\$", "", pkg_sha) == "Format:%h") {
         pkg_sha <- system("git rev-parse --short HEAD", intern=TRUE)
diff --git a/vignettes/articles/introduction_normal.Rmd b/vignettes/articles/introduction_normal.Rmd
index bf49a9e..b102d97 100644
--- a/vignettes/articles/introduction_normal.Rmd
+++ b/vignettes/articles/introduction_normal.Rmd
@@ -237,11 +237,13 @@ ocI <- rbind(data.frame(cfb_truth=cfb_truth, typeI=typeI1,
              data.frame(cfb_truth=cfb_truth, typeI=typeI4,
                         design="40:20 with robust prior for placebo")
 )
-qplot(cfb_truth, typeI, data=ocI, colour=design, geom="line", main="Type I Error") +
+ggplot(ocI, aes(cfb_truth, typeI, colour=design)) +
+    geom_line() +
+    ggtitle("Type I Error") +
     xlab(expression(paste('True value of change from baseline ', mu[act] == mu[pbo]))) +
-        ylab('Type I error') +
-            coord_cartesian(ylim=c(0,0.2)) +
-  theme(legend.justification=c(1,1),legend.position=c(0.95,0.85))
+    ylab('Type I error') +
+    coord_cartesian(ylim=c(0,0.2)) +
+    theme(legend.justification=c(1,1),legend.position=c(0.95,0.85))
 ```
 
 
@@ -274,17 +276,17 @@ ocP <- rbind(data.frame(cfb_truth1=cfb_truth1, cfb_truth2=cfb_truth2,
                         design="40:20 with non-robust prior for placebo"),
              data.frame(cfb_truth1=cfb_truth1, cfb_truth2=cfb_truth2,
                         delta=delta, power=power4,
-                        design="40:20 with robust prior for placebo")
-)
-
-qplot(delta, power, data=ocP, colour=design, geom="line", main="Power") +
-  xlab('True value of difference (act - pbo)')+ ylab('Power') +
-  scale_y_continuous(breaks=c(seq(0,1,0.2),0.9)) +
-  scale_x_continuous(breaks=c(seq(-80,0,20),-70)) +
-  geom_hline(yintercept=0.9,linetype=2) + 
-  geom_vline(xintercept=-70,linetype=2) +
-  theme(legend.justification=c(1,1),legend.position=c(0.95,0.85))
- 
+                        design="40:20 with robust prior for placebo"))
+
+ggplot(ocP, aes(delta, power, colour=design)) +
+    geom_line() +
+    ggtitle("Power") + 
+    xlab('True value of difference (act - pbo)')+ ylab('Power') +
+    scale_y_continuous(breaks=c(seq(0,1,0.2),0.9)) +
+    scale_x_continuous(breaks=c(seq(-80,0,20),-70)) +
+    geom_hline(yintercept=0.9,linetype=2) + 
+    geom_vline(xintercept=-70,linetype=2) +
+    theme(legend.justification=c(1,1),legend.position=c(0.95,0.85)) 
 ```
 
 
diff --git a/vignettes/introduction.Rmd b/vignettes/introduction.Rmd
index 2d811c3..5e6cbf9 100644
--- a/vignettes/introduction.Rmd
+++ b/vignettes/introduction.Rmd
@@ -214,8 +214,8 @@ ess_weight <- rbind(ess_weight,
                     data.frame(weight=c(0, 1),
                                ess=c(ess(map), ess(mixbeta(c(1,1,1))))))
 
-qplot(weight, ess, data=ess_weight, geom=c("point", "line"),
-      main="ESS of robust MAP for varying weight of robust component") +
+ggplot(ess_weight, aes(weight, ess)) + geom_point() + geom_line() +
+    ggtitle("ESS of robust MAP for varying weight of robust component") +
     scale_x_continuous(breaks=seq(0,  1, by=0.1)) +
     scale_y_continuous(breaks=seq(0, 40, by=5))
 ```
@@ -290,7 +290,7 @@ ocI <- rbind(data.frame(theta=theta, typeI=typeI_robust,    prior="robust"),
              data.frame(theta=theta, typeI=typeI_classic,   prior="uniform 24:24")
              )
 
-qplot(theta, typeI, data=ocI, colour=prior, geom="line", main="Type I Error")
+ggplot(ocI, aes(theta, typeI, colour=prior)) + geom_line() + ggtitle("Type I Error")
 ```
 
 Note that observing response rates greater that 50% is highly
@@ -304,8 +304,9 @@ Hence, it is resonable to restrict the response rates $\theta$ for
 which we evaluate the type I error to a a range of plausible values:
 
 ```{r}
-qplot(theta, typeI, data=subset(ocI, theta < 0.5), colour=prior, geom="line",
-      main="Type I Error - response rate restricted to plausible range")
+ggplot(ocI, aes(theta, typeI, colour=prior)) + geom_line() +
+    ggtitle("Type I Error - response rate restricted to plausible range") +
+    coord_cartesian(xlim=c(0, 0.5))
 ```
 
 ### Power
@@ -331,7 +332,8 @@ ocP <- rbind(data.frame(theta_active, theta_control, delta=delta, power=power_ro
              data.frame(theta_active, theta_control, delta=delta, power=power_classic,   prior="uniform 24:24")
              )
 
-qplot(delta, power, data=ocP, colour=prior, geom="line", main="Power")
+ggplot(ocP, aes(delta, power, colour=prior)) + geom_line() +
+    ggtitle("Power")
 
 ```
 
@@ -378,7 +380,7 @@ ocC <- rbind(data.frame(y2=treat_y2, y1_crit=crit_robust(treat_y2),    prior="ro
              data.frame(y2=treat_y2, y1_crit=crit_uniform(treat_y2),   prior="uniform")
              )
 
-qplot(y2, y1_crit, data=ocC, colour=prior, geom="step", main="Critical values y1(y2)")
+ggplot(ocC, aes(y2, y1_crit, colour=prior)) + geom_step() + ggtitle("Critical values y1(y2)")
 ```
 
 The graph shows that the decision will always be negative if there are