bugGNU Scientific Library - Bugs: bug #47646, gsl_ran_beta returns NaN for small...


You are not allowed to post comments on this tracker with your current authentication level.

bug #47646: gsl_ran_beta returns NaN for small arguments

Submitted by:  Alexey Radul <axch>
Submitted on:  Thu 07 Apr 2016 10:15:31 PM UTC  
Category:  Accuracy problem Severity:  3 - Normal
Operating System:  Ubuntu 14.04 Status:  Fixed
Assigned to:  None Open/Closed:  Open
Release:  1.16

Wed 01 Nov 2017 09:14:52 AM UTC, comment #5:

[From Yu]

I was thinking if it's worthwhile to allow the use of close-form cdf for some functions. Currently, the test routine calculates the expected probability in each bin using numerical integration of pdf. In this particular case, it blows up. Using cdf in stead of integration can avoid the problem, but I am not sure if this is appropriate. There are other functions with close-form cdf as well, and I am wondering whether it can be used to replace integration as well.
In a modified test routine, I added a flag to indicate whether the input function ``pdf" is actually a cdf, and if that is the case, the probability in each bin would be simply cdf[x+dx] - cdf[x]. Please see the relevant lines here https://github.com/ohliumliu/gsl-playground/blob/issue2/scratchpad/test_beta_small.c#L461-#L464 One obvious issue is that in this function, the input ``pdf" could mean pdf or cdf depending on the value of the flag. This may be confusing.
Without much experience in numerical computation, maybe I am proposing something crazy.


Patrick Alken <psa>
Project Administrator
Thu 15 Sep 2016 08:20:53 PM UTC, comment #4:

I'm also not sure how to properly test for this case, as the PDF integration test fails for these small arguments.

Patrick Alken <psa>
Project Administrator
Thu 15 Sep 2016 08:19:41 PM UTC, comment #3:

Patch applied in commit 05c5b5179d5f12e893198e54741cb4b8c289ad09

Patrick Alken <psa>
Project Administrator
Fri 15 Jul 2016 09:46:49 PM UTC, comment #2:

Following up based on numpy issues
(please let me know if this kind of information should be sent over the mailing list.)

numpy fixed this problem around March 2015
It was related to a Dirichlett distribution bug that is still open.

As for the code, the while loop is there to handle the very RARE cases of `u` and `v` being zero, because `rk_double` could give a zero. But with `gsl_ran_uniform_pos` this loop is redundant.

In terms of test, `numpy` code only check if there is `nan` generated. See numpy #5858 above.

Yu Liu <yuliu>
Fri 15 Jul 2016 06:10:22 AM UTC, comment #1:

The diff file is a straight implementation of the numpy implementation.

  • Testing for mean and variance (bernoulli_test.c) shows some small difference from expectation. Further tests are probably required.
  • I am still working on the tests, trying to use the test methods shipped with GSL. Because cdf seems to diverge at 0 and 1, testCDF errors out due to integration convergence. I am now trying to use beta_cdf directly in testPDF routine. Please let me know if there is an appropriate approach.
Yu Liu <yuliu>
Thu 07 Apr 2016 10:15:31 PM UTC, original submission:

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).

Alexey Radul <axch>


Attached Files
file #37915:  test_beta_small.c added by yuliu (42KiB - application/octet-stream - A more stringent test using cdf instead of integration to get expected counts.)
file #37884:  beta_distribution.diff added by yuliu (1KiB - application/octet-stream - Hello,)
file #37885:  bernoulli_test.c added by yuliu (2KiB - application/octet-stream - Hello,)


Depends on the following items: None found

Items that depend on this one: None found


Carbon-Copy List
  • -email is unavailable- added by psa (Posted a comment)
  • -email is unavailable- added by yuliu (Updated the item)
  • -email is unavailable- added by axch (Submitted the item)



    Follow 6 latest changes.

    Date Changed by Updated Field Previous Value => Replaced by
    2017-11-01 psa Open/ClosedClosed => Open
    2016-09-15 psa StatusNone => Fixed
        Open/ClosedOpen => Closed
    2016-07-16 yuliu Attached File- => Added test_beta_small.c, #37915
    2016-07-15 yuliu Attached File- => Added beta_distribution.diff, #37884
        Attached File- => Added bernoulli_test.c, #37885

    Back to the top

    Powered by Savane 3.3