bugGNU Octave - Bugs: bug #36607, norminv doesn't work for very...

 
 

bug #36607: norminv doesn't work for very small values

Submitter:  None
Submitted:  Thu 07 Jun 2012 01:48:52 PM UTC
   
 
Category:  Libraries Severity:  3 - Normal
Priority:  5 - Normal Item Group:  Inaccurate Result
Status:  Fixed Assigned to:  None
Originator Name:  Originator Email:  -email is unavailable-
Open/Closed:  * Closed Release:  * 3.6.1
Operating System:  * GNU/Linux Fixed Release:  None
Planned Release:  None
* Mandatory Fields

Add a New Comment Rich Markup
   

Jump to the original submission

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.

Axel Mathéi <amathei>
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.

Rik <rik5>
Group administrator
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)

Axel Mathéi <amathei>
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.

Rik <rik5>
Group administrator
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)

Axel Mathéi <amathei>
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?

Rik <rik5>
Group administrator
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...

Jordi Gutiérrez Hermoso <jordigh>
Group Member
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)

Axel Mathéi <amathei>
Thu 07 Jun 2012 07:30:10 PM UTC, comment #2: 

In Matlab, I get:


erfinv (-1 + [1e-16 1e-17 1e-18 1e-19 1e-20])
ans = -5.8639      -Inf      -Inf      -Inf      -Inf

norminv ([1e-16 1e-17 1e-18 1e-19 1e-20])
ans = -8.2221   -8.4938   -8.7573   -9.0133   -9.2623


norminv can also be calculated using erfcinv, as


norminv = - sqrt (2) * erfcinv (2 * x);


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.

Anonymous
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


inv = sqrt (2) * erfinv (2 * x - 1);


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


erfinv (-1) = -Inf


Could you run the following code in Matlab to help elucidate things?


erfinv (-1 + [1e-16 1e-17 1e-18 1e-19 1e-20])
norminv ([1e-16 1e-17 1e-18 1e-19 1e-20])



Rik <rik5>
Group administrator
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

norminv(10^(-16))
ans = -8.2095
norminv(10^(-17))
ans = -Inf


while in Matlab I get

norminv(10^(-16))
ans = -8.2221
norminv(10^(-17))
ans = -8.4938


I am using Octave 3.6.1 compiled from source on Ubuntu 11.10, and Matlab R2010a.


Anonymous

 

(Note: upload size limit is set to 16384 kB, after insertion of the required escape characters.)

Attach Files:
   
   
Comment:
   

Attached Files
file #26036:  erfcinvcleanup.diff added by amathei (676B - text/x-patch)
file #26033:  erfcinv.diff added by amathei (10KiB - text/x-patch - erfcinv implemented)
file #26030:  erfcinv.patch added by amathei (8KiB - text/x-patch - erfcinv implemented without calling erfinv; erfcinv used in norminv instead of erfinv)

 

Depends on the following items: None found

Items that depend on this one: None found

 

Carbon-Copy List
  • -email is unavailable- added by jordigh (Posted a comment)
  • -email is unavailable- added by rik5 (Posted a comment)
  •  

    There are 0 votes so far. Votes easily highlight which items people would like to see resolved in priority, independently of the priority of the item set by tracker managers.

    Only group members can vote.

     

    Follow 9 latest changes.

    Date Changed by Updated Field Previous Value => Replaced by
    2012-06-15 amathei Attached File- Added erfcinvcleanup.diff, #26036
    2012-06-15 rik5 StatusPatch Submitted Fixed
        Open/ClosedOpen Closed
    2012-06-14 amathei Attached File- Added erfcinv.diff, #26033
    2012-06-14 jordigh StatusConfirmed Patch Submitted
    2012-06-14 amathei Attached File- Added erfcinv.patch, #26030
    2012-06-07 rik5 CategoryNone Libraries
        Item GroupNone Inaccurate Result
        StatusNone Confirmed

    Back to the top

    Powered by Savane 3.13-02a9.
    Corresponding source code