Mon 08 Feb 2016 04:08:51 AM UTC, comment #14:
A brief look at the hygeXXX case convinced me that they won't be easy to vectorise efficiently. It requires simultaneously working on PDFs of different lengths. If each is made a column of a 2D array, then there can be excessive padding of the short columns. If they are concatenated to a single long array, there can be catastrophic loss of precision when taking cumsum.
One approach would be to code these entirely in C++. A less drastic change would be to implement a "cumsum with reset points" in C++ and use that.
All supplied patches have been submitted, so I'm changing the status back to "none", and making the summary specific to hypergeometric distributions. I'm also tentatively lowering the severity and priority.

Mon 12 Oct 2015 04:25:58 AM UTC, comment #13:
I reviewed and checked in binoinv.m (http://hg.savannah.gnu.org/hgweb/octave/rev/9d2023d1a63c) and nbininv.m (http://hg.savannah.gnu.org/hgweb/octave/rev/d6d04088ac9e). I had to update some of the code to follow Octave coding conventions. In particular, using a space after a function and before the opening parenthesis as in "sin (x)" rather than "sin(x)". The new BIST tests also didn't run at first and required some debugging.
The speedup for binoinv was 28X and an incredible 85X for nbininv.
For gaminv and betainv, I rewrote these to use dountil loops with a maximum loop count of 40, rather than 10,000. If the results ever do fail to converge within 40 cycles a warning will be issued. But at least these won't "hang" anymore while looping from 40 to 10,000.
I'll leave this bug report open for the hygeXXX.m functions which remain.
As an example,
which shows that the scalar case is much, much faster than the vector case.

Tue 15 Sep 2015 11:28:47 AM UTC, comment #12:
Here are drafts for new binoinv.m and nbininv.m.
In version 4.0.0, binoinv.m didn't have a while(1) loop and so was "safe". It would run at most n iterations. However, the current version is about 10 times faster, on my tests.
These are both probably overkill given how often they are used, but it was fun writing them :)
I chose not to use Gile's algorithm for these. The quick read I had of the paper Rik found suggested that it is only useful for n>500. That paper was by a summer student  presumably an undergrad getting a feel for research  and so I think it was a case of "Look! We can apply this technique more broadly", rather than aiming for something practical.
Regarding the other distributions, I think betainv.m and gaminv.m are doing a Newton search, which is pretty efficient. It can't be vectorised because each iteration relies on the result of the previous. I recommend we keep them as they are. (They may have been updated since the original bug report.)
I don't understand the hygeXXX.m code well enough to comment.
(file #34889, file #34890)

Tue 18 Aug 2015 07:32:41 PM UTC, comment #11:
The two remaining functions which are critical to upgrade are
binoinv
nbininv
Both of these use while (1) loops which may not converge.
Interestingly, there seems to have been recent progress on this exact issue. There is an academic project from 2013 which calculates new ways of calculating the inverse binomial function (http://www.lancaster.ac.uk/postgrad/moorhead/Work/SummerProject13.pdf). What is doubly interesting, is that the research is based on Mike Gile's approximation to the inverse Poission distribution which Lachlan just implemented. I didn't have time to read and digest the whole paper, but Appendix B contains Matlab code for the calculation already. Presumably since this is an academic publication we might ask the author for permission to use and then just change the coding style to match Octave conventions.
Appendix C contains Matlab code for the general calculation of the inverse of a CDF. If this is accurate, we might be able to extend this to all of the other missing inverse functions.

Tue 18 Aug 2015 06:40:45 PM UTC, comment #10:
@Lachlan: Thanks, the new version works perfectly. Where the old and new implementations disagree, I prefer the new implementation as it is closer to the expected value.
The test script I uploaded also had benchmarking code in it. I find that for scalar lambda the speedup is dramatic, ~6X the old implementation or 80% of the run time. For lambda as an array the speedup is smaller because there is still a for loop involved. Nevertheless, the improvement is still 33%.
I checked in the changes here http://hg.savannah.gnu.org/hgweb/octave/rev/b3f56ed8d1f9.

Tue 18 Aug 2015 09:21:13 AM UTC, comment #9:
Lachlan: I looked through your code and just want
to congratulate you on a job welldone.
And, of course, is there any chance that you would
be willing, and have the time, to work on another
of the "needy" functions that Rik listed?
By the way, the documentation of this function in the
Manual is very brief. Feel free to add some detail if
you think that it would help.
In any case, thanks for your help.
Michael

Tue 18 Aug 2015 06:46:59 AM UTC, comment #8:
Thanks for the extra info, Rik.
The bug is fixed. The problem was that I was assuming X to be a vector, not matrix. I used if (any (k)) instead of if (any (k(:))).
Previously, I only tested it with vectors. I tried your script with a 100x100x100 case, and it gave eight outbyone errors. I've also tried it when the mean is a matrix and X is a scalar, and both are matrices, and both are three dimensional.
(file #34672)

Mon 17 Aug 2015 03:01:33 AM UTC, comment #7:
@Lachlan: Here is a sample debugging session showing an offby2 error. What's very interesting is that if poissinv is called standalone on the failing example it works. This seems to indicate that there is some state involved.

Sun 16 Aug 2015 07:45:36 AM UTC, comment #6:
Thanks, Rik.
I had already spotted discrepancies by 1. I should have commented on them. The are an offbyepsilon error, as follows.
The Poisson distribution is discrete. The answer is always an integer, and so if the answer is wrong at all, it will be wrong by an integer.
These functions are finding n = argmin (F(n) < x). The error occurs when F(n)=x+epsilon or F(n)=xepsilon, and the n we return is off by one.
This "epsilon" isn't the machine eps. It is the error associated with the asympmtotic approximation used by the algorithm, which will be a bit larger.
I'll check out the cases with a discrepancy of 2 when I get a chance.

Sun 16 Aug 2015 03:32:54 AM UTC, comment #5:
I reviewed and made some minor changes to follow Octave coding guidelines. Appendix C and D of the documentation have some of the guidelines we use. Things I changed
1) Use '##' for comments that standalone on a line
2) Use '#' for comments that are appended to a line of code
3) Use parentheses around the test clause in an if statement
4) Use space between function name and opening parenthesis: 'find (x)' rather than 'find(x)'
I've uploaded the latest version as poissinv3.m. I've also uploaded a file I was using to test the new implementation against the old one. I found that they occasionally disagree by an integer (1 or 2) and not just a small number on the order of eps. Could you investigate? I'm pretty sure it's something small having to do with the cutover between the limit and the analytic approximation, but I can't check the new implementation in until I'm certain it is producing correct results.
(file #34654, file #34655)

Sat 15 Aug 2015 04:10:11 AM UTC, comment #4:
Thanks, Rik. There's no rush to review the patch, and thanks for checking Google Scholar; I'm trying to help, not make work :)

