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

Anne Archibald peridot.faceted@gmail....
Tue Dec 22 15:49:39 CST 2009


2009/12/22 Warren Weckesser <warren.weckesser@enthought.com>:
> Robert Kern wrote:
>> On Tue, Dec 22, 2009 at 14:01, Anne Archibald
>> <aarchiba@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.

That does look pretty dodgy, since the fault that happens is that one
of the bounds is returned, but it's the wrong one (top bound rather
than bottom bound). I'll experiment with it some more once I can get
scipy to compile again (the new lambertw won't compile; I think it's a
cython version issue).

Anne


More information about the SciPy-Dev mailing list