test the computation of Davies' terms

JuliaStats · Nov 29, 2023 · 3ab5b06 · 3ab5b06
1 parent 251c061
commit 3ab5b06
Showing 1 changed file with 12 additions and 6 deletions.
diff --git a/test/univariate/continuous/generalizedchisq.jl b/test/univariate/continuous/generalizedchisq.jl
@@ -3,8 +3,14 @@ using Test
 using Random
 
 @testset "Generalized chi-squared" begin
-    rng = Random.MersenneTwister(0)
     Computations = Distributions.GChisqComputations
+    # Function obtain the cf from Davies' terms, given:
+    # exp(θ*im)/ρ = exp(-u*x*im)*cf(u)
+    function cfdavies(d, u, x=randn())
+        θ, ρ = Computations.daviesterms(d, u, x)
+        return exp(θ*im + u*x*im)/ρ
+    end
+    rng = Random.MersenneTwister(0)
     quadgk = Distributions.quadgk
     # General case
     gx2 = GeneralizedChisq([1,-1], [1,2], [1.5, 1.5], 10, 1)
@@ -13,7 +19,7 @@ using Random
     @test maximum(gx2) == Inf
     @test insupport(gx2, rand(gx2))
     @test let t = randn()
-        cf(gx2, t) ≈ Computations.cf_inherit(gx2, t)
+        cf(gx2, t) ≈ Computations.cf_inherit(gx2, t) ≈ cfdavies(gx2, t)
     end
     @test cdf(gx2, 15) - cdf(gx2, 5) ≈ first(quadgk(x->pdf(gx2, x), 5, 15))
     @test pdf(gx2, 15) - pdf(gx2, 5) ≈ first(quadgk(x->Computations.daviesdpdf(gx2, x), 5, 15))
@@ -34,7 +40,7 @@ using Random
     @test all(≥(10), rand(gx2, 10))
     @test minimum(gx2) == 10
     @test let t = randn()
-        cf(gx2, t) ≈ Computations.cf_inherit(gx2, t)
+        cf(gx2, t) ≈ Computations.cf_inherit(gx2, t) ≈ cfdavies(gx2, t)
     end
     @test cdf(gx2, 15) - cdf(gx2, 5) ≈ first(quadgk(x->pdf(gx2, x), 5, 15)) ≈ first(quadgk(x->pdf(gx2, x), 10, 15))
     @testset for p in 0:0.05:1
@@ -46,7 +52,7 @@ using Random
     @test all(≤(10), rand(gx2, 10))
     @test maximum(gx2) == 10
     @test let t = randn()
-        cf(gx2, t) ≈ Computations.cf_inherit(gx2, t)
+        cf(gx2, t) ≈ Computations.cf_inherit(gx2, t) ≈ cfdavies(gx2, t)
     end
     @test cdf(gx2, 15) - cdf(gx2, 5) ≈ first(quadgk(x->pdf(gx2, x), 5, 15)) ≈ first(quadgk(x->pdf(gx2, x), 5, 10))
     @testset for p in 0:0.05:1
@@ -58,7 +64,7 @@ using Random
     gx2b = GeneralizedChisq(Int[], Int[], Int[], 10, 1)
     normalequiv = Normal(gx2.μ, gx2.σ)
     @test let t = randn()
-        cf(gx2, t) ≈ cf(gx2b, t) ≈ Computations.cf_inherit(gx2, t) ≈ cf(normalequiv, t)
+        cf(gx2, t) ≈ cf(gx2b, t) ≈ Computations.cf_inherit(gx2, t) ≈ cfdavies(gx2, t) ≈ cf(normalequiv, t)
     end
     let x = 10 + randn()
         @test cdf(gx2, x) ≈ cdf(gx2b, x) ≈ Computations.daviescdf(gx2, x) ≈ cdf(normalequiv, x)
@@ -73,7 +79,7 @@ using Random
     gx2 = GeneralizedChisq([0,0], [1,2], [0,1], 10, 0)
     gx2b = GeneralizedChisq(Int[], Int[], Int[], 10, 0)
     @test let t = randn()
-        cf(gx2, t) ≈ cf(gx2b, t) ≈ Computations.cf_inherit(gx2, t)
+        cf(gx2, t) ≈ cf(gx2b, t) ≈ Computations.cf_inherit(gx2, t) ≈ cfdavies(gx2, t)
     end
     let x = 9
         @test cdf(gx2, x) == cdf(gx2b, x) == 0