Tue 11 Aug 2015 06:58:50 PM UTC, comment #3:
Thank you Lachlan for the patch. I hope I can review it within the next two weeks (awfully busy and all of Octave is volunteerrun).
Since Octave is used widely we generally try to write basic good code that will perform well under a range of circumstances. This includes doing things like eliminating for loops which is good practice regardless of the particular value of lambda. A direct answer to your question is that I don't know what value of lamda is more common. I tried doing a search on scholar.google.com to sample a number of academic papers and just see what values of lamda are common in the literature, but it didn't work very well.

Mon 10 Aug 2015 02:16:57 AM UTC, comment #2:
Update: The file poissinv_newnew.m is an improvement on poissinv_new.m with the following changes:
 For scalar lambda, it avoids the loop for small cases
 it is optimised for small lambda instead of large lambda, by doing the explicit summation before doing the asymptotic approximation.
It is about twice as fast as poissinv_new.m when the answer is less than 20, but slower when the answer is greater than 20 (but is still several orders of magnitude faster than poissinv.m in that case).
Do we have any usage data to say whether small or large lambda is more common?
Do we have any usage data to say which of these functions is used most often?
(file #34618)

Sat 08 Aug 2015 02:13:50 PM UTC, comment #1:
The file poissinv_new.m is a replacement for poissinv.m that is much faster for large lambda, and has no infinite loop.
(file #34596)

Fri 23 Sep 2011 04:50:54 AM UTC, original submission:
The following functions in scripts/statistics/distributions are implemented with loops which are slow in Octave.
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(1eps,4)
will create an infinite loop that can only be broken with CtrlC.
At the very least, the 'while (1)' construction should be replaced with a for loop which will eventually terminate.
The best possible outcome would be to recode these entirely to get rid of the loops.
