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

 
 

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

Submitted by:  None
Submitted on:  Thu 07 Jun 2012 01:48:52 PM UTC  
 
Category: LibrariesSeverity: 3 - Normal
Priority: 5 - NormalItem Group: Inaccurate Result
Status: FixedAssigned to: None
Originator Name: Originator Email: -unavailable-
Open/Closed: ClosedRelease: 3.6.1
Operating System: GNU/Linux

Add a New Comment(Rich Markup)
   

You are not logged in

Please log in, so followups can be emailed to you.

 

(Jump to the original submission 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>
Project 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>
Project 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>
Project 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>
Project Administrator
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:

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.

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

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?

Rik <rik5>
Project 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

while in Matlab I get

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 File(s):
   
   
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
  • -unavailable- added by jordigh (Posted a comment)
  • -unavailable- added by amathei (Updated the item)
  • -unavailable- added by rik5 (Posted a comment)
  •  

    Do you think this task is very important?
    If so, you can click here to add your encouragement to it.
    This task has 0 encouragements so far.

    Only project members can vote.

     

    Please enter the title of George Orwell's famous dystopian book (it's a date):

     

     

    Follow 9 latest changes.

    Date Changed By Updated Field Previous Value => Replaced By
    Fri 15 Jun 2012 07:01:29 AM UTCamatheiAttached File-=>Added erfcinvcleanup.diff, #26036
    Fri 15 Jun 2012 05:37:49 AM UTCrik5StatusPatch Submitted=>Fixed
      Open/ClosedOpen=>Closed
    Thu 14 Jun 2012 09:59:13 PM UTCamatheiAttached File-=>Added erfcinv.diff, #26033
    Thu 14 Jun 2012 03:32:09 PM UTCjordighStatusConfirmed=>Patch Submitted
    Thu 14 Jun 2012 03:03:57 PM UTCamatheiAttached File-=>Added erfcinv.patch, #26030
    Thu 07 Jun 2012 03:31:31 PM UTCrik5CategoryNone=>Libraries
      Item GroupNone=>Inaccurate Result
      StatusNone=>Confirmed

    Back to the top


    Powered by Savane 3.1-cleanup1