Solving for shape1 and shape2 in qbeta() given its quantiles in R

Question

Background:

qbeta(p, shape1, shape2) is an built-in base R function. In this function, shape1 and shape2 are always > 0 (for my purposes always > 1). Also, p is probability thus 0 <= p <= 1.

Question:

Suppose I know qbeta(p = c(.025, .975), shape1 = x, shape2 = y ) = c(.3, .8 ), THEN, is there a way I can find out what x and y could reasonably be (i.e., within some margin of error)?

Note: There are, I believe, some relevant comments and codes provided in HERE. So, in case one of the colleagues can makes sense or use of those codes and comments and apply them here, I'll deeply appreciate it.

Finally, my ultimate goal is to wrap a function out of this x & y finding process.

Answer 1

Here is a brute force method. I'm hoping someone else can come up with a smarter way to do this.

d = expand.grid(shape1 = 1:100/10, shape2 = 1:100/10)

for (i in 1:NROW(d)){
    x = round(qbeta(p = c(0.025, 0.975), shape1 = d[i,1], shape2 = d[i,2]),2)
    if(x[1] > 0.28 & x[1] < 0.32 & x[2] > 0.78 & x[2] < 0.82){
        print(d[i,])
        break
    }
}

#     shape1 shape2
#5368    6.8    5.4

qbeta(p = c(.025, .975), shape1 = 6.8, shape2 = 5.4)
#[1] 0.2860268 0.8109389

I am providing whuber's solution from here below

################################
#REQUIRED FUNCTIONS
# Logistic transformation of the Beta CDF.
f.beta <- function(alpha, beta, x, lower=0, upper=1) {
    p <- pbeta((x-lower)/(upper-lower), alpha, beta)
    log(p/(1-p))
}

# Sums of squares.
delta <- function(fit, actual){
    sum((fit-actual)^2)
}

# The objective function handles the transformed parameters `theta` and
# uses `f.beta` and `delta` to fit the values and measure their discrepancies.
objective <- function(theta, x, prob, ...) {
    ab <- exp(theta) # Parameters are the *logs* of alpha and beta
    fit <- f.beta(ab[1], ab[2], x, ...)
    return (delta(fit, prob))
}
################################

#INPUT
x <- c(.3, .8)
x.p <- (function(p) log(p/(1-p)))(c(0.025, 0.975))

#SOLVE
start <- log(c(1e1, 1e1))
sol <- nlm(objective, start, x=x, prob=x.p, lower=0, upper=1,
                        typsize=c(1,1), fscale=1e-12, gradtol=1e-12)
parms <- exp(sol$estimate)

#ANSWER
parms
#[1] 7.597429 6.016240

#TEST
qbeta(p = c(.025, .975), parms[1], parms[2])
#[1] 0.3 0.8