Fast Implementation of Inverse Incomplete Beta Function in C++

Question

I'm looking for an implementation of the Inverse Incomplete Beta Function, possibly already written in C++ or easy to implement myself. However, I need it to be FAST! As in, I'm going to be running this in the inner loop of an optimizer, so it would hopefully take under a couple hundred clock cycles.

There are already a couple of threads here, but in this case I'm willing to throw away a lot of accuracy for speed. Also, the domain is somewhat restricted, as I'm only using integer values for a and b.

More background on the problem: I'm giving an integer number of trials n and an integer k <= n of these trials that were successful. I'm assuming that the background distribution for the underlying probability of a successful trial is uniform in [0,1], so given that i've seen some number of trials and successes my posterior distribution should be a beta distribution. In a Bayesian model I'm essentially trying to find the pth percentile of likely underlying probabilities.

Again, I don't need this to be extremely accurate, just fast. I can deal with up to +/- 1% inaccuracy. However, it can't be extremely inaccurate for small numbers: my inputs range from nearly zero to tens of thousands.

Thanks in advance! If any clarification is needed let me know.

Answer 1

One approach is to make a table. If you need to differentiate you'll need to interpolate it. This is probably only an alternative for you if you keep all but one of the parameters fixed, but I think you do? (OP says no)

You may want to make the table binning non-linear to get a good accuracy at low x as you asked for. Try bin size proportional to x, x^2 etc.

2. Use a simple search method such as http://en.wikipedia.org/wiki/Secant_method to find your value as function of a power expansion of the (non inverse) incomplete beta function. http://en.wikipedia.org/wiki/Beta_function#Incomplete_beta_function . This works only if it is monotonic.

3. First make sure you actually need to make your own method. Perhaps try these first? http://www.boost.org/doc/libs/1_35_0/libs/math/doc/sf_and_dist/html/math_toolkit/special/sf_beta/ibeta_inv_function.html