[Bug-gsl] [bug #47646] gsl_ran_beta returns NaN for small arguments

[Alexey Radul]

URL:
  <http://savannah.gnu.org/bugs/?47646>

                 Summary: gsl_ran_beta returns NaN for small arguments
                 Project: GNU Scientific Library
            Submitted by: axch
            Submitted on: Thu 07 Apr 2016 10:15:31 PM GMT
                Category: Accuracy problem
                Severity: 3 - Normal
        Operating System: Ubuntu 14.04
                  Status: None
             Assigned to: None
             Open/Closed: Open
                 Release: 1.16
         Discussion Lock: Any

    _______________________________________________________

Details:

Calling, for example, gsl_ran_beta(rng, 1e-5, 1e-5), returns NaN most of the
time.  Expected behavior is to return 0.0 or 1.0 most of the time (equal
probability of each).

Presumably this is due to both calls to gsl_ran_gamma(rng, 1e-5, 1.0) in the
body of gsl_ran_beta returning 0 fairly often, which is not unreasonable
behavior for sampling a gamma with so small of a shape parameter.

For an example alternative, numpy has a different control path when a and b
are both <= 1:
https://github.com/numpy/numpy/blob/master/numpy/random/mtrand/distributions.c#L188-L218

This came up in a computational investigation of a Polya-tree style
representation of a Dirichlet process (context: Bayesian non-parametrics).




    _______________________________________________________

Reply to this item at:

  <http://savannah.gnu.org/bugs/?47646>

_______________________________________________
  Message sent via/by Savannah
  http://savannah.gnu.org/