[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