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

josef.pktd@gmai... josef.pktd@gmai...
Wed Sep 16 19:59:54 CDT 2009


On Wed, Sep 16, 2009 at 8:46 PM,  <josef.pktd@gmail.com> wrote:
> On Wed, Sep 16, 2009 at 8:14 PM, David Warde-Farley <dwf@cs.toronto.edu> wrote:
>>
>> 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.
>>>
>> 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')
>

I checked the t distribution again. It is still missing the same trick for _isf.
David, could you check and, if correct, commit the following addition
to the t distribution.

2850	    def _ppf(self, q, df):
2851	        return special.stdtrit(df, q)

+	    def _isf(self, q, df):
+	        return -special.stdtrit(df, q)

Thanks,

Josef


More information about the Scipy-dev mailing list