# How can I calculate the α and β parameters for a Beta distribution if I know the mean and 95% upper boundary I want the distribution to have? [duplicate]

As an input I have the mean (expected value) of a probability and the 95% upper bound. Now, I would like to define the α and β parameters for a Beta distribution such that I can do Monte Carlo simulations. I know the expressions for these parameters in terms of E[X] and Var[X], but do not have Var[X] directly available. Is it possible to find Var[X] with the given upper bound, or are there other ways to find α and β?

• My first choice would be numerical optimization to find the best fitting distribution. Commented Oct 14, 2022 at 9:13

Here's an R solution using a root solver that exploits the relationship between the mean and the parameters. Because $$\bar{x}=\frac{\alpha}{\alpha + \beta}$$, once we know $$\alpha$$, we also know $$\beta$$ through $$\beta = \frac{\alpha - \bar{x}\alpha}{\bar{x}}$$.

# Setting the true parameters for illustration
alpha <- 2
beta <- 46
q <- 0.975 # The quantile

# The mean and quantile
bmean <- alpha/(alpha + beta)
qval <- qbeta(q, alpha, beta)

# Function to find alpha using root solver
find_alpha <- function(alpha, mean, qval, q) {
beta <- (alpha - alpha*mean)/mean # Beta can be calculated from the known mean and the estimated alpha
qbeta(q, alpha, beta) - qval
}

# Estimate alpha using a root solver
alphaest <- uniroot(find_alpha, interval = c(0, 1000), mean = bmean, qval = qval, q = q, extendInt = "upX")$root # Calculate beta exploiting the relationship between alpha and the known mean betaest <- (alphaest - alphaest*bmean)/bmean # The final estimates c(alphaest, betaest) [1] 2.00000 46.00001  Using optimisation: target_mean = .3 target_upper_ci = .7 try_parameter = function(a, return_info = F){ # Since dist_mean = a / (a + b), find value of b that achieves desired mean b = (a - a*target_mean) / target_mean # Find 95% upper boundary for these parameters upper_ci = qbeta(.975, a, b) if(return_info){ data.frame(a, b, mean = a / (a + b), upper_ci) } else { # Return error for optimisation abs(target_upper_ci - upper_ci) } } a_values = seq(0, 10, .01) plot(a_values, try_parameter(a_values), 'l')  result = optimise(try_parameter, interval = c(.1, 10)) result #> $$minimum #> [1] 1.638568 #> #>$$objective #> [1] 7.156458e-07 try_parameter(result$minimum, return_info = T)
#>          a        b mean  upper_ci
#> 1 1.638568 3.823326  0.3 0.7000007