[SciPy-dev] scipy.stats: sf for symmetric distributions.

[David Warde-Farley]

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