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?

  • 1
    $\begingroup$ Do you want to know the mathematics or simply a code solution in R or some such? $\endgroup$ – John Aug 24 '10 at 10:40
  • $\begingroup$ I need to implement this in C#, so the mathematics would be good. A code sample would be fine too, if it doesn't rely on some builtin R/Matlab/Mathematica function I can't translate to C#. $\endgroup$ – nikie Aug 24 '10 at 10:55
  • $\begingroup$ PDF, CDF or inverse CDF? $\endgroup$ – J. M. is not a statistician Aug 24 '10 at 12:29
  • $\begingroup$ If you do not insist on Beta, you can use Kumaraswamy distribution that is very similar and has much simpler algebraic form: en.wikipedia.org/wiki/Kumaraswamy_distribution $\endgroup$ – Tim Sep 1 '16 at 21:16

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 at http://www.danielsoper.com/statcalc/calc37.aspx. 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 (www.nr.com). 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).

  • $\begingroup$ Perfect! I had the NR book on my desk all the time, but never thought to look there. Thanks a lot. $\endgroup$ – nikie Aug 24 '10 at 17:34

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, 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.

  • $\begingroup$ Stupid question: How did you do that graphical experiment? I tried to plot the distribution for alpha/beta around 100, but I couldn't see anything due to underflow errors. $\endgroup$ – nikie Aug 24 '10 at 15:57
  • $\begingroup$ You don't want to plot the integrand: you want to plot the integral. However, you can get the integrand in many ways. One is to enter "plot D (beta(x, 1000000, 2000000), x) / beta(1, 1000000, 2000000) from 0.3325 to 0.334" at the Wolfram Alpha site. The integral itself is seen with "Plot beta(x, 1000000, 2000000) / beta(1, 1000000, 2000000) from 0.3325 to 0.334". $\endgroup$ – whuber Aug 24 '10 at 17:27
  • $\begingroup$ I plotted the integrand, i.e. the pdf of the beta distribution, in Stata - it has a builtin function for the pdf. For large alpha and beta you need to restrict the range of the plot to see it's close to normal. If I was programming it myself i'd compute its logarithm then exponentiate at the end. That should help with the underflow problems. The beta function in the denominator is defined in terms of gamma functions, equivalent to factorials for integer alpha and beta, and many packages/libraries include lngamma() or lnfactorial() instead/as well as gamma() and factorial() functions. $\endgroup$ – onestop Aug 25 '10 at 12:51

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.

  • $\begingroup$ When alpha and beta are not too far apart (i.e., alpha/beta are bounded above and below), the SD of Beta[alpha, beta] is proportional to 1/Sqrt(alpha). E.g., for alpha = beta = 10^6, the SD is very close to 1/Sqrt(8) / 1000. I think there will be no problem with the representation of l and r even if you're only using single precision floats. $\endgroup$ – whuber Aug 24 '10 at 15:25
  • $\begingroup$ which is to say that $10^6$ is not 'sufficiently large' ;) $\endgroup$ – shabbychef Aug 24 '10 at 16:25
  • 1
    $\begingroup$ Yeah, it's a crazy number for a beta application. BTW, those inequalities won't produce good intervals at all, because they are extremes over all distributions (satisfying certain constraints). $\endgroup$ – whuber Aug 24 '10 at 17:19
  • $\begingroup$ @whuber: You're right, they are crazy numbers. With my naive algorithm, the "sane" numbers were easy and worked well, but I couldn't imagine how to calculate it for "crazy" parameters. Hence the question. $\endgroup$ – nikie Aug 24 '10 at 19:52
  • 2
    $\begingroup$ OK, you're right: once alpha+beta exceeds 10^30 or so, you will have difficulties with doubles :-). (But if you represent l and r as differences from the mean of alpha/(alpha+beta), you'll be fine until alpha or beta exceed about 10^303.) $\endgroup$ – whuber Aug 25 '10 at 15:34

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.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.