Last active
August 16, 2020 09:33
-
-
Save daob/1422e978ff98bdf466fbcb4d9bf3e53e to your computer and use it in GitHub Desktop.
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
# 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.
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
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 :)