GNU 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
Submitter:  Alexey Radul <axch>  
Submitted:  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: 
Patrick Alken <psa> 
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> 
Thu 15 Sep 2016 08:19:41 PM UTC, comment #3: Patch applied in commit 05c5b5179d5f12e893198e54741cb4b8c289ad09 
Patrick Alken <psa> 
Fri 15 Jul 2016 09:46:49 PM UTC, comment #2: Following up based on numpy issues
numpy fixed this problem around March 2015
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.

Yu Liu <yuliu> 
Thu 07 Apr 2016 10:15:31 PM UTC, original submission:
Calling, for example, gsl_ran_beta(rng, 1e5, 1e5), 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, 1e5, 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#L188L218
This came up in a computational investigation of a Polyatree style representation of a Dirichlet process (context: Bayesian nonparametrics). 
Alexey Radul <axch> 
Depends on the following items: None found
Items that depend on this one: None found
Follow 6 latest changes.
Date  Changed by  Updated Field  Previous Value  =>  Replaced by 

20171101  psa  Open/Closed  Closed  Open  
20160915  psa  Status  None  Fixed  
Open/Closed  Open  Closed  
20160716  yuliu  Attached File    Added test_beta_small.c, #37915  
20160715  yuliu  Attached File    Added beta_distribution.diff, #37884  
Attached File    Added bernoulli_test.c, #37885 
[From Yu]
Patrick,
I was thinking if it's worthwhile to allow the use of closeform 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 closeform 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/gslplayground/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