[SciPy-dev] scipy.special.bdtrik bug (ticket #1076)

[Warren Weckesser]

Robert Kern wrote:
> On Tue, Dec 22, 2009 at 14:01, Anne Archibald
> <aarchiba at physics.mcgill.ca> wrote:
>   
>> Hi,
>>
>> I was recently doing some calculations with scipy.stats.binom().ppf
>> and found a nasty bug ( http://projects.scipy.org/scipy/ticket/1076 ).
>> If the binomial probability is tiny, totally wrong answers emerge. The
>> problem turns out to be in the function scipy.special.bdtrik. There's
>> no documentation anywhere about what this is supposed to do, but from
>> context it's pretty clear it exists to calculate this value. It's a
>> cephes function, and I got a little lost trying to track down its
>> implementation. Maybe someone who's more familiar with cephes could
>> point me to the code, and how to put it in __all__?
>>     
>
> [~]$ cd svn/scipy/scipy/special
> [special]$ grin -i bdtrik
> ./_cephesmodule.c:
>   827 :         f = PyUFunc_FromFuncAndData(cephes3_functions,
> cdfbin2_data, cephes_4_types, 2, 3, 1, PyUFunc_None, "bdtrik", "", 0);
>   828 :         PyDict_SetItemString(dictionary, "bdtrik", f);
>   902 :         f = PyUFunc_FromFuncAndData(cephes3_functions,
> cdfnbn2_data, cephes_4_types, 2, 3, 1, PyUFunc_None, "nbdtrik", "",
> 0);
>   903 :         PyDict_SetItemString(dictionary, "nbdtrik", f);
> ....
>
> [special]$ grin cdfbin2
> ./_cephesmodule.c:
>   208 : static void * cdfbin2_data[] = {(void *)cdfbin2_wrap, (void
> *)cdfbin2_wrap};
>   827 :         f = PyUFunc_FromFuncAndData(cephes3_functions,
> cdfbin2_data, cephes_4_types, 2, 3, 1, PyUFunc_None, "bdtrik", "", 0);
> ./cdf_wrappers.c:
>    93 : double cdfbin2_wrap(double p, double xn, double pr) {
> ./cdf_wrappers.h:
>    18 : extern double cdfbin2_wrap(double p, double xn, double pr);
>
> [special]$ less cdf_wrappers.c
> # I see that cdfbin2_wrap() wraps the Fortran subroutine CDFBIN. This
> tells me that it's from the cdflib collection of functions, not the
> Cephes library itself.
>
>   

A quick look at the wrappers and the fortran function makes me think the 
bug is in the wrappers.  If the fortran function CDFBIN returns with 
STATUS == 1 or STATUS == 2, the wrapper returns BOUND, and the caller 
would only know something was wrong if  
scipy_special_print_error_messages is not zero.


Warren