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

[Charles R Harris]

On Tue, Dec 22, 2009 at 2:49 PM, Anne Archibald
<peridot.faceted at gmail.com>wrote:

> 2009/12/22 Warren Weckesser <warren.weckesser at enthought.com>:
> > 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.
>
> 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).
>
>
There is a recent (about 1 mo old) version of cython out. Maybe we should
regenerate everything, although I'm not sure what the new release buys us, I
didn't see any release notes.

Chuck
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mail.python.org/pipermail/scipy-dev/attachments/20091222/ac96c600/attachment.html>