X-Git-Url: https://git.nmode.ca/Rnaught/blobdiff_plain/59de2553220ffff0a62ed3e454876df2bb93916b..HEAD:/R/seq_bayes.R diff --git a/R/seq_bayes.R b/R/seq_bayes.R index 2f87121..b25d273 100644 --- a/R/seq_bayes.R +++ b/R/seq_bayes.R @@ -45,15 +45,14 @@ #' @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. @@ -74,22 +73,43 @@ #' # 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) { - if (any(cases == 0)) { - times <- which(cases > 0) - if (length(times) < 2) { - stop("Vector of case counts must contain at least two positive integers.") - } - cases <- cases[times] - } else { - times <- seq_along(cases) + 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 + ) + } + + times <- which(cases > 0) + if (length(times) < 2) { + stop("Case counts must contain at least two positive integers.", + call. = FALSE + ) } + cases <- cases[times] support <- seq(0, kappa, 0.01) tau <- diff(times)