I'm using pbeta(), the cumulative probability distribution fuction in R, to find the probabilities of input values (p) from 0-1. I'm interested in just a subset of those input values that give probabilities between 0.025 and 0.975 (this is a confidence interval for my data).
My initial approach was to list 99 values between 0 and 1 and run pbeta() on them, and then extract those values which had outputs in the right range:
alpha <- 50 # Set values for the other two pbeta() parameters.
beta <- 100
p <- c(seq(from=0.01, to=0.99, by=0.01)) # Inputs must be between 0.01 and 0.99 inclusive, so evaluate 99 points in that range.
pbeta_outputs <- c(rep(0,99)) # A vector to hold the pbeta() output.
for (i in seq(1:length(p)) ) { # Fill pbeta_outputs using pbeta().
pbeta_outputs[i] <- pbeta(p[i], shape1=alpha, shape2=beta)
}
probabilities <- data.frame(p, pbeta_outputs) # Arrange inputs and outputs in a data frame.
probabilities_keep <- subset(probabilities, pbeta_outputs>=0.025 & pbeta_outputs<=0.975) # Subset to keep only inputs with outputs in the right range.
But, this is imprecise. Instead, I'd like to directly calculate the range of input values that would produce pbeta() output in the right range.
Something like:
successful_inputs <- find_input_range(pbeta(x, shape1=alpha, shape2=beta), 0.025:0.975)
Does a function like this exist?
It seems that what you are looking for is qbeta?
alpha <- 50
beta <- 100
# get lower bound
(lb <- qbeta(0.025, shape1 = alpha, shape2 = beta))
# [1] 0.2603782
# get upper bound
(ub <- qbeta(0.975, shape1 = alpha, shape2 = beta))
# [1] 0.4104985
pbeta(ub, shape1 = alpha, shape2 = beta) - pbeta(lb, shape1 = alpha, shape2 = beta)
# [1] 0.95