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

[Patrick Alken]

Follow-up Comment #5, bug #47646 (project gsl):

[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

___

Reply to this item at:

  

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

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

[Patrick Alken]

Update of bug #47646 (project gsl):

 Open/Closed:  Closed => Open   


___

Reply to this item at:

  

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

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

[Yu Liu]

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


On Thu, Sep 15, 2016 at 1:20 PM, Patrick Alken 
wrote:

> Follow-up Comment #4, bug #47646 (project gsl):
>
> I'm also not sure how to properly test for this case, as the PDF
> integration
> test fails for these small arguments.
>
> ___
>
> Reply to this item at:
>
>   
>
> ___
>   Message sent via/by Savannah
>   http://savannah.gnu.org/
>
>

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

[Patrick Alken]

Follow-up Comment #4, bug #47646 (project gsl):

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

___

Reply to this item at:

  

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

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

[Patrick Alken]

Update of bug #47646 (project gsl):

  Status:None => Fixed  
 Open/Closed:Open => Closed 

___

Follow-up Comment #3:

Patch applied in commit 05c5b5179d5f12e893198e54741cb4b8c289ad09

___

Reply to this item at:

  

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

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

[Yu Liu]

Follow-up Comment #2, bug #47646 (project gsl):

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.

___

Reply to this item at:

  

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

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

[Alexey Radul]

URL:
  

 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:

  

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