2
$\begingroup$

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 β?

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

2 Answers 2

4
$\begingroup$

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
$\endgroup$
3
$\begingroup$

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
$\endgroup$

Not the answer you're looking for? Browse other questions tagged or ask your own question.