R/fit_beta_1exp.R

In tipmap: Tipping Point Analysis for Bayesian Dynamic Borrowing

#' Fit beta distribution for one expert
#'
#' @description
#' Fit beta distribution to weights elicited from one expert via the roulette method.
#'
#' @param df A dataframe generated by \code{get_model_input_1exp}.
#'
#' @details This function is based on \code{SHELF::fitdist} and yields identical results.
#'
#' @return Parameters of beta fit.
#'
#' @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)
{
  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)
}