Skip to content

Instantly share code, notes, and snippets.

@daob
Last active August 16, 2020 09:33
Show Gist options
  • Star 0 You must be signed in to star a gist
  • Fork 0 You must be signed in to fork a gist
  • Save daob/1422e978ff98bdf466fbcb4d9bf3e53e to your computer and use it in GitHub Desktop.
Save daob/1422e978ff98bdf466fbcb4d9bf3e53e to your computer and use it in GitHub Desktop.
# 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])
@softloud
Copy link

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

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment