R/fit_beta_1exp.R

In tipmap: Tipping Point Analysis for Bayesian Dynamic Borrowing

#' @title
#' Fit beta distribution for one expert
#'
#' @description
#' Fit beta distribution to data elicited from one expert via the roulette method.
#'
#' @param df A dataframe generated by \code{\link{get_model_input_1exp}}.
#'
#' @details This function is based on \code{SHELF::fitdist} and yields identical results.
#'
#' @return Parameters (\code{alpha} and \code{beta}) of a beta fit.
#' 
#' @export
#' 
#' @seealso [get_model_input_1exp()] and [fit_beta_mult_exp()]. 
#' 
#' @examples
#' chips <- c(0, 2, 3, 2, 1, 1, 1, 0, 0, 0)
#' x <- get_cum_probs_1exp(chips)
#' y <- get_model_input_1exp(x)
#' fit_beta_1exp(df = y)["par"] 
#'
fit_beta_1exp <- function(df) {
  assert_that(is.data.frame(df), msg = "`df` must be a data frame")
  assert_that(all(names(df) == c("w", "cum_probs")), msg = "`df` must contain exactly `w` and `cum_probs`")
  assert_that(is.numeric(df$w), msg = "`df$w` must be numeric")
  assert_that(is.numeric(df$cum_probs), msg = "`df$cum_probs` must be numeric")
  assert_that(all(is.finite(df$w)), msg = "`df$w` must be finite")
  assert_that(all(is.finite(df$cum_probs)), msg = "`df$cum_probs` must be finite")
  assert_that(all(df$cum_probs >= 0 & df$cum_probs <= 1), msg = "`df$cum_probs` must lie in [0, 1]")
  assert_that(all(diff(df$cum_probs) >= 0), msg = "`df$cum_probs` must be non-decreasing")
  
  w <- df$w
  cum_probs <- df$cum_probs
  inc <- (cum_probs > 0) & (cum_probs < 1)
  
  assert_that(sum(inc) >= 2, msg = "Need at least two cumulative probabilities strictly between 0 and 1")
  
  min_cum_probs <- min(cum_probs[inc])
  max_cum_probs <- max(cum_probs[inc])
  min_weight <- min(w[inc])
  max_weight <- max(w[inc])
  
  min_q <- stats::qnorm(min_cum_probs)
  max_q <- stats::qnorm(max_cum_probs)
  
  m <- (min_weight * max_q - max_weight * min_q) / (max_q - min_q)
  v <- ((max_weight - min_weight) / (max_q - min_q))^2
  
  alpha <- abs(m^3 / v * (1 / m - 1) - m)
  beta <- abs(alpha / m - alpha)
  
  if (identical(cum_probs[inc], w[inc])) {
    alpha <- 1
    beta <- 1
  }
  
  assert_that(is.finite(alpha), msg = "Failed to derive a finite starting value for `alpha`")
  assert_that(is.finite(beta), msg = "Failed to derive a finite starting value for `beta`")
  assert_that(alpha > 0, msg = "Starting value for `alpha` must be positive")
  assert_that(beta > 0, msg = "Starting value for `beta` must be positive")
  
  beta_error <- function(params, weight, probs) {
    sum((stats::pbeta(q = weight, shape1 = exp(params[1]), shape2 = exp(params[2])) - probs)^2)
  }
  
  beta_fit <- stats::optim(
    par = c(log(alpha), log(beta)),
    fn = beta_error,
    weight = w[inc],
    probs = cum_probs[inc]
  )
  
  if (beta_fit$convergence != 0) {
    warning("Algorithm did not converge.")
  }
  
  beta_fit$par <- exp(beta_fit$par)
  names(beta_fit$par) <- c("alpha", "beta")
  beta_fit
}