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

[josef.pktd at gmail.com]

On Wed, Sep 16, 2009 at 6:24 PM, David Warde-Farley <dwf at cs.toronto.edu> wrote:
> On 16-Sep-09, at 5:25 PM, Pauli Virtanen wrote:
>
>>> sage: thepdf
>>> 1/4*sqrt(2)/(1/2*x^2 + 1)^(3/2)
>>> sage: n(2 * integral(thepdf, x, -infinity, -1), prec=100)
>>> 0.42264973081037423549085121950
>>
>> Btw, this integral evaluates to
>>
>>       1 - 1/sqrt(3)
>>
>> Now, depending on how you evaluate this in floating point, you get
>> different results:
>>
>>       1. - 1./sqrt(3.)        = 0.42264973081037416
>>       (sqrt(3.) - 1)/sqrt(3.) = 0.42264973081037421
>>       (3. - sqrt(3.))/3.      = 0.42264973081037427
>>
>> Perhaps the 0.###27 was obtained like this...
>
> It was obtained via 1 - scipy.special(2, 1),
>
>> It's nice if you can actually maintain this level of accuracy in
>> scipy.stats in a way that works on all platforms -- doesn't work in
>> scipy.special as the boost tests show.
>
> Well, I guess we're still in trouble (depending on the special
> function) but perhaps less trouble since we're not doing an
> unnecessary subtraction.
>
> Anyway, the ###21 result is closer than ###27 so I'm taking pleasure
> in (*extremely*) small things.
>
> David
> _______________________________________________
> Scipy-dev mailing list
> Scipy-dev at scipy.org
> http://mail.scipy.org/mailman/listinfo/scipy-dev
>

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.

Thanks,

Josef