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

[josef.pktd at gmail.com]

On Wed, Sep 16, 2009 at 8:14 PM, David Warde-Farley <dwf at cs.toronto.edu> wrote:
>
> On 16-Sep-09, at 6:40 PM, josef.pktd at 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.
>>
> Apparently the 'mpmath' package (easy_install'able) does arbitrary
> precision, and has a built-in normal cdf:
>
> Out[34]:
> mpf
> ('7.6198530241605260659733432515993083635040332779569605780353524e-24')
>
> In [35]: mpmath.functions.ncdf(-10.2)
> Out[35]:
> mpf
> ('9.9136251225600715794942192071914050771727823887948971717481189e-25')
>
> In [36]: mpmath.functions.ncdf(-20)
> Out[36]:
> mpf
> ('2.7536241186062336950756227808574653328074977347593305676993728e-89')
>
> In [37]: mpmath.functions.npdf(-20.2) / mpmath.functions.ncdf(-20)
> Out[37]:
> mpf('0.35995251388497550054560007939076980260015643935329865252047705')
>
>
> Looks pretty good to me.
>
> David

Thanks, I have mpmath installed for a while but haven't done much with
it. I thought initially it might be a good benchmark for testing
scipy.special.

The normal cdf at the lower tail in scipy.special looks really useful.
It would be possible to build a truncated normal tail distribution
with it. The conditional density looks precise enough. And it's
possible to use symmetry to build the upper tail.

more fun with distributions ?

Josef

>>> mpmath.functions.npdf(-20.2) / mpmath.functions.ncdf(-20)
mpf('0.3599525138849755005456000793907698026001564393532986525204771176091088952825768140613952325481008790977')
>>> from scipy import stats
>>> stats.norm.pdf(-20.2) / stats.norm.cdf(-20)
0.35995251388498711
>>> stats.norm.pdf(-30.2) / stats.norm.cdf(-20)
1.3002766730385739e-110
>>> mpmath.functions.npdf(-30.2) / mpmath.functions.ncdf(-20)
mpf('1.300276673038541314317767156123948362259642213373797768631623741891659293976807121550129087917318179814e-110')
>>> stats.norm.pdf(-21.2) / stats.norm.cdf(-20)
3.684252913218099e-010
>>> mpmath.functions.npdf(-21.2) / mpmath.functions.ncdf(-20)
mpf('0.0000000003684252913218045966433287833404555465054012445296205154822975165088200949391675713389170126995681295188')
>>> stats.norm.pdf(-20.02) / stats.norm.cdf(-20)
13.437063718396516
>>> mpmath.functions.npdf(-20.02) / mpmath.functions.ncdf(-20)
mpf('13.43706371839603326063511929219320799476797708512756936002002696323573503280240684590166610062394498302')