#' @return If `post` is identical to `TRUE`, a list containing the following
#' components is returned:
#' * `supp` - the support of the posterior distribution of R0
-#' * `pmf` - the probability mass function of the posterior distribution
+#' * `pmf` - the probability mass function of the posterior distribution of R0
#'
#' Otherwise, if `post` is identical to `FALSE`, only the estimate of R0 is
#' returned. Note that the estimate is equal to `sum(supp * pmf)` (i.e., the
#' posterior mean).
#'
-#' @references
-#' [Bettencourt and Riberio (PloS One, 2008)](
-#' https://doi.org/10.1371/journal.pone.0002185)
+#' @references [Bettencourt and Riberio (PloS One, 2008)](
+#' https://doi.org/10.1371/journal.pone.0002185)
+#'
+#' @seealso `vignette("seq_bayes_post", package = "Rnaught")` for examples of
+#' using the posterior distribution.
#'
#' @export
#'
#' # believed to be at most 4.
#' estimate <- seq_bayes(cases, mu = 1, kappa = 4)
#'
-#' # Same as above, but return the posterior distribution instead of the
+#' # Same as above, but return the posterior distribution of R0 instead of the
#' # estimate.
#' posterior <- seq_bayes(cases, mu = 1, kappa = 4, post = TRUE)
#'
+#' # Display the support and probability mass function of the posterior.
+#' posterior$supp
+#' posterior$pmf
+#'
#' # Note that the following always holds:
#' estimate == sum(posterior$supp * posterior$pmf)
seq_bayes <- function(cases, mu, kappa = 20, post = FALSE) {
+ validate_cases(cases, min_length = 2, min_count = 0)
+ if (!is_real(mu) || mu <= 0) {
+ stop("The serial interval (`mu`) must be a number greater than 0.",
+ call. = FALSE
+ )
+ }
+ if (!is_real(kappa) || kappa < 1) {
+ stop(
+ paste("The largest value of the uniform prior (`kappa`)",
+ "must be a number greater than or equal to 1."
+ ), call. = FALSE
+ )
+ }
+ if (!identical(post, TRUE) && !identical(post, FALSE)) {
+ stop("The posterior flag (`post`) must be set to `TRUE` or `FALSE`.",
+ call. = FALSE
+ )
+ }
+
if (any(cases == 0)) {
times <- which(cases > 0)
if (length(times) < 2) {
- stop("Vector of case counts must contain at least two positive integers.")
+ stop("Case counts must contain at least two positive integers.",
+ call. = FALSE
+ )
}
cases <- cases[times]
} else {