rabauke/trng4

Error in Gamma_dist with large kappa

Closed this issue · 10 comments

Deleted user commented

Hi,

First of all, awesome project, nice job!

I've found a bug that if we set the first parameter in Gamms_dist be large, the generator will generate 'nan'. The numbers I used are trng::gamma_dist<> g(1000.,1000.);

Thanks!

Deleted user commented

I read through gamma_dist.hpp, the possbile reason is that
math::pow(y, p.kappa()-1) goes to inf when kappa is large
Maybe changing the calculation to log scale will help?

rabauke commented

Thanks for your bug report. I will try to fix this for the next release.

Regards,
Heiko

Deleted user commented

Hi Heiko,

Thanks for replying.

Because of an urgency to use your library in my project, I tried modifying two lines by myself in gamma_dist.hpp this morning. Now the problem seems fixed. I paste them here just FYI:

//result_type f1=math::pow(y, p.kappa()-1)_math::exp(-y-ln_Gamma_kappa);
result_type f1=math::exp( (p.kappa()-1)_math::ln(y) - y - ln_Gamma_kappa);

//return math::pow(x, p.kappa()-1)/
// (math::exp(x+math::ln_Gamma(p.kappa()))_p.theta());
return math::exp( (p.kappa()-1)_math::ln(x)- x-math::ln_Gamma(p.kappa())) /
(p.theta());

Thanks!
Leo

rabauke commented

Thanks for your patch. I will incorporate it in the next release.

Regards,
Heiko

rabauke commented

This has been fixed in version 4.16.

oelhammouchi commented

Hi @rabauke, thank you very much for your work on this amazing package! I'm using it for an R package of my own, but I'm running into the same problem described here. In my case I'm using shape = 68628592.392190069 and scale = 0.00018490835317560905. This package is meant to deal with extreme scenario's, hence the outlying values.

Some quick troubleshooting in the debugger shows that the inverse CDF routine for the gamma distribution does not converge, returning a NaN. Would you be so kind as to look into this?

rabauke commented

Dear @OthmanElHammouchi, as mentioned above this has been addressed in version 4.16 some time ago. Which version are you using? Knowing the affected version will help to debug the issue.

oelhammouchi commented

@rabauke I'm using TRNG through the R interface package, I assume the versions are the same. As you can see in the output below, it's version 4.23.1. I've also checked the C++ source files and they do contain the fix proposed above @ghost above, so I think this version is patched. Thanks a lot for your help!

R version 4.3.2 (2023-10-31)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Debian GNU/Linux 12 (bookworm)

Matrix products: default
BLAS:   /usr/lib/x86_64-linux-gnu/openblas-pthread/libblas.so.3
LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/libopenblasp-r0.3.21.so;  LAPACK version 3.11.0

locale:
 [1] LC_CTYPE=en_GB.UTF-8       LC_NUMERIC=C               LC_TIME=en_GB.UTF-8
 [4] LC_COLLATE=en_GB.UTF-8     LC_MONETARY=en_GB.UTF-8    LC_MESSAGES=en_GB.UTF-8
 [7] LC_PAPER=en_GB.UTF-8       LC_NAME=C                  LC_ADDRESS=C
[10] LC_TELEPHONE=C             LC_MEASUREMENT=en_GB.UTF-8 LC_IDENTIFICATION=C

time zone: Europe/Amsterdam
tzcode source: system (glibc)

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base

other attached packages:
[1] rTRNG_4.23.1-2

loaded via a namespace (and not attached):
[1] compiler_4.3.2     Rcpp_1.0.12        codetools_0.2-19   RcppParallel_5.1.7
oelhammouchi commented

FYI: I've found a temporary fix by generating a uniform sample from TRNG and passing this to the inverse CDF from R (which can be accessed in C++ via the Rcpp package) and it seems to work fine.

rabauke commented

@OthmanElHammouchi I identified very slow numerical convergence in the function GammaP_ser for some extreme parameters as the root-cause and opened a new ticket. See #31.

  • ❤️1