[SciPy-dev] scipy.stats: sf for symmetric distributions.
David Warde-Farley
dwf@cs.toronto....
Wed Sep 16 19:02:38 CDT 2009
On 16-Sep-09, at 6:40 PM, josef.pktd@gmail.com wrote:
> David,
> Since you have multi-precision available with sage, could you check
> the precision of the normal cdf in the lower tail?
>
> In http://projects.scipy.org/scipy/ticket/962 there was a discussion
> about using the truncated normal in the tails. In the last comment
> there, I mentioned that it would be possible to use the lower tail if
> the precision of scipy.special in this case is high enough. But I
> wasn't able to verify how good the normal cdf really is.
>
> These numbers are in the range 1e-20 to 1e-80, so here precision
> really matters.
Sadly, Sage works out the integral involving the erf() function, which
is only implemented for double-precision.
I did however do it with numerical_integral which provides an error
bound (not a very good one, it turns out):
sage: numerical_integral(1/sqrt(2*pi) * exp(-x^2/2), -infinity, -10.2)
(9.9136251237219932e-25, 6.0437083422386946e-28)
In [3]: norm.cdf(-10.2)
Out[3]: 9.9136251225599674e-25
sage: numerical_integral(1/sqrt(2*pi) * exp(-x^2/2), -infinity, -20)
(2.7536241192643532e-89, 1.2525306363698311e-92)
In [4]: norm.cdf(-20)
Out[4]: 2.7536241186061549e-89
sage: numerical_integral(1/sqrt(2*pi) * exp(-x^2/2), -infinity, -10)
(7.6198530279086935e-24, 5.1663456001768262e-27)
In [5]: norm.cdf(-10)
Out[5]: 7.619853024160474e-24
We seem to be getting about ten digits of agreement, but that could be
an artifact.
sage: numerical_integral(1/sqrt(2*pi) * exp(-x^2/2), -infinity, -30)
(4.906713926779639e-198, 3.680071161187036e-199)
In [6]: norm.cdf(-30)
Out[6]: 4.906713927147908e-198
(gotta love that bound)
Maybe GSL has an arbitrary precision erf taht we could compare against?
David
More information about the Scipy-dev
mailing list