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

Finally tore myself away from videogaming to take a look at this. Thank you for such helpful commenting, really appreciate it.

So glad I posted the beta of it (and am finding errors before submission). It's used in one line of code in a behemoth of a simulation algorithm I'm building in another package, and my collaborator mentioned he'd make use of it, too, so I thought I'd spin it off for him to use. Did not expect so much interest! I thought It would have an audience of n = 1.

So nice of you to put together such a clear solution for me :)

@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