diff --git a/R/dbetabinom.R b/R/dbetabinom.R
index 3c92f18e..ff9c9de2 100644
--- a/R/dbetabinom.R
+++ b/R/dbetabinom.R
@@ -19,8 +19,6 @@
 #'  whether to return the log density value (not default).
 #' @return The density values of the beta-binomial distribution at `x`.
 #'
-#' @note `x`, `a` and `b` can be vectors.
-#'
 #' @example examples/dbetabinom.R
 #' @export
 dbetabinom <- function(x, m, a, b, log = FALSE) {
@@ -57,8 +55,6 @@ dbetabinom <- function(x, m, a, b, log = FALSE) {
 #'  whether to return the log density value (not default).
 #' @return The (log) density values of the mixture of beta-binomial distributions at `x`.
 #'
-#' @note `x` can be a vector.
-#'
 #' @example examples/dbetabinomMix.R
 #' @export
 dbetabinomMix <- function(x, m, par, weights, log = FALSE) {
@@ -89,6 +85,7 @@ dbetabinomMix <- Vectorize(dbetabinomMix, vectorize.args = "x")
 #'
 #' @example examples/getBetamixPost.R
 #' @export
+
 getBetamixPost <- function(x, n, par, weights) {
   ## check the format
   stopifnot(
@@ -125,6 +122,7 @@ getBetamixPost <- function(x, n, par, weights) {
 }
 
 
+
 #' Beta-mixture density function
 #'
 #' Note that `x` can be a vector.
diff --git a/R/postprob.R b/R/postprob.R
index 4fbb71b9..716f0730 100644
--- a/R/postprob.R
+++ b/R/postprob.R
@@ -35,7 +35,6 @@ postprobBeta <- function(x, n, p, a = 1, b = 1) {
   stats::pbeta(p, a + x, b + n - x, lower.tail = FALSE)
 }
 
-
 #' Posterior Probability of Efficacy Given Beta-Mixture Prior
 #'
 #' @description `r lifecycle::badge("experimental")`
diff --git a/R/postprobDist.R b/R/postprobDist.R
index 728752e5..26e96a28 100644
--- a/R/postprobDist.R
+++ b/R/postprobDist.R
@@ -2,135 +2,194 @@
 #' @include postprob.R
 NULL
 
-#' Compute the posterior probability with beta prior on SOC
-#'
-#' Using the approach by Thall and Simon (Biometrics, 1994), evaluate the
-#' posterior probability of having Pr(P_E > P_S + delta | data) (but see below
-#' for relative delta margin). Both for the new treatment E as well as for the
-#' SOC S data might be available. However the default is that no data is
-#' available for the SOC, corresponding to the single arm trial situation. Note
-#' that a uniform prior is the useful default for the treatment proportion,
-#' while in the single arm trial an informative prior on the SOC proportion is
-#' useful.
-#'
-#' Beta mixture prior can be specified for the treatment (\code{parE}
-#' and \code{weights} parameters) and control proportion (\code{parS} and
-#' \code{weightsS} parameters), see \code{\link{postprob}} for details. Note
-#' that being able to specify a beta mixture prior also on the control
-#' treatment is e.g. important for the futility decision making (see the
-#' \code{\link{oc2}} code).
-#'
-#' @param x number of successes (in the treatment group). Note that \code{x}
-#' can be a vector.
-#' @param n number of patients (in the treatment group)
-#' @param xS number of successes in the SOC group (default: 0)
-#' @param nS number of patients in the SOC group (default: 0)
-#' @param delta margin by which the response rate in the treatment group should
-#' be better than in the SOC group (default: 0)
-#' @param relativeDelta should the delta be relative? (not default). If this is
-#' \code{TRUE}, then a relative delta is used. This means we want to have
-#' response at least in delta proportion of the SOC non-responding patients.
-#' Non-responding patients rate is 1 - P_S, and we want to have P_S + (1 - P_S)
-#' * delta response rate (at least) in the treatment. That is, we evaluate the
-#' posterior probability Pr(P_E > P_S + (1 - P_S) * delta | data).
-#' @param parE the beta parameters matrix, with K rows and 2 columns,
-#' corresponding to the beta parameters of the K components. default is a
-#' uniform prior.
-#' @param weights the mixture weights of the beta mixture prior. Default are
-#' uniform weights across mixture components.
-#' @param parS beta parameters for the SOC group (default: uniform)
-#' @param weightsS weights for the SOC group (default: uniform)
-#' @return the posterior probability
+#' The Posterior Beta Mixture Integrand when `delta` is relative
+#'
+#' The helper function to generate Integrand function when `relative Delta = TRUE`.
+#'
+#' @typed delta : number
+#'  margin by which the response rate in the treatment group should
+#'  be better than in the SOC group. Must be >= `0`. See note.
+#' @typed p_s : number
+#'  probability of success or response rate of standard of care or `SOC` group.
+#' @typed activeBetamixPost : list
+#'  a list of posterior parameters of a beta-mixture-binomial distribution with generic names
+#'  `par` and `weights`. See `[getBetaMix()]`.
+#' @typed controlBetamixPost : list
+#'  a list of posterior parameters of a beta-mixture-binomial distribution with generic names
+#'  `par` and `weights`. See `[getBetaMix()]`.
+#'
+#' @return Function that is an argument for `[stats::integrate()]`.
+#'
+#' @keywords internal
+#'
+h_integrand_relDelta <- function(p_s, delta, activeBetamixPost, controlBetamixPost) {
+  cdf <- postprob(
+    x = 0, # Needed for Vectorize()
+    p = (1 - p_s) * delta + p_s,
+    betamixPost = activeBetamixPost
+  )
+  pdf <- with(
+    controlBetamixPost,
+    dbetaMix(x = p_s, par = par, weights = weights)
+  )
+  cdf * pdf
+}
+
+#' The Posterior Beta Mixture Integrand when Delta is absolute
+#'
+#' The helper function to generate Integrand function when `relative Delta = FALSE`,
+#' a default setting.
+#'
+#' @inheritParams h_integrand_relDelta
+#'
+#' @return Function that is an argument for `[stats::integrate()]`.
+#'
+#' @keywords internal
+#'
+h_integrand <- function(p_s, delta, activeBetamixPost, controlBetamixPost) {
+  cdf <- postprob(
+    x = 0, # Needed for Vectorize()
+    p = p_s + delta,
+    betamixPost = activeBetamixPost
+  )
+  pdf <- with(
+    controlBetamixPost,
+    dbetaMix(x = p_s, par = par, weights = weights)
+  )
+  cdf * pdf
+}
+
+#' Generating bounds for the Integration of Beta Mixture Posterior
+#'
+#' Using the quantile of the Beta Mixture Distribution from parameters given by standard of care `SOC` or
+#' experimental group `E` to determine bounds as inputs to `[stats::integrate()]`.
+#'
+#' @inheritParams h_integrand_relDelta
+#' @return Integrand function
+#'
+#' @keywords internal
+#'
+h_get_bounds <- function(controlBetamixPost) {
+  epsilon <- .Machine$double.xmin
+  with(
+    controlBetamixPost,
+    qbetaMix(
+      p = c(epsilon, 1 - epsilon),
+      par = par,
+      weights = weights
+    )
+  )
+}
+
+#' Compute the Posterior Probability with Beta Prior on `SOC`
+#'
+#' @description `r lifecycle::badge("experimental")`
+#'
+#' Using the approach by Thall and Simon (Biometrics, 1994), this evaluates the
+#' posterior probability of achieving superior response rate in the treatment group compared to standard of care (SOC).
+#' See note below for two formulations of the difference in response rates.
+#'
+#' @inheritParams h_integrand_relDelta
+#' @typed x : numeric
+#'  number of success counts in the treatment group.
+#' @typed n : number
+#'  number of patients in the treatment group.
+#' @typed xS : number
+#'  number of success counts in the SOC group.
+#' @typed nS : number
+#'  number of patients in the SOC group.
+#' @typed relativeDelta : flag
+#'  If `TRUE`, then a `relativeDelta` is used. Represents that a minimum
+#'  response rate in magnitude of `delta` of the SOC non-responding patients. See note.
+#' @typed parE : "`numeric` or `matrix`"
+#'  parameters for beta distribution. If it is a matrix, it needs to have 2 columns,
+#'  and each row corresponds to each component of a beta-mixture distribution
+#'  for the `E` group. See details.
+#' @typed weights : numeric
+#'  the non-negative mixture weights of the beta mixture prior for group `E`. See details.
+#' @typed parS : "`numeric` or `matrix`"
+#'  parameters for beta distribution. If it is a matrix, it needs to have 2 columns,
+#'  and each row corresponds to each component of a beta-mixture distribution
+#'  for the `S` group. See details.
+#' @typed weightsS : numeric
+#'  the non-negative mixture weights of the beta mixture prior for group `S`. See details.
+#' @typed epsilon : number
+#'  the smallest non-negative floating number to represent the lower bound for
+#'  the interval of integration.
+#' @return The posterior probability
+#'
+#' @note
+#'
+#' ## Delta :
+#'
+#' The desired improvement is denoted as `delta`. There are two options in using `delta`.
+#' The absolute case when `relativeDelta = FALSE` and relative as when `relativeDelta = TRUE`.
+#'
+#' 1. The absolute case is when we define an absolute delta, greater than `P_S`,
+#' the response rate of the `SOC` group such that
+#' the posterior is `Pr(P_E > P_S + delta | data)`.
+#'
+#' 2. In the relative case, we suppose that the treatment group's
+#' response rate is assumed to be greater than `P_S + (1-P_S) * delta` such that
+#' the posterior is `Pr(P_E > P_S + (1 - P_S) * delta | data)`.
+#'
+#' @details
+#'
+#' The beta mixture prior for the `E` arm requires argument `parE` and `weights`.
+#' The beta mixture prior for the `S` arm requires argument `parS` and `weightsS`.
+#' See `[postprob()]` for details.
+#'
+#' If a beta-mixture is used, by default, the weights are uniform across the components.
+#' Weights can exceed 1, to which the algorithm will normalize the weights such that all weights sum to 1.
 #'
 #' @example examples/postprobDist.R
 #' @export
-postprobDist <- function(x, n,
-                         xS = 0, nS = 0,
-                         delta = 0,
+postprobDist <- function(x,
+                         n,
+                         xS = 0,
+                         nS = 0,
+                         delta,
                          relativeDelta = FALSE,
                          parE = c(1, 1),
                          weights,
                          parS = c(1, 1),
                          weightsS) {
-  ## if parE is a vector => situation where there is only one component
   if (is.vector(parE)) {
-    ## check that it has exactly two entries
-    stopifnot(identical(length(parE), 2L))
-
-    ## and transpose to matrix with one row
+    assert_true(identical(length(parE), 2L))
     parE <- t(parE)
   }
-
-  ## if prior weights of the beta mixture are not supplied
-  if (missing(weights)) {
-    weights <- rep(1, nrow(parE))
-  }
-
-  ## if parS is a vector => situation where there is only one component
   if (is.vector(parS)) {
-    ## check that it has exactly two entries
-    stopifnot(identical(length(parS), 2L))
-
-    ## and transpose to matrix with one row
+    assert_true(identical(length(parS), 2L))
     parS <- t(parS)
   }
-
-  ## if prior weights of the beta mixture are not supplied
+  if (missing(weights)) {
+    weights <- rep(1, nrow(parE))
+  }
   if (missing(weightsS)) {
     weightsS <- rep(1, nrow(parS))
   }
-
-  ## compute updated beta parameters
+  assert_numeric(x, lower = 0, upper = n, finite = TRUE)
+  assert_number(n, lower = x, finite = TRUE)
+  assert_number(xS, lower = 0, upper = n, finite = TRUE)
+  assert_numeric(nS, lower = 0, finite = TRUE)
+  assert_number(delta, lower = 0, finite = TRUE)
+  assert_flag(relativeDelta)
+  assert_numeric(weights, lower = 0, finite = TRUE)
+  assert_numeric(weightsS, lower = 0, finite = TRUE)
+  assert_numeric(parE, lower = 0, finite = TRUE)
+  assert_numeric(parS, lower = 0, finite = TRUE)
   activeBetamixPost <- getBetamixPost(x = x, n = n, par = parE, weights = weights)
   controlBetamixPost <- getBetamixPost(x = xS, n = nS, par = parS, weights = weightsS)
-
-  ## use numerical integration to compute this probability, as given on p.338
-  ## in the article by Thall and Simon (1994):
-  integrand <-
-    if (relativeDelta) {
-      function(p) {
-        cdf <- postprob(
-          x = x,
-          p = (1 - delta) * p + delta,
-          betamixPost = activeBetamixPost
-        )
-
-        pdf <- with(
-          controlBetamixPost,
-          dbetaMix(x = p, par = par, weights = weights)
-        )
-
-        return(cdf * pdf)
-      }
-    } else {
-      function(p) {
-        cdf <- postprob(
-          x = x,
-          p = p + delta,
-          betamixPost = activeBetamixPost
-        )
-
-        pdf <- with(
-          controlBetamixPost,
-          dbetaMix(x = p, par = par, weights = weights)
-        )
-
-        return(cdf * pdf)
-      }
-    }
-
-  ## do the integration. be careful to cover the region where there can
-  ## really be any non-zero values. I.e. only integrate over the region where
-  ## the beta density of the control is non-zero.
-  epsilon <- 1e-13
-  bounds <- with(
-    controlBetamixPost,
-    qbetaMix(
-      p = c(epsilon, 1 - epsilon),
-      par = par,
-      weights = weights
-    )
-  )
+  assert_names(names(activeBetamixPost), identical.to = c("par", "weights"))
+  assert_names(names(controlBetamixPost), identical.to = c("par", "weights"))
+  if (relativeDelta) {
+    epsilon <- .Machine$double.xmin
+    integrand <- h_integrand_relDelta
+  } else {
+    epsilon <- .Machine$double.xmin
+    integrand <- h_integrand
+  }
+  bounds <- h_get_bounds(controlBetamixPost = controlBetamixPost)
   intRes <- integrate(
     f = integrand,
     lower =
@@ -142,11 +201,13 @@ postprobDist <- function(x, n,
       min(
         ifelse(relativeDelta, 1, 1 - delta),
         bounds[2]
-      )
+      ),
+    delta = delta,
+    activeBetamixPost = activeBetamixPost,
+    controlBetamixPost = controlBetamixPost
   )
-
   if (intRes$message == "OK") {
-    return(intRes$value)
+    intRes$value
   } else {
     stop(intRes$message)
   }
diff --git a/examples/postprobDist.R b/examples/postprobDist.R
index a7dac8a9..0c132b84 100644
--- a/examples/postprobDist.R
+++ b/examples/postprobDist.R
@@ -1,26 +1,37 @@
-## example similar to Lee and Liu:
-postprobDist(x = 16, n = 23, parE = c(0.6, 0.4), parS = c(0.6, 0.4))
+# An example similar to Lee and Liu (2008).
+postprobDist(
+  x = 16,
+  n = 23,
+  parE = c(0.6, 0.4),
+  parS = c(0.6, 0.4),
+  delta = 0.1,
+  relativeDelta = FALSE
+)
 
-## these two should give the same result:
+# For a sequence of success outcomes for Experimental arm.
 postprobDist(
-  x = 27, n = 34,
-  xS = 0, nS = 0,
-  delta = 0.15,
-  parE = c(1, 1),
-  parS = c(50007530, 49924090)
+  x = c(16, 17),
+  n = 23,
+  parE = c(0.6, 0.4),
+  parS = c(0.6, 0.4),
+  delta = 0.1,
+  relativeDelta = FALSE
 )
 
-postprob(x = 27, n = 34, p = 0.65, parE = c(1, 1))
-## ok, almost
+# When we use a relative difference and look at several possible number of responses.
+postprobDist(
+  x = 1:23,
+  n = 23,
+  parE = c(0.2, 0.8),
+  parS = c(0.6, 0.4),
+  delta = 0.1,
+  relativeDelta = TRUE
+)
 
-## try out mixtures:
-## play around with the beta parameters and weights to
-## get a feeling.
-## Note that very extreme beta parameters do no longer increase
-## the return value, because then that mixture component is too
-## unlikely a posteriori
+# When we use beta mixtures for both the Experimental and SOC arms.
 postprobDist(
-  x = 16, n = 23,
+  x = 16,
+  n = 23,
   parE =
     rbind(
       c(0.6, 0.4),
@@ -31,5 +42,44 @@ postprobDist(
       c(0.6, 0.4),
       c(10, 10)
     ),
-  weightsS = c(1, 3)
+  weightsS = c(1, 3),
+  delta = 0.1
+)
+
+# Experimental arm only (strictly single arm trial), uniform prior in Experimental arm.
+# Non-uniform Prior used for SOC arm as no precedent data.
+postprobDist(
+  x = 16,
+  n = 23,
+  xS = 0,
+  nS = 0,
+  delta = 0,
+  relativeDelta = FALSE,
+  parS = c(2, 3),
+  weightsS = 1
+)
+# Experimental arm and SOC, uniform prior in both E and S arms, default setting used.
+postprobDist(
+  x = 16,
+  n = 20,
+  xS = 10,
+  nS = 20,
+  delta = 0,
+  relativeDelta = FALSE,
+  weightsS = 1
+)
+
+# Experimental and SOC arm, with beta mix prior for both arms.
+# For each of the SOC arm is of 3 priors, therefore 3 sets of beta parameters, and 3 weights.
+postprobDist(
+  x = 16,
+  n = 20,
+  xS = 10,
+  nS = 20,
+  delta = 0.1,
+  relativeDelta = TRUE,
+  parE = rbind(c(1, 1), c(3, 4), c(8, 9)),
+  weights = c(5, 3, 2),
+  parS = rbind(c(4, 5), c(2, 3), c(4, 4)),
+  weightsS = c(2, 5, 3)
 )
diff --git a/inst/WORDLIST b/inst/WORDLIST
index 223e850f..1db7dca6 100644
--- a/inst/WORDLIST
+++ b/inst/WORDLIST
@@ -74,6 +74,7 @@ dw
 dX
 dY
 dZ
+ee
 eq
 estunated
 ExpectedN
@@ -88,6 +89,7 @@ hyoithesis
 hyperprior
 increasement
 inefficacious
+Integrand
 integrations
 LciU
 ldots
@@ -105,6 +107,7 @@ nn
 nnE
 nnF
 nnr
+normalised
 nS
 oc
 ocPostProb
diff --git a/man/dbetabinom.Rd b/man/dbetabinom.Rd
index 1fb32cee..b59b273f 100644
--- a/man/dbetabinom.Rd
+++ b/man/dbetabinom.Rd
@@ -28,9 +28,6 @@ Calculates the density function of the beta-binomial distribution.
 The beta-binomial density function has the following form:
 \verb{p(x) = (m! / (x!*(m-x)!)) * Beta(x+a,m-x+b) / Beta(a,b)}
 }
-\note{
-\code{x}, \code{a} and \code{b} can be vectors.
-}
 \examples{
 dbetabinom(x = 2, m = 29, a = 0.2, b = 0.4, log = FALSE)
 
diff --git a/man/dbetabinomMix.Rd b/man/dbetabinomMix.Rd
index 6ec5b36e..67f3d116 100644
--- a/man/dbetabinomMix.Rd
+++ b/man/dbetabinomMix.Rd
@@ -26,9 +26,6 @@ The (log) density values of the mixture of beta-binomial distributions at \code{
 
 Calculates the density function for a mixture of beta-binomial distributions.
 }
-\note{
-\code{x} can be a vector.
-}
 \examples{
 dbetabinomMix(x = 2, m = 29, par = rbind(c(0.2, 0.4)), weights = 1)
 
diff --git a/man/h_get_bounds.Rd b/man/h_get_bounds.Rd
new file mode 100644
index 00000000..f1ffc1a4
--- /dev/null
+++ b/man/h_get_bounds.Rd
@@ -0,0 +1,20 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/postprobDist.R
+\name{h_get_bounds}
+\alias{h_get_bounds}
+\title{Generating bounds for the Integration of Beta Mixture Posterior}
+\usage{
+h_get_bounds(controlBetamixPost)
+}
+\arguments{
+\item{controlBetamixPost}{(\code{list}):\cr a list of posterior parameters of a beta-mixture-binomial distribution with generic names
+\code{par} and \code{weights}. See \verb{[getBetaMix()]}.}
+}
+\value{
+Integrand function
+}
+\description{
+Using the quantile of the Beta Mixture Distribution from parameters given by standard of care \code{SOC} or
+experimental group \code{E} to determine bounds as inputs to \verb{[stats::integrate()]}.
+}
+\keyword{internal}
diff --git a/man/h_integrand.Rd b/man/h_integrand.Rd
new file mode 100644
index 00000000..7b9c2ef1
--- /dev/null
+++ b/man/h_integrand.Rd
@@ -0,0 +1,28 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/postprobDist.R
+\name{h_integrand}
+\alias{h_integrand}
+\title{The Posterior Beta Mixture Integrand when Delta is absolute}
+\usage{
+h_integrand(p_s, delta, activeBetamixPost, controlBetamixPost)
+}
+\arguments{
+\item{p_s}{(\code{number}):\cr probability of success or response rate of standard of care or \code{SOC} group.}
+
+\item{delta}{(\code{number}):\cr margin by which the response rate in the treatment group should
+be better than in the SOC group. Must be >= \code{0}. See note.}
+
+\item{activeBetamixPost}{(\code{list}):\cr a list of posterior parameters of a beta-mixture-binomial distribution with generic names
+\code{par} and \code{weights}. See \verb{[getBetaMix()]}.}
+
+\item{controlBetamixPost}{(\code{list}):\cr a list of posterior parameters of a beta-mixture-binomial distribution with generic names
+\code{par} and \code{weights}. See \verb{[getBetaMix()]}.}
+}
+\value{
+Function that is an argument for \verb{[stats::integrate()]}.
+}
+\description{
+The helper function to generate Integrand function when \verb{relative Delta = FALSE},
+a default setting.
+}
+\keyword{internal}
diff --git a/man/h_integrand_relDelta.Rd b/man/h_integrand_relDelta.Rd
new file mode 100644
index 00000000..1cc98f90
--- /dev/null
+++ b/man/h_integrand_relDelta.Rd
@@ -0,0 +1,27 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/postprobDist.R
+\name{h_integrand_relDelta}
+\alias{h_integrand_relDelta}
+\title{The Posterior Beta Mixture Integrand when \code{delta} is relative}
+\usage{
+h_integrand_relDelta(p_s, delta, activeBetamixPost, controlBetamixPost)
+}
+\arguments{
+\item{p_s}{(\code{number}):\cr probability of success or response rate of standard of care or \code{SOC} group.}
+
+\item{delta}{(\code{number}):\cr margin by which the response rate in the treatment group should
+be better than in the SOC group. Must be >= \code{0}. See note.}
+
+\item{activeBetamixPost}{(\code{list}):\cr a list of posterior parameters of a beta-mixture-binomial distribution with generic names
+\code{par} and \code{weights}. See \verb{[getBetaMix()]}.}
+
+\item{controlBetamixPost}{(\code{list}):\cr a list of posterior parameters of a beta-mixture-binomial distribution with generic names
+\code{par} and \code{weights}. See \verb{[getBetaMix()]}.}
+}
+\value{
+Function that is an argument for \verb{[stats::integrate()]}.
+}
+\description{
+The helper function to generate Integrand function when \verb{relative Delta = TRUE}.
+}
+\keyword{internal}
diff --git a/man/postprobDist.Rd b/man/postprobDist.Rd
index 4658f604..348f3e74 100644
--- a/man/postprobDist.Rd
+++ b/man/postprobDist.Rd
@@ -2,14 +2,14 @@
 % Please edit documentation in R/postprobDist.R
 \name{postprobDist}
 \alias{postprobDist}
-\title{Compute the posterior probability with beta prior on SOC}
+\title{Compute the Posterior Probability with Beta Prior on \code{SOC}}
 \usage{
 postprobDist(
   x,
   n,
   xS = 0,
   nS = 0,
-  delta = 0,
+  delta,
   relativeDelta = FALSE,
   parE = c(1, 1),
   weights,
@@ -18,83 +18,103 @@ postprobDist(
 )
 }
 \arguments{
-\item{x}{number of successes (in the treatment group). Note that \code{x}
-can be a vector.}
+\item{x}{(\code{numeric}):\cr number of success counts in the treatment group.}
+
+\item{n}{(\code{number}):\cr number of patients in the treatment group.}
 
-\item{n}{number of patients (in the treatment group)}
+\item{xS}{(\code{number}):\cr number of success counts in the SOC group.}
 
-\item{xS}{number of successes in the SOC group (default: 0)}
+\item{nS}{(\code{number}):\cr number of patients in the SOC group.}
 
-\item{nS}{number of patients in the SOC group (default: 0)}
+\item{delta}{(\code{number}):\cr margin by which the response rate in the treatment group should
+be better than in the SOC group. Must be >= \code{0}. See note.}
 
-\item{delta}{margin by which the response rate in the treatment group should
-be better than in the SOC group (default: 0)}
+\item{relativeDelta}{(\code{flag}):\cr If \code{TRUE}, then a \code{relativeDelta} is used. Represents that a minimum
+response rate in magnitude of \code{delta} of the SOC non-responding patients. See note.}
 
-\item{relativeDelta}{should the delta be relative? (not default). If this is
-\code{TRUE}, then a relative delta is used. This means we want to have
-response at least in delta proportion of the SOC non-responding patients.
-Non-responding patients rate is 1 - P_S, and we want to have P_S + (1 - P_S)
-\itemize{
-\item delta response rate (at least) in the treatment. That is, we evaluate the
-posterior probability Pr(P_E > P_S + (1 - P_S) * delta | data).
-}}
+\item{parE}{(\code{numeric} or \code{matrix}):\cr parameters for beta distribution. If it is a matrix, it needs to have 2 columns,
+and each row corresponds to each component of a beta-mixture distribution
+for the \code{E} group. See details.}
 
-\item{parE}{the beta parameters matrix, with K rows and 2 columns,
-corresponding to the beta parameters of the K components. default is a
-uniform prior.}
+\item{weights}{(\code{numeric}):\cr the non-negative mixture weights of the beta mixture prior for group \code{E}. See details.}
 
-\item{weights}{the mixture weights of the beta mixture prior. Default are
-uniform weights across mixture components.}
+\item{parS}{(\code{numeric} or \code{matrix}):\cr parameters for beta distribution. If it is a matrix, it needs to have 2 columns,
+and each row corresponds to each component of a beta-mixture distribution
+for the \code{S} group. See details.}
 
-\item{parS}{beta parameters for the SOC group (default: uniform)}
+\item{weightsS}{(\code{numeric}):\cr the non-negative mixture weights of the beta mixture prior for group \code{S}. See details.}
 
-\item{weightsS}{weights for the SOC group (default: uniform)}
+\item{epsilon}{(\code{number}):\cr the smallest non-negative floating number to represent the lower bound for
+the interval of integration.}
 }
 \value{
-the posterior probability
+The posterior probability
 }
 \description{
-Using the approach by Thall and Simon (Biometrics, 1994), evaluate the
-posterior probability of having Pr(P_E > P_S + delta | data) (but see below
-for relative delta margin). Both for the new treatment E as well as for the
-SOC S data might be available. However the default is that no data is
-available for the SOC, corresponding to the single arm trial situation. Note
-that a uniform prior is the useful default for the treatment proportion,
-while in the single arm trial an informative prior on the SOC proportion is
-useful.
+\ifelse{html}{\href{https://lifecycle.r-lib.org/articles/stages.html#experimental}{\figure{lifecycle-experimental.svg}{options: alt='[Experimental]'}}}{\strong{[Experimental]}}
+
+Using the approach by Thall and Simon (Biometrics, 1994), this evaluates the
+posterior probability of achieving superior response rate in the treatment group compared to standard of care (SOC).
+See note below for two formulations of the difference in response rates.
 }
 \details{
-Beta mixture prior can be specified for the treatment (\code{parE}
-and \code{weights} parameters) and control proportion (\code{parS} and
-\code{weightsS} parameters), see \code{\link{postprob}} for details. Note
-that being able to specify a beta mixture prior also on the control
-treatment is e.g. important for the futility decision making (see the
-\code{\link{oc2}} code).
+The beta mixture prior for the \code{E} arm requires argument \code{parE} and \code{weights}.
+The beta mixture prior for the \code{S} arm requires argument \code{parS} and \code{weightsS}.
+See \verb{[postprob()]} for details.
+
+If a beta-mixture is used, by default, the weights are uniform across the components.
+Weights can exceed 1, to which the algorithm will normalize the weights such that all weights sum to 1.
+}
+\note{
+\subsection{Delta :}{
+
+The desired improvement is denoted as \code{delta}. There are two options in using \code{delta}.
+The absolute case when \code{relativeDelta = FALSE} and relative as when \code{relativeDelta = TRUE}.
+\enumerate{
+\item The absolute case is when we define an absolute delta, greater than \code{P_S},
+the response rate of the \code{SOC} group such that
+the posterior is \code{Pr(P_E > P_S + delta | data)}.
+\item In the relative case, we suppose that the treatment group's
+response rate is assumed to be greater than \code{P_S + (1-P_S) * delta} such that
+the posterior is \code{Pr(P_E > P_S + (1 - P_S) * delta | data)}.
+}
+}
 }
 \examples{
-## example similar to Lee and Liu:
-postprobDist(x = 16, n = 23, parE = c(0.6, 0.4), parS = c(0.6, 0.4))
+# An example similar to Lee and Liu (2008).
+postprobDist(
+  x = 16,
+  n = 23,
+  parE = c(0.6, 0.4),
+  parS = c(0.6, 0.4),
+  delta = 0.1,
+  relativeDelta = FALSE
+)
 
-## these two should give the same result:
+# For a sequence of success outcomes for Experimental arm.
 postprobDist(
-  x = 27, n = 34,
-  xS = 0, nS = 0,
-  delta = 0.15,
-  parE = c(1, 1),
-  parS = c(50007530, 49924090)
+  x = c(16, 17),
+  n = 23,
+  parE = c(0.6, 0.4),
+  parS = c(0.6, 0.4),
+  delta = 0.1,
+  relativeDelta = FALSE
 )
 
-postprob(x = 27, n = 34, p = 0.65, parE = c(1, 1))
-## ok, almost
+# When we use a relative difference and look at several possible number of responses.
+postprobDist(
+  x = 1:23,
+  n = 23,
+  parE = c(0.2, 0.8),
+  parS = c(0.6, 0.4),
+  delta = 0.1,
+  relativeDelta = TRUE
+)
 
-## try out mixtures:
-## play around with the beta parameters and weights to
-## get a feeling.
-## Note that very extreme beta parameters do no longer increase
-## the return value, because then that mixture component is too
-## unlikely a posteriori
+# When we use beta mixtures for both the Experimental and SOC arms.
 postprobDist(
-  x = 16, n = 23,
+  x = 16,
+  n = 23,
   parE =
     rbind(
       c(0.6, 0.4),
@@ -105,6 +125,45 @@ postprobDist(
       c(0.6, 0.4),
       c(10, 10)
     ),
-  weightsS = c(1, 3)
+  weightsS = c(1, 3),
+  delta = 0.1
+)
+
+# Experimental arm only (strictly single arm trial), uniform prior in Experimental arm.
+# Non-uniform Prior used for SOC arm as no precedent data.
+postprobDist(
+  x = 16,
+  n = 23,
+  xS = 0,
+  nS = 0,
+  delta = 0,
+  relativeDelta = FALSE,
+  parS = c(2, 3),
+  weightsS = 1
+)
+# Experimental arm and SOC, uniform prior in both E and S arms, default setting used.
+postprobDist(
+  x = 16,
+  n = 20,
+  xS = 10,
+  nS = 20,
+  delta = 0,
+  relativeDelta = FALSE,
+  weightsS = 1
+)
+
+# Experimental and SOC arm, with beta mix prior for both arms.
+# For each of the SOC arm is of 3 priors, therefore 3 sets of beta parameters, and 3 weights.
+postprobDist(
+  x = 16,
+  n = 20,
+  xS = 10,
+  nS = 20,
+  delta = 0.1,
+  relativeDelta = TRUE,
+  parE = rbind(c(1, 1), c(3, 4), c(8, 9)),
+  weights = c(5, 3, 2),
+  parS = rbind(c(4, 5), c(2, 3), c(4, 4)),
+  weightsS = c(2, 5, 3)
 )
 }
diff --git a/tests/testthat/test-postprob.R b/tests/testthat/test-postprob.R
index e9381886..38d8732c 100644
--- a/tests/testthat/test-postprob.R
+++ b/tests/testthat/test-postprob.R
@@ -19,7 +19,7 @@ test_that("postprob gives the correct number result", {
 })
 
 test_that("postprob gives the correct number result", {
-  # 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
+  # 2 component beta mixture prior, i.e., P_E ~ 0.5*beta(0.6,0.4) + 0.5*beta(1,1) and Pr(P_E > p | data) = 0.05559802
   result <- postprob(
     x = 10,
     n = 23,
diff --git a/tests/testthat/test-postprobDist.R b/tests/testthat/test-postprobDist.R
new file mode 100644
index 00000000..15cae106
--- /dev/null
+++ b/tests/testthat/test-postprobDist.R
@@ -0,0 +1,208 @@
+# postprobDist ----
+
+test_that("postprobDist gives the correct number result.", {
+  result <- postprobDist(x = 16, n = 23, parE = c(0.6, 0.4), parS = c(0.6, 0.4), delta = 0.1)
+  expect_equal(result, 0.4431067, tolerance = 1e-5)
+})
+
+test_that("postprobDist gives incrementally higher values with larger x", {
+  is_lower <- postprobDist(x = 16, n = 23, delta = 0.1, parE = c(0.6, 0.4), parS = c(0.6, 0.4))
+  is_higher <- postprobDist(x = 20, n = 23, delta = 0.1, parE = c(0.6, 0.4), parS = c(0.6, 0.4))
+  expect_true(is_lower < is_higher)
+})
+
+test_that("postprobDist gives higher values with larger x and returns identical result when x is a vector", {
+  expected_lower <- postprobDist(x = 16, n = 23, delta = 0.1, parE = c(0.6, 0.4), parS = c(0.6, 0.4))
+  expected_higher <- postprobDist(x = 20, n = 23, delta = 0.1, parE = c(0.6, 0.4), parS = c(0.6, 0.4))
+  result <- postprobDist(x = c(16, 20), n = 23, delta = 0.1, parE = c(0.6, 0.4), parS = c(0.6, 0.4))
+  expect_identical(result, c(expected_lower, expected_higher))
+})
+
+test_that("postprobDist gives the correct result for a beta-mixture", {
+  result <- postprobDist(
+    x = 10,
+    n = 23,
+    delta = 0.1,
+    parE = rbind(
+      c(0.6, 0.4),
+      c(1, 1)
+    ),
+    parS = rbind(
+      c(0.6, 0.4),
+      c(1, 1)
+    )
+  )
+  expect_equal(result, 0.3143941, tolerance = 1e-5)
+})
+
+test_that("postprobDist gives incrementally higher values with larger x for a beta-mixture", {
+  is_lower <- postprobDist(
+    x = 10,
+    n = 23,
+    delta = 0.1,
+    parE = rbind(
+      c(0.6, 0.4),
+      c(1, 1)
+    ),
+    parS = rbind(
+      c(0.6, 0.4),
+      c(1, 1)
+    )
+  )
+  is_higher <- postprobDist(
+    x = 16,
+    n = 23,
+    delta = 0.1,
+    parE = rbind(
+      c(0.6, 0.4),
+      c(1, 1)
+    ),
+    parS = rbind(
+      c(0.6, 0.4),
+      c(1, 1)
+    )
+  )
+  expect_true(is_lower < is_higher)
+})
+
+test_that("postprobDist gives the correct result for a weighted beta-mixture", {
+  result <- postprobDist(
+    x = 10,
+    n = 23,
+    delta = 0.1,
+    parE = rbind(
+      c(0.6, 0.4),
+      c(1, 1)
+    ),
+    parS = rbind(
+      c(0.6, 0.4),
+      c(1, 1)
+    ),
+    weights = c(0.5, 0.5),
+    weightsS = c(0.3, 0.7)
+  )
+  expect_equal(result, 0.3248885, tolerance = 1e-4)
+})
+
+test_that("postprobDist gives an error when n is not a number.", {
+  expect_error(
+    results <- postprobDist(
+      x = c(16, 17),
+      n = c(23, 20),
+      parE = c(0.6, 0.4),
+      parS = c(0.6, 0.4),
+      delta = 0.1,
+      relativeDelta = FALSE
+    ), "number of items to replace is not a multiple of replacement length"
+  )
+})
+
+test_that("postprobDist gives an error when xS and nS are not numbers", {
+  expect_error(
+    results <- postprobDist(
+      x = c(10, 16),
+      n = 23,
+      xS = c(9, 10),
+      nS = c(20, 22),
+      delta = 0.1,
+      parE = c(0.6, 0.4),
+      parS = c(0.6, 0.4),
+      weights = c(0.5),
+      weightsS = c(0.3),
+    ), "Must have length 1"
+  )
+})
+
+# h_integrand_relDelta ----
+test_that("h_integrand_relDelta gives the correct numerical result for a beta-mixture.", {
+  x <- 16
+  n <- 23
+  xS <- 10
+  nS <- 20
+  parE <- t(c(1, 3))
+  parS <- t(c(1, 1))
+  weights <- 1
+  weightsS <- 1
+  p_s <- 0.1
+  delta <- 0.1
+  relativeDelta <- TRUE
+  activeBetamixPost <- getBetamixPost(x = x, n = n, par = parE, weights = weights)
+  controlBetamixPost <- getBetamixPost(x = xS, n = nS, par = parS, weights = weightsS)
+  results <- h_integrand_relDelta(
+    p_s = p_s,
+    delta = delta,
+    activeBetamixPost = activeBetamixPost,
+    controlBetamixPost = controlBetamixPost
+  )
+  expect_equal(results, 0.0001352829, tolerance = 1e-4)
+})
+
+test_that("h_integrand_relDelta gives the correct numerical result for a weighted beta-mixture.", {
+  x <- 16
+  n <- 23
+  xS <- 10
+  nS <- 20
+  parE <- rbind(c(1, 3), c(2, 3))
+  parS <- rbind(c(1, 1), c(3, 4))
+  weights <- c(5, 10)
+  weightsS <- c(3, 4)
+  p_s <- 0.1
+  delta <- 0.1
+  relativeDelta <- TRUE
+  activeBetamixPost <- getBetamixPost(x = x, n = n, par = parE, weights = weights)
+  controlBetamixPost <- getBetamixPost(x = xS, n = nS, par = parS, weights = weightsS)
+  results <- h_integrand_relDelta(
+    p_s = p_s,
+    delta = delta,
+    activeBetamixPost = activeBetamixPost,
+    controlBetamixPost = controlBetamixPost
+  )
+  expect_equal(results, 6.498862e-05, tolerance = 1e-4)
+})
+
+# h_integrand ----
+test_that("h_integrand gives the correct numerical result for a beta-mixture", {
+  x <- 16
+  n <- 23
+  xS <- 10
+  nS <- 20
+  parE <- t(c(1, 3))
+  parS <- t(c(1, 1))
+  weights <- 1
+  weightsS <- 1
+  p_s <- 0.1
+  delta <- 0.1
+  relativeDelta <- TRUE
+  activeBetamixPost <- getBetamixPost(x = x, n = n, par = parE, weights = weights)
+  controlBetamixPost <- getBetamixPost(x = xS, n = nS, par = parS, weights = weightsS)
+  results <- h_integrand(
+    p_s = p_s,
+    delta = delta,
+    activeBetamixPost = activeBetamixPost,
+    controlBetamixPost = controlBetamixPost
+  )
+  expect_equal(results, 0.0001352828, tolerance = 1e-4)
+})
+
+test_that("h_integrand works as expected for a weighted beta-mixture.", {
+  x <- 16
+  n <- 23
+  xS <- 10
+  nS <- 20
+  parE <- rbind(c(1, 3), c(2, 3))
+  parS <- rbind(c(1, 1), c(3, 4))
+  weights <- c(5, 10)
+  weightsS <- c(3, 4)
+  p_s <- 0.1
+  delta <- 0.1
+  relativeDelta <- FALSE
+  activeBetamixPost <- getBetamixPost(x = x, n = n, par = parE, weights = weights)
+  controlBetamixPost <- getBetamixPost(x = xS, n = nS, par = parS, weights = weightsS)
+  results <- h_integrand(
+    p_s = p_s,
+    delta = delta,
+    activeBetamixPost = activeBetamixPost,
+    controlBetamixPost = controlBetamixPost
+  )
+  expect_equal(results, 6.498861e-05, tolerance = 1e-4)
+})