beta.R

# Function that takes desrired mean, distance, and probability, and outputs
#   another function to be optimized.
get_objective <- function(desired_mean, desired_dist, desired_mass) {
  # Return a function that can be passed to optim()
  function(shape1) {
    # Enforce desired mean:
    shape2 <- shape1 * ((1 / desired_mean) - 1)
    
    # Use R internals to get the definite integral:
    current_mass <- pbeta(q = desired_mean + desired_dist, 
                          shape1 = shape1, shape2 = shape2) - 
                    pbeta(q = desired_mean - desired_dist, 
                          shape1 = shape1, shape2 = shape2)
    
    # Loss is squared difference between desired and obtained measure:
    (current_mass - desired_mass)^2 
  }
}

# Get objective and optimize it to get params 
#  (not really necessary with Beta, but useful for generality)
est_parameters <- function(desired_mean, desired_dist, desired_mass) {
  my_objective <- get_objective(desired_mean = desired_mean, desired_dist = desired_dist, desired_mass = desired_mass)
  
  res <- optim(1, my_objective, method = "Brent", lower = 1e-3, upper = 10)
  
  shape1_est <- res$par
  shape2_est <- shape1_est * ((1/desired_mean) - 1)
  
  c(shape1_est, shape2_est)
}

# Get the parameters for Beta corresponding to mean, distance and mass
(pars <- est_parameters(0.3, 0.2, 0.8))

# Check that all is well
mean(rbeta(1e5, pars[1], pars[2])) # Mean is 0.3
# 80% of probability is between 0.1 and 0.5:
pbeta(q = 0.5, pars[1], pars[2]) - pbeta(q = 0.1, pars[1], pars[2])

Please don't forget to pr to add yourself to the contributors.