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

 
 

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

Submitter:  Guido Walter Pettinari <coccoinomane>
Submitted:  Tue 13 Apr 2010 03:06:32 PM UTC
   
 
Category:  None Severity:  3 - Normal
Priority:  5 - Normal Item Group:  Incorrect Result
Status:  Fixed Assigned to:  highegg
Originator Name:  Open/Closed:  * Closed
Release:  * 3.2.3 Operating System:  * Mac OS
Fixed Release:  None Planned Release:  None
* Mandatory Fields

Add a New Comment Rich Markup
   

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>
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 -email is 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>
Group 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 -email is 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>
Group 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 Files:
   
   
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
  • -email is unavailable- added by highegg (Posted a comment)
  • -email is unavailable- added by dbateman (Posted a comment)
  • -email is unavailable- added by coccoinomane (Submitted the item)
  •  

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

    Date Changed by Updated Field Previous Value => Replaced by
    2010-04-15 highegg StatusWont Fix Fixed
        Assigned toNone highegg
        Open/ClosedOpen Closed
    2010-04-13 dbateman StatusNone Wont Fix
    2010-04-13 coccoinomane Attached File- Added gaussian_cdf.zip, #20216

    Back to the top

    Powered by Savane 3.13-02a9.
    Corresponding source code