Error in Gamma_dist with large kappa about trng4 HOT 10 CLOSED

@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.

from trng4.

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?

from trng4.

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

Regards,
Heiko

from trng4.

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

from trng4.

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

Regards,
Heiko

from trng4.

This has been fixed in version 4.16.

from trng4.

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?

from trng4.

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.

from trng4.

@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

from trng4.

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.

from trng4.