How can I (numerically) approximate values for a beta distribution with large alpha & beta 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?
 A: A quick graphical experiment suggests that the beta distribution looks very like a normal distribution when alpha and beta are both very large. By googling "beta distribution limit normal" i found http://nrich.maths.org/discus/messages/117730/143065.html?1200700623$^\dagger$, which gives a handwaving 'proof'.
The wikipedia page for the beta distribution gives its mean, mode (v close to mean for large alpha and beta) and variance, so you could use a normal distribution with the same mean & variance to get an approximation. Whether it's a good enough approximation for your purposes depends on what your purposes are.

$^\dagger$ The link is broken.
A: I am going to infer you want an interval $[l,r]$ such that the probability that a random draw from the Beta RV is in the interval with probability 0.99, with bonus points for $l$ and $r$ being symmetric around the mode. By Gauss' Inequality or the Vysochanskii-Petunin inequality, you can construct intervals that contain the interval $[l,r]$, and would be fairly decent approximations. For sufficiently large $\alpha, \beta$, you will have numerical underflow issues in even representing $l$ and $r$ as distinct numbers, so this route may be good enough.
A: A normal approximation only gets better, if you work on the logit scale (i.e. after transforming the probabilities using the logit-function $\text{logit}(\text{probability}) := \log \text{probability} - \log(1-\text{probability})$. You then use a normal distribution with mean $\psi(\alpha)-\psi(\beta)$, where $\psi()$ denotes the digamma-function and with variance $\psi^{(1)}(\alpha)+\psi^{(1)}(\beta)$, where $\psi^{(1)}()$ denotes the trigamma function (see Wikipedia).
If you want xx% distribution percentiles, you take those of this normal distribution and then backtransform them to the proportion scale using the $\text{expit}$ function (i.e. $\frac{e^x}{1+e^x}$, aka the inverse $\text{logit}$ function).
A: A Normal approximation works extremely well, especially in the tails.  Use a mean of $\alpha/(\alpha+\beta)$ and a variance of $\frac{\alpha\beta}{(\alpha+\beta)^{2} (1+\alpha+\beta)}$.  For example, the absolute relative error in the tail probability in a tough situation (where skewness might be of concern) such as $\alpha = 10^6, \beta = 10^8$ peaks around $0.00026$ and is less than $0.00006$ when you're more than 1 SD from the mean.  (This is not because beta is so large: with $\alpha = \beta = 10^6$, the absolute relative errors are bounded by $0.0000001$.)  Thus, this approximation is excellent for essentially any purpose involving 99% intervals.
In light of the edits to the question, note that one does not compute beta integrals by actually integrating the integrand: of course you'll get underflows (although they don't really matter, because they don't contribute appreciably to the integral). There are many, many ways to compute the integral or approximate it, as documented in Johnson & Kotz (Distributions in Statistics). An online calculator is found here.  You actually need the inverse of this integral.  Some methods to compute the inverse are documented on the Mathematica site at http://functions.wolfram.com/GammaBetaErf/InverseBetaRegularized/.  Code is provided in Numerical Recipes. A really nice online calculator is the Wolfram Alpha site (www.wolframalpha.com ): enter inverse beta regularized (.005, 1000000, 1000001) for the left endpoint and inverse beta regularized (.995, 1000000, 1000001) for the right endpoint ($\alpha=1000000, \beta=1000001$, 99% interval).
A: If $p$ is a beta distributed variable, then it is the log-odds of $p$ (ie: $log(p/(1-p))$ that is approximately normally distributed. This is true even for highly skewed beta distributions as along as $min(\alpha,\beta) > 100$
For example
f <- function(n, a, b) {
    p <- rbeta(n, a, b)
    lor <- log(p/(1-p))
    ks.test(lor, 'pnorm', mean(lor), sd(lor))$p.value
}
summary(replicate(50, f(10000, 100, 1000000)))

typically produces an output like

summary(replicate(50, f(10000, 100, 1000000)))
     Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  0.01205 0.10870 0.18680 0.24810 0.36170 0.68730 

ie typical p-values are around 0.2.
So even with 10000 samples the Kolmogorov-Smirnov test lacks the power to distinguish the log odds ratio transformation of a highly skewed beta distributed variable with $\alpha=100, \beta=100000$.
However a similar test on the distribution of $p$ itself
f2 <- function(n, a, b) {
    p <- rbeta(n, a, b)
    ks.test(p, 'pnorm', mean(p), sd(p))$p.value
}
summary(replicate(50, f2(10000, 100, 1000000)))

produces something like
summary(replicate(50, f2(10000, 100, 1000000)))
     Min.   1st Qu.    Median      Mean   3rd Qu.      Max. 
2.462e-05 3.156e-03 7.614e-03 1.780e-02 1.699e-02 2.280e-01 

with typical p-values around 0.01
The R qqnorm function also gives a helpful visualization, producing  a very straight looking plot for the log-odds distribution indicating approximate normality the distribution of the beta dsitribute variable produces a distinctive curve indicating non normality
Therefore it is reasonable to use a Gaussian approximation in log-odds space
even for highly skewed $\alpha,\beta$ values as long as both are over 100.
