[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
## [Help-gsl] gamma distribution

**From**: |
Sam Mason |

**Subject**: |
[Help-gsl] gamma distribution |

**Date**: |
Thu, 2 May 2013 18:02:02 +0100 |

Hi,
I've been using the GSL for a while and have just realised that I'm
working around bugs in the Gamma distribution's RNG and PDF. Some
parameters (i.e. small alpha and beta) cause the RNG to return zero.
As far as I understand this is incorrect mathematically and is infact
due to rounding errors. I'm OK with this "feature", but when you feed
this zero back to gsl_ran_gamma_pdf() (with the same parameters) you
get zero probabilityâ€“which feels wrong.
The implementation in R returns +Inf when alpha < 1, which seems more useful.
I'm currently working around this by adding DBL_MIN onto the value I
get back from gsl_ran_gamma(), but this doesn't seem very nice.
Would a patch that adjusts the behaviour of gsl_ran_gamma_pdf(), i.e.
changing the second "if" statement to:
else if (x == 0)
{
if (a == 1)
return 1/b ;
else if (a < 1)
return GSL_POSINF ;
else
return 0 ;
}
do the right thing in everybody's eyes?
Thanks,
Sam

[Prev in Thread] |
**Current Thread** |
[Next in Thread] |

**[Help-gsl] gamma distribution**,
*Sam Mason* **<=**