Sat 16 Jun 2012 07:50:47 AM UTC, comment #10:
Yes, it is similar to erfinv. I adapted the erfinv function to the domain of erfcinv.
|
Fri 15 Jun 2012 03:27:45 PM UTC, comment #9:
Axel, do you have a reference for the algorithm that I should put at the head of the code? It looks very similar to the existing algorithm from P.J. Acklam which is used for erfinv. If it is, then I can change the note there to say that the citation covers both erfinv and erfcinv.
|
Fri 15 Jun 2012 07:01:29 AM UTC, comment #8:
You're welcome. Sorry for the cleanup you had to do, here is another useless line I just spotted.
(file #26036)
|
Fri 15 Jun 2012 05:37:49 AM UTC, comment #7:
I committed your patch here (http://hg.savannah.gnu.org/hgweb/octave/rev/cb85e836d035). I also added your name to the list of contributors. Thanks for this work.
|
Thu 14 Jun 2012 09:59:13 PM UTC, comment #6:
Here is a better version of the patch: complies the guidelines, NEWS and unimplemented.m updated, additional tests on mappers.cc.
The algorithm is conceptually similar to erfinv, but no subtraction is performed on small input values.
(file #26033)
|
Thu 14 Jun 2012 03:54:51 PM UTC, comment #5:
This is a good patch. In addition, it adds the erfcinv function which Octave has been missing. Could you also remove 'erfcinv' from the list of missing functions in unimplemented.m and add it to the list of new functions present in 3.8.0 in the NEWS file?
Finally, do you have a reference source for the code which implements erfcinv? Is it sufficiently accurate for double precision IEEE-854 floating point?
|
Thu 14 Jun 2012 03:32:09 PM UTC, comment #4:
Thanks for the patch.
Can you produce it according to these guidelines?
http://www.octave.org/wiki/index.php?title=Commit_message_guidelines
http://jordi.inversethought.com/blog/how-to-write-a-patch-for-octave/
If you can't do this, I or someone else may do it, but someone has to do it. It's a bit more trouble for us to do it, and since this bug doesn't affect me personally, I have a little less motivation do this, and I might not apply it in a timely manner.
I hope you don't find this extra work too onerous...
|
Thu 14 Jun 2012 03:03:57 PM UTC, comment #3:
Here is a patch where I implemented erfcinv without calling erfinv, and I applied it to norminv as stated in comment #2. It solves this bug.
(file #26030)
|
Thu 07 Jun 2012 07:30:10 PM UTC, comment #2:
In Matlab, I get:
norminv can also be calculated using erfcinv, as
which avoids the subtraction you mention.
Unfortunately, base Octave doesn't implement erfcinv, while the one implemented in specfun simply defines it in terms of erfinv, hence it recreates the aforementioned problem.
|
Thu 07 Jun 2012 03:31:31 PM UTC, comment #1:
Confirmed. The problem can be traced to the erfinv() function.
Octave attempts to compute the inverse of the Gaussian using
On machines using IEEE-854 8-byte double precision floating point numbers the machine precision (eps()) is 2.2e-16. When you have a calculation, such as subtraction, where one nomber is of order 1 and the other is below eps, such as 1e-17, then the smaller number is effectively rounded to zero. This leads to
Could you run the following code in Matlab to help elucidate things?
|
Thu 07 Jun 2012 01:48:52 PM UTC, original submission:
For small inputs, norminv in Octave gives a different result from Matlab.
For very small inputs, norminv in Octave gives -Inf, while Matlab gives a valid answer.
For example, in Octave I get
while in Matlab I get
I am using Octave 3.6.1 compiled from source on Ubuntu 11.10, and Matlab R2010a.
|