R/pois.R

In fcci: Feldman-Cousins Confidence Intervals

#' Poisson rate confidence intervals
#' @md
#' @description
#' Build confidence intervals for the rate of a Poisson random variable using
#' the Feldman-Cousins construction (Feldman and Cousins, 1998).
#' @param n number of events. Numeric vector of length one.
#' @param b expected number of background events. Numeric vector of length one
#' or three (see Details).
#' @param cl confidence level for the returned confidence interval.
#' A number between zero and one.
#' @param acc required accuracy in the confidence interval endpoints.
#' A positive number.
#' @param lambda_min,lambda_max window of Poisson rate
#' values for whom the Neyman belt construction is performed. Numeric vectors of
#' length one or \code{NULL} (see Details). For advanced users only, use the
#' default value (\code{NULL}) if unsure.
#' @return a numeric vector of length two, containing the confidence interval
#' endpoints.
#' @details The Feldman-Cousins confidence intervals construction
#' (Feldman and Cousins, 1998) provides a unified treatment for the estimation
#' of Poisson rates, which
#' produces consistent results for both the common case, where the number of
#' observed events \code{n} is greater than the number of background events
#' \code{b}, and the boundary cases, when \code{n} is equal or even less than
#' \code{b}. In these last cases, a naive treatment can lead to either
#' overcovering (that is, unnecessarily large confidence intervals;
#' this is the case for, *e.g.*, the intervals returned from
#' \code{stats::poisson.test} for zero events and no background) or,
#' even worse, undercovering and/or non-sensical confidence intervals
#' (relevant examples can be found in (Feldman and Cousins, 1998)).
#'
#' In order to find confidence intervals, the standard Neyman construction
#' is performed over a grid of values for the estimated rate. The search
#' interval can be controlled through the parameters \code{lambda_min} and
#' \code{lambda_max}; if these are left as default (\code{NULL}), they are
#' estimated automatically based on a rough estimate, which should be adequate
#' in a large majority of cases.
#'
#' The argument \code{b} controls the expected number of background events,
#' which would differ from zero if some of the event counts are expected
#' to correspond to spurious events. In other words, the number of total counts
#' \code{n} is assumed to follow a Poisson distribution with rate
#' \code{lambda + b}, with \code{b} known and \code{lambda} unknown and to be
#' estimated.
#'
#' The original reference discusses a minor correction to confidence intervals,
#' which ensures that the produced upper limits are decreasing functions of
#' \code{b}.
#' This is computationally intensive and not necessary for the correctness of
#' Feldman-Cousins confidence intervals, but may facilitate interpretation for
#' e.g. sensitivity analyses.
#' The correction can be implemented by providing a length three numeric to the
#' \code{b} argument, such that \code{b[[1]]} is the expected background rate as
#' before, while \code{b[[2]]} and \code{b[[3]]} represent a maximum value
#' for \code{b} and a grid step. The construction is repeated over the whole
#' grid of \code{b} values and the actual upper limit is chosen as the
#' greatest upper limit encountered in the sequence of intervals so constructed.
#' @references
#' Feldman, Gary J. and Cousins, Robert D.
#' "Unified approach to the classical statistical analysis of small signals"
#' *Phys. Rev. D* **57**, issue 7 (1998): 3873-3889.
#' @examples
#' confint_pois(n = 10, cl = 0.68)
#' confint_pois(n = 10, b = 5, cl = 0.68)
#'
#' # Compare:
#' confint_pois(n = 0, cl = 0.95)
#' # with:
#' poisson.test(x = 0, conf.level = 0.95)$conf.int
#' @export
confint_pois <- function(
	n, b = 0, cl = 0.95, acc = 1e-3 * (max(n - b, 0) + 1),
	lambda_min = NULL, lambda_max = NULL
	)
{
	check_args_pois(n, b, cl, acc, lambda_min, lambda_max)

	grid <- pois_lambda_grid(n, b, cl, acc, lambda_min, lambda_max)

	if (length(b) == 3)
		res <- confint_pois_adj_cpp(
			n, b[[1]], cl,
			grid[["min"]], grid[["max"]], grid[["step"]],
			b[[2]], b[[3]]
			)
	else
		res <- confint_pois_cpp(
			n, b, cl, grid[["min"]], grid[["max"]], grid[["step"]]
			)

	check_truncation_pois(res, grid, acc)

	attr(res, "cl") <- cl

	return(res)
}


#-------------------------------- Helpers -------------------------------------#
check_args_pois <- function(n, b, cl, acc, lambda_min, lambda_max)
{
	tryCatch(
		# try
		{
		assertthat::assert_that(
			is_event_count(n),
			is.numeric(b), length(b) %in% c(1, 3),
			is_non_negative(b[[1]]),
			length(b) == 1 || b[[2]] > b[[1]],
			length(b) == 1 || is_positive(b[[3]]),
			is_probability(cl), cl > 0, cl < 1,
			is_positive(acc)
			)
		if (!is.null(lambda_min))
			assertthat::assert_that(is_non_negative(lambda_min))
		if (!is.null(lambda_max))
			assertthat::assert_that(is_non_negative(lambda_max))
		if (!is.null(lambda_min) && !is.null(lambda_max))
			assertthat::assert_that(lambda_min < lambda_max)
		}
		,
		# catch
		error = function(cnd)
			rlang::abort(cnd$message, class = "domain_error")
		)
}

pois_lambda_grid <- function(n, b, cl, acc, lambda_min, lambda_max)
{
	if (is.null(lambda_min) || is.null(lambda_max)) {
		z <- -qnorm((1 - cl) / 2)
		hw <- 2 * z * sqrt(n + 10)
	}

	if (is.null(lambda_min))
		lambda_min <- max(0, n - b - hw)

	if (is.null(lambda_max))
		lambda_max <- max(n - b, 0) + hw

	lambda_step <- min(acc, lambda_max - lambda_min) / 2

	return( list(min = lambda_min, max = lambda_max, step = lambda_step) )
}

check_truncation_pois <- function(res, grid, acc) {
	if (grid[["min"]] + acc > res[[1]] & res[[1]] > 0) {
		h <- "Lower endpoint truncation"
		x <- "Lower limit might be overestimated due to grid size."
		i <- paste0("Try reducing 'lambda_min' or improving 'acc',",
			    "until this warning stops to appear.")
		rlang::warn(c(h, x = x, i = i), class = "truncation_warning")
	}

	if (res[[2]] + acc > grid[["max"]]) {
		h <- "Upper endpoint truncation"
		x <- "Upper limit might be underestimated due to grid size."
		i <- paste0("Try increasing 'lambda_max' or improving 'acc',",
			    "until this warning stops to appear.")
		rlang::warn(c(h, x = x, i = i), class = "truncation_warning")
	}
}