Re: Overhaul of statistical distribution functions

[Rik]

On 09/20/2011 11:38 PM, Michael D Godfrey wrote:

Rik,

I only looked quickly at binopdf, but it appeared to me that
you may not have seen the attached paper and algorithm.
This algorithm is pretty widely used, including in R.

9/21/11

Michael,

I cleaned up a lot, but wasn't always trying to re-implement algorithms which looked reasonably solid.  The Loader paper looks good though, and particularly for large N as you said.  Are there any copyright issues with filing a bug report on the Octave bug tracker and attaching the Loader paper?  I just don't have time to implement it myself.

In this regard, there are still several .m functions which use iterative algorithms which are going to be slow in Octave.  These are:

binoinv
nbininv
poissinv

betainv
gaminv

hygecdf
hygeinv
hygernd

The first three are problematic because they use a 'while (1)' construction and the stopping condition is not always met. For example,

poissinv(1-eps,4)

creates an infinite loop that requires a Ctrl-C to kill.  At the very least, it would be good to change these to for loops so that they would stop eventually. 

I would rather rewrite those entirely to avoid the iterative nature, but to do so I need a function to compute the inverse of the incomplete Beta function.  Matlab has such a function called betaincinv (http://www.mathworks.com/help/techdoc/ref/betaincinv.html).  Do you have reference papers or code for such a function?  A little bit of searching found this http://people.sc.fsu.edu/~jburkardt/f_src/asa109/asa109.html which has implementations in C++, Fortran, or Matlab.  Given the .edu address it may be possible to contact the author and get agreement to use the code.  I'm out of time to pursue that angle though which is why I was going to file an issue report about those functions.

--Rik