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 problemSeverity: 3 - Normal
Operating System: Ubuntu 14.04Status: Fixed
Assigned to: NoneOpen/Closed: Open
Release: 1.16

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

[From Yu]

Patrick,
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.

Best,
yu

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
https://github.com/numpy/numpy/pull/5858
It was related to a Dirichlett distribution bug that is still open.
https://github.com/numpy/numpy/issues/5851

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
  • -unavailable- added by psa (Posted a comment)
  • -unavailable- added by yuliu (Updated the item)
  • -unavailable- added by axch (Submitted the item)
  •  

    Please enter the title of George Orwell's famous dystopian book (it's a date):

     

     

    Follow 6 latest changes.

    Date Changed By Updated Field Previous Value => Replaced By
    Wed 01 Nov 2017 09:15:14 AM UTCpsaOpen/ClosedClosed=>Open
    Thu 15 Sep 2016 08:19:41 PM UTCpsaStatusNone=>Fixed
      Open/ClosedOpen=>Closed
    Sat 16 Jul 2016 09:36:48 PM UTCyuliuAttached File-=>Added test_beta_small.c, #37915
    Fri 15 Jul 2016 05:58:53 AM UTCyuliuAttached File-=>Added beta_distribution.diff, #37884
      Attached File-=>Added bernoulli_test.c, #37885

    Back to the top


    Powered by Savane 3.1-cleanup1