[Scipy-tickets] [SciPy] #1614: logcdf function of normal distribution (scipy.stats) can not handle a wide enough range of values

SciPy Trac scipy-tickets@scipy....
Wed Mar 7 18:37:25 CST 2012

#1614: logcdf function of normal distribution (scipy.stats) can not handle a wide
enough range of values
 Reporter:  andrewschein             |       Owner:  somebody   
     Type:  defect                   |      Status:  new        
 Priority:  normal                   |   Milestone:  Unscheduled
Component:  scipy.stats              |     Version:  devel      
 Keywords:  normal distribution cdf  |  

Comment(by andrewschein):

 A cursory glance of the R source code indicates that the log.p is
 implemented by taking the log()
 of the CDF (as opposed to some direct computation).  The comments in the R
 code state:


  *  The "_both" , lower, upper, and log_p  variants were added by
  *  Martin Maechler, Jan.2000;
  *  as well as log1p() and similar improvements later on.


 It is possible that more than mere interface changes went into the work of
 implementing log_p, and this could explain the results that follow.

 Since it appears that the cdflib directory is not fully incorporated into
 the scipy build system, I copied a handful of the necessary .f files into
 a test directory and rolled my own module for the purpose of testing the
 cdfnor code.


 import cdfnor
 import numpy as np

 def wrapper(theta) :
     return cdfnor.cdfnor(1,0,theta,0,1,0,0)

 thetas = [0,-1,-2,-4,-8,-32,-64]
 for theta in thetas :
     print theta, wrapper(theta), np.log(wrapper(theta))


 produces the following output:


 ais@albus:temp[1]$ python test.py
 0 0.5 -0.69314718056
 -1 0.158655253931 -1.84102164501
 -2 0.0227501319482 -3.78318433368
 -4 3.16712418331e-05 -10.3601014865
 -8 6.22096057427e-16 -35.0134371599
 -32 5.45208060351e-225 -516.385648626
 Warning: divide by zero encountered in log
 -64 0.0 -inf

 So as written the cdfnor.f implementation does not match the precision of
 the R implementation.

Ticket URL: <http://projects.scipy.org/scipy/ticket/1614#comment:3>
SciPy <http://www.scipy.org>
SciPy is open-source software for mathematics, science, and engineering.

More information about the Scipy-tickets mailing list