Is there a numerically stable way to calculate values of a beta distribution for large integer alpha, beta (e.g. alpha,beta > 1000000)?
Actually, I only need a 99% confidence interval around the mode, if that somehow makes the problem easier.
Add: I'm sorry, my question wasn't as clearly stated as I thought it was. What I want to do is this: I have a machine that inspects products on a conveyor belt. Some fraction of these products is rejected by the machine. Now if the machine operator changes some inspection setting, I want to show him/her the estimated reject rate and some hint about how reliable the current estimate is.
So I thought I treat the actual reject rate as a random variable X, and calculate the probability distribution for that random variable based on the number of rejected objects N and accepted objects M. If I assume a uniform prior distribution for X, this is a beta distribution depending on N and M. I can either display this distribution to the user directly or find an interval [l,r] so that the actual reject rate is in this interval with p >= 0.99 (using shabbychef's terminology) and display this interval. For small M, N (i.e. immediately after the parameter change), I can calculate the distribution directly and approximate the interval [l,r]. But for large M,N, this naive approach leads to underflow errors, because x^N*(1-x)^M is to small to be represented as a double precision float.
I guess my best bet is to use my naive beta-distribution for small M,N and switch to a normal distribution with the same mean and variance as soon as M,N exceed some threshold. Does that make sense?