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 \code{\link{get_model_input_1exp}} and \code{\link{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) {
  # check inputs
  assert_that(is.data.frame(df))
  assert_that(all(names(df) == c("w", "cum_probs")))
  assert_that(is.numeric(df$w))
  assert_that(is.numeric(df$cum_probs))
  # beta fit function
  w <- df$w
  cum_probs <- df$cum_probs
  ## exclude cum_probs that are 0 or 1
  inc <- (cum_probs > 0) & (cum_probs < 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 <- qnorm(min_cum_probs)
  max_q <- 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
  ## set starting values
  alpha <- abs(m^3/v * (1/m - 1) - m)
  beta <- abs(alpha/m - alpha)
  if(identical(cum_probs[inc], w[inc])){
    alpha <- beta <- 1
  }
  ## function to be optimized
  beta_error <- function(params, weight, probs) {
    sum((pbeta(q = weight, shape1 = exp(params[1]), shape2 = exp(params[2])) -
           probs)^2)
  }
  ## optimization
  beta_fit <- 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")
  return(beta_fit)
}