I'm trying to implement a function while programming that finds an exact binomial test confidence given a probability of success, p, number of trials, k, and value of alpha. So for example given p=0.05, k=100 and alpha = 0.05 I want to retrieve something like (5,20) since $P(5 \geq X \geq 20)$ is ~0.04 which is about as close as you can get to alpha=0.05. Now there are several intervals which can be found that fit the test but you want one that makes the probability as close as possible to alpha (but not greater than) while also making sure that the probability on either side of the interval is roughly the same (so about $\frac{alpha}{2}$ on either side). Is there a good method to find the "best" (or at least a good) fit which can be easily implemented with some code (for example, while condition try these intervals).
EDIT: I came up with this solution, is it sound?
std::pair<unsigned long, unsigned long> BinomialConfidenceInterval(double p, unsigned long trials, double alpha)
{
unsigned long LowerBound = 0;
unsigned long UpperBound = trials;
unsigned long currentLowerBound = LowerBound;
unsigned long currentUpperBound = UpperBound;
while (CumulativeBinomProbability(currentLowerBound, currentUpperBound, p, trials) > (1-alpha))
{
if (CumulativeBinomProbability(LowerBound, currentLowerBound + 1, p, trials) > CumulativeBinomProbability(currentUpperBound - 1, UpperBound, p, trials))
currentUpperBound--;
else
currentLowerBound++;
}
std::pair<unsigned long, unsigned long> interval{ currentLowerBound , currentUpperBound };
return interval;
}