fix numerical edge cases in Devroye algorithm

wzhorton · May 10, 2023 · de0320c · de0320c
1 parent 61e0bbc
commit de0320c
Show file tree

Hide file tree

Showing 2 changed files with 3 additions and 3 deletions.
diff --git a/Project.toml b/Project.toml
@@ -1,7 +1,7 @@
 name = "PolyaGammaHybridSamplers"
 uuid = "c636ee4f-4591-4d8c-9fae-2dea21daa433"
 authors = ["W.Z. Horton"]
-version = "1.0.3"
+version = "1.0.4"
 
 [deps]
 Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"

diff --git a/src/rand_pgdevroye.jl b/src/rand_pgdevroye.jl
@@ -37,7 +37,7 @@ function rand_Jstar(z::Real, rng::AbstractRNG)
     t = 0.64 # paper recommends this constant
     K = 0.125*π^2 + 0.5*z^2
     p = 0.5*π*inv(K) * exp(-t*K)
-    q = 2*exp(-z) * cdf(InverseGaussian(inv(z), 1.0), t)
+    q = 2*exp(-z + logcdf(InverseGaussian(inv(z), 1.0), t))
     while true
         # Generate X
         if rand(rng, Uniform(0.0, 1.0)) < p/(p+q) #U ~ Uniform(0,1)
@@ -60,7 +60,7 @@ function rand_Jstar(z::Real, rng::AbstractRNG)
             n += 1
             if isodd(n)
                 S -= a_coefs(n, X, t)
-                if Y < S
+                if Y <= S
                     return X # accept X and exit
                 end
             else