dgpd and pgpd return wrong value with xi->0

[avehtari]

dgpd and pgpd are not handling xi=0 special case correctly

dgpd(1) [1] 0 dgpd(1, xi=0) [1] 0 dgpd(1, xi=-1e-10) [1] 0.3678794 dgpd(1, xi=1e-10) [1] 0.3678794 pgpd(1) [1] 0 pgpd(1, xi=0) [1] 0 pgpd(1, xi=-1e-10) [1] 0.6321206 pgpd(1, xi=1e-10) [1] 0.6321206

Also cases with xi close to 0 return non-smooth results

dgpd(1, xi=1e-16) [1] 1 dgpd(1, xi=1e-15) [1] 0.3294855 dgpd(1, xi=1e-14) [1] 0.3681736 dgpd(1, xi=1e-13) [1] 0.3681736 dgpd(1, xi=1e-12) [1] 0.3678467 dgpd(1, xi=1e-11) [1] 0.3678794 dgpd(1, xi=1e-10) [1] 0.3678794

The correct results with xi->0 can be obtained from dexp and pexp

dexp(1) [1] 0.3678794 pexp(1) [1] 0.6321206

[twolodzko]

Thanks, I'll look at it.

[avehtari]

Here's a link to a Stan code which uses numerically stable way to compute these https://github.com/avehtari/BDA_R_demos/blob/master/demos_rstan/gpareto_functions/gpareto.stan Even if you don't know Stan the code should be easily readable.

And thanks for making the extraDistr package, it's beeen useful.