bugGNU Octave - Bugs: bug #29543, normcdf lacks precision on the...

 
 

bug #29543: normcdf lacks precision on the tails of the distribution

Submitted by:  Guido Walter Pettinari <coccoinomane>
Submitted on:  Tue 13 Apr 2010 03:06:32 PM UTC  
 
Category: NoneSeverity: 3 - Normal
Priority: 5 - NormalItem Group: Incorrect Result
Status: FixedAssigned to: Jaroslav Hajek <highegg>
Originator Name: Open/Closed: Closed
Release: 3.2.3Operating System: Mac OS

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)

Thu 15 Apr 2010 07:12:18 PM UTC, comment #6:

The accuracy of our erf is fine, in fact it's normally handled by libc. The problem in fact is that the underlying formula
(1 + erf(x/sqrt(2)))/2 suffers from severe cancellation for larger negative x. All that is needed is to employ erfc instead.

http://hg.savannah.gnu.org/hgweb/octave/rev/e74bff13aa26

cheers

Jaroslav Hajek <highegg>In charge of this item.
Tue 13 Apr 2010 11:11:13 PM UTC, comment #5:

Well it depends.. normcdf depends on the accuracy of the underlying implementation of the erf function. Frankly we should fix that and the precision of normcdf would improve.. Again another need for GSL ;-) .... This really is a good topic to discuss on -unavailable- has there are many issues to resolve for the inclusion of GSL in octave and how to do that should be discussed before hand

Cheers
David

David Bateman <dbateman>
Project Member
Tue 13 Apr 2010 10:35:41 PM UTC, comment #4:

Dear David,

thank you for your answer.

I understand your point about GSL. Nevertheless, I suspect that 'normcdf' needs to be fixed anyway. Isn't that a bug the fact that 'normcdf' gives wrong results when the argument is farther than ~7-8 sigmas from the mean?

Cheers,

Guido

Guido Walter Pettinari <coccoinomane>
Tue 13 Apr 2010 10:27:08 PM UTC, comment #3:

This should be discussed on -unavailable- as part of a larger inclusion of GSL in Octave to replace many other functions and for this reason I'm flagging this as won't fix. There has been discussions a couple of times, but as

http://n4.nabble.com/precision-loss-with-besselh-bessely-for-order-n-16-td1672999.html#a1673001

points out GSL is double precision only and this would require extra manipulation to handle Octaves single precision type, which is perhaps undesirable.

In any case perhaps this should be added to the GSL bindings on octave-forge is this function isn't already included.

D.

David Bateman <dbateman>
Project Member
Tue 13 Apr 2010 06:10:08 PM UTC, comment #2:

UPDATE:
I wrote an oct-file function that wraps the GSL function to compute the normal CDF, i.e. "gsl_cdf_gaussian_P".

It works fine (it agrees with Mathematica and Matlab for all the cases I considered) and it is from 1.2 to 5 times faster than Octave's "normcdf".

You can find the relevant code attached; I compiled it with the following command:

mkoctfile gaussian_cdf.cc gaussian_pdfcdf.cc -L/opt/local/lib -lgsl -lgslcblas

Obviously you should specify the directory that contains GSL libraries instead of /opt/local/lib.

I hope it can be useful.

Cheers,

Guido

(file #20216)

Guido Walter Pettinari <coccoinomane>
Tue 13 Apr 2010 03:17:20 PM UTC, comment #1:

I am afraid I quoted a wrong value for normcdf(-9). The correct one according to both Mathematica & Matlab is:
normcdf(-9) = 1.12859e-19.

On the other hand, normcdf(-8) as calculated by both Mathematica & Matlab is:
normcdf(-8) = 6.22096e-16.

Octave's result is:
normcdf(-8) = 6.1062e-16

which is quite different from the correct one.

Cheers,

Guido

Guido Walter Pettinari <coccoinomane>
Tue 13 Apr 2010 03:06:32 PM UTC, original submission:

Hello,

results from 'normcdf' have a very low precision when given arguments that are more than 9 standard deviations far from the mean.

Example:

>> normcdf(-8, 0, 1)

ans = 6.1062e-16

>> normcdf(-9, 0, 1)

ans = 0

but normcdf(-9) = 6.221e-16, which is far away from float and double limits on my machine.

The same happens (i) with different combinations of mean and standard deviation and (ii) when calculating 1 - normcdf(9, 0, 1).

On the other hand, 'normpdf' seems to work properly.

Cheers,

Guido

Guido Walter Pettinari <coccoinomane>

 

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

Attach File(s):
   
   
Comment:
   

Attached Files
file #20216:  gaussian_cdf.zip added by coccoinomane (3KiB - application/zip - Substitution of normcdf - wrapper to GSL function "gsl_cdf_gaussian_P")

 

Depends on the following items: None found

Items that depend on this one: None found

 

Carbon-Copy List
  • -unavailable- added by highegg (Posted a comment)
  • -unavailable- added by dbateman (Posted a comment)
  • -unavailable- added by coccoinomane (Submitted the item)
  •  

    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 5 latest changes.

    Date Changed By Updated Field Previous Value => Replaced By
    Thu 15 Apr 2010 07:12:18 PM UTChigheggStatusWont Fix=>Fixed
      Assigned toNone=>highegg
      Open/ClosedOpen=>Closed
    Tue 13 Apr 2010 10:27:08 PM UTCdbatemanStatusNone=>Wont Fix
    Tue 13 Apr 2010 06:10:08 PM UTCcoccoinomaneAttached File-=>Added gaussian_cdf.zip, #20216

    Back to the top


    Powered by Savane 3.1-cleanup1