120 issues

R/dbetabinom.R

@@ -83,9 +83,20 @@ h_getBetamixPost <- function(x, n, par, weights) {
   assert_number(x, lower = 0, upper = n, finite = TRUE)
   assert_number(n, lower = 0, finite = TRUE)
   assert_matrix(par, min.rows = 1, max.cols = 2, mode = "numeric")
-  assert_numeric(weights, min.len = 0, len = nrow(par), finite = TRUE)
-  # We renormalize weights.
-  weights <- weights / sum(weights)
+  assert_numeric(
+    weights,
+    min.len = 0,
+    len = nrow(par),
+    finite = TRUE
+  )
+  assert_true(all(weights >= 0))
+  # Correcting weights that do not sum to 1
+  if (sum(weights) != 1) {
+    warning("Weights have been corrected. Advise to review allocated weights")
+    # We renormalize weights.
+    weights <- weights / sum(weights)
+  }
+
   # We now compute updated parameters.
   postPar <- par
   postPar[, 1] <- postPar[, 1] + x

R/postprob.R

@@ -68,7 +68,14 @@ postprobBeta <- function(x, n, p, a = 1, b = 1) {
 #'
 #' @example examples/postprob.R
 #' @export
-postprob <- function(x, n, p, parE = c(1, 1), weights, betamixPost, log.p = FALSE) {
+postprob <- function(
+    x,
+    n,
+    p,
+    parE = c(1, 1),
+    weights,
+    betamixPost,
+    log.p = FALSE) {
   if (missing(betamixPost)) {
     assert_flag(log.p)
     if (is.vector(parE)) {
@@ -79,7 +86,7 @@ postprob <- function(x, n, p, parE = c(1, 1), weights, betamixPost, log.p = FALS
     }
     assert_matrix(parE)
     if (missing(weights)) {
-      weights <- rep(1, nrow(parE))
+      weights <- rep(1 / nrow(parE), nrow(parE))
     }
     betamixPost <- lapply(
       x,

R/predprob.R

@@ -45,8 +45,7 @@
 #'
 #' @example examples/predprob.R
 #' @export
-predprob <- function(x, n, Nmax, p, thetaT, parE = c(1, 1),
-                     weights) {
+predprob <- function(x, n, Nmax, p, thetaT, parE = c(1, 1), weights) {
   # m = Nmax - n future observations
   m <- Nmax - n
   if (is.vector(parE)) {
@@ -55,18 +54,29 @@ predprob <- function(x, n, Nmax, p, thetaT, parE = c(1, 1),
     parE <- t(parE)
   }
   if (missing(weights)) {
-    weights <- rep(1, nrow(parE))
+    weights <- rep(1 / nrow(parE), nrow(parE))
   }
   betamixPost <- h_getBetamixPost(x = x, n = n, par = parE, weights = weights)
-
   density <- with(
     betamixPost,
     dbetabinomMix(x = 0:m, m = m, par = par, weights = weights)
   )
-  assert_numeric(density, lower = 0, upper = 1 + .Machine$double.eps, finite = TRUE, any.missing = FALSE)
+  assert_numeric(
+    density,
+    lower = 0,
+    upper = 1 + .Machine$double.eps,
+    finite = TRUE,
+    any.missing = FALSE
+  )
   assert_number(thetaT, lower = 0, upper = 1, finite = TRUE)
   # posterior probabilities to be above threshold p
-  posterior <- postprob(x = x + c(0:m), n = Nmax, p = p, parE = parE, weights = weights)
+  posterior <- postprob(
+    x = x + c(0:m),
+    n = Nmax,
+    p = p,
+    parE = parE,
+    weights = weights
+  )
   list(
     result = sum(density * (posterior > thetaT)),
     table = data.frame(

examples/postprob.R

@@ -10,21 +10,23 @@ postprob(x = 16, n = 23, p = 0.60, par = c(0.6, 0.4))
 # 2 component beta mixture prior :
 # i.e., P_E ~ 0.6*beta(0.6,0.4) + 0.4*beta(1,1) and Pr(P_E > p | data) = 0.823
 postprob(
-  x = 16, n = 23, p = 0.60,
-  par =
-    rbind(
-      c(0.6, 0.4),
-      c(1, 1)
-    ),
+  x = 16,
+  n = 23,
+  p = 0.60,
+  par = rbind(
+    c(0.6, 0.4),
+    c(1, 1)
+  ),
   weights = c(0.6, 0.4)
 )
 
 postprob(
-  x = 0:23, n = 23, p = 0.60,
-  par =
-    rbind(
-      c(0.6, 0.4),
-      c(1, 1)
-    ),
+  x = 0:23,
+  n = 23,
+  p = 0.60,
+  par = rbind(
+    c(0.6, 0.4),
+    c(1, 1)
+  ),
   weights = c(0.6, 0.4)
 )

examples/predprob.R

@@ -2,25 +2,41 @@
 # Nmax = 40, x = 16, n = 23, beta(0.6,0.4) prior distribution,
 # thetaT = 0.9. The control response rate is 60%:
 predprob(
-  x = 16, n = 23, Nmax = 40, p = 0.6, thetaT = 0.9,
+  x = 16,
+  n = 23,
+  Nmax = 40,
+  p = 0.6,
+  thetaT = 0.9,
   parE = c(0.6, 0.4)
 )
 
 # Lowering/Increasing the probability threshold thetaT of course increases
 # /decreases the predictive probability of success, respectively:
 predprob(
-  x = 16, n = 23, Nmax = 40, p = 0.6, thetaT = 0.8,
+  x = 16,
+  n = 23,
+  Nmax = 40,
+  p = 0.6,
+  thetaT = 0.8,
   parE = c(0.6, 0.4)
 )
 
 predprob(
-  x = 16, n = 23, Nmax = 40, p = 0.6, thetaT = 0.95,
+  x = 16,
+  n = 23,
+  Nmax = 40,
+  p = 0.6,
+  thetaT = 0.95,
   parE = c(0.6, 0.4)
 )
 
 # Mixed beta prior
 predprob(
-  x = 20, n = 23, Nmax = 40, p = 0.6, thetaT = 0.9,
+  x = 20,
+  n = 23,
+  Nmax = 40,
+  p = 0.6,
+  thetaT = 0.9,
   parE = rbind(c(1, 1), c(25, 15)),
-  weights = c(3, 1)
+  weights = c(0.5, 0.5)
 )