diff --git a/src/univariate/continuous/chernoff.jl b/src/univariate/continuous/chernoff.jl index 3d319db25..2d89c3546 100644 --- a/src/univariate/continuous/chernoff.jl +++ b/src/univariate/continuous/chernoff.jl @@ -149,15 +149,9 @@ module ChernoffComputations _pdf(x::Real) = g(x)*g(-x)*0.5 _cdf(x::Real) = (x < 0.0) ? _cdfbar(-x) : 0.5 + quadgk(_pdf,0.0,x)[1] _cdfbar(x::Real) = (x < 0.0) ? _cdf(x) : quadgk(_pdf, x, Inf)[1] -end - -pdf(d::Chernoff, x::Real) = ChernoffComputations._pdf(x) -logpdf(d::Chernoff, x::Real) = log(ChernoffComputations.g(x))+log(ChernoffComputations.g(-x))+log(0.5) -cdf(d::Chernoff, x::Real) = ChernoffComputations._cdf(x) -function quantile(d::Chernoff, tau::Real) # some commonly used quantiles were precomputed - precomputedquants=[ + const precomputedquants=[ 0.0 -Inf; 0.01 -1.171534341573129; 0.025 -0.9981810946684274; @@ -178,46 +172,10 @@ function quantile(d::Chernoff, tau::Real) 0.99 1.17153434157313; 1.0 Inf ] - (0.0 <= tau && tau <= 1.0) || throw(DomainError(tau, "illegal value of tau")) - present = searchsortedfirst(precomputedquants[:, 1], tau) - if present <= size(precomputedquants, 1) - if tau == precomputedquants[present, 1] - return precomputedquants[present, 2] - end - end - - # one good approximation of the quantiles can be computed using Normal(0.0, stdapprox) with stdapprox = 0.52 - stdapprox = 0.52 - dnorm = Normal(0.0, 1.0) - if tau < 0.001 - return -newton(x -> tau - ChernoffComputations._cdfbar(x), ChernoffComputations._pdf, quantile(dnorm, 1.0 - tau)*stdapprox) - - end - if tau > 0.999 - return newton(x -> 1.0 - tau - ChernoffComputations._cdfbar(x), ChernoffComputations._pdf, quantile(dnorm, tau)*stdapprox) - end - return newton(x -> ChernoffComputations._cdf(x) - tau, ChernoffComputations._pdf, quantile(dnorm, tau)*stdapprox) # should consider replacing x-> construct for speed -end -minimum(d::Chernoff) = -Inf -maximum(d::Chernoff) = Inf -insupport(d::Chernoff, x::Real) = isnan(x) ? false : true - -mean(d::Chernoff) = 0.0 -var(d::Chernoff) = 0.26355964132470455 -modes(d::Chernoff) = [0.0] -mode(d::Chernoff) = 0.0 -skewness(d::Chernoff) = 0.0 -kurtosis(d::Chernoff) = -0.16172525511461888 -kurtosis(d::Chernoff, excess::Bool) = kurtosis(d) + (excess ? 0.0 : 3.0) -entropy(d::Chernoff) = -0.7515605300273104 - -### Random number generation -rand(d::Chernoff) = rand(default_rng(), d) -function rand(rng::AbstractRNG, d::Chernoff) # Ziggurat random number generator --- slow in the tails - # constants needed for the Ziggurat algorithm - A = 0.03248227216266608 - x = [ + # constants needed for the Ziggurat random sampling algorithm + const randA = 0.03248227216266608 + const randx = [ 1.4765521793744492 1.3583996502410562 1.2788224934376338 @@ -251,7 +209,7 @@ function rand(rng::AbstractRNG, d::Chernoff) # Ziggurat random n 0.2358977457249061 1.0218214689661219e-7 ] - y=[ + const randy=[ 0.02016386420423385 0.042162593823411566 0.06607475557706186 @@ -285,22 +243,63 @@ function rand(rng::AbstractRNG, d::Chernoff) # Ziggurat random n 1.3789927085182485 1.516689116183566 ] +end + +pdf(d::Chernoff, x::Real) = ChernoffComputations._pdf(x) +logpdf(d::Chernoff, x::Real) = log(ChernoffComputations.g(x))+log(ChernoffComputations.g(-x))+log(0.5) +cdf(d::Chernoff, x::Real) = ChernoffComputations._cdf(x) + +function quantile(d::Chernoff, tau::Real) + precomputedquants = ChernoffComputations.precomputedquants + (0 <= tau && tau <= 1) || throw(DomainError(tau, "illegal value of tau")) + present = searchsortedfirst(precomputedquants[:, 1], tau) + if present <= size(precomputedquants, 1) + if tau == precomputedquants[present, 1] + return precomputedquants[present, 2] + end + end + + # one good approximation of the quantiles can be computed using Normal(0.0, stdapprox) with stdapprox = 0.52 + stdapprox = 0.52 + dnorm = Normal(0.0, 1.0) + return quantile_newton(d, tau, quantile(dnorm, tau)*stdapprox) +end + +minimum(d::Chernoff) = -Inf +maximum(d::Chernoff) = Inf +insupport(d::Chernoff, x::Real) = isnan(x) ? false : true + +mean(d::Chernoff) = 0.0 +var(d::Chernoff) = 0.26355964132470455 +modes(d::Chernoff) = [0.0] +mode(d::Chernoff) = 0.0 +skewness(d::Chernoff) = 0.0 +kurtosis(d::Chernoff) = -0.16172525511461888 +kurtosis(d::Chernoff, excess::Bool) = kurtosis(d) + (excess ? 0.0 : 3.0) +entropy(d::Chernoff) = -0.7515605300273104 + +### Random number generation + +function rand(rng::AbstractRNG, d::Chernoff) + A = ChernoffComputations.randA + x = ChernoffComputations.randx + y = ChernoffComputations.randy n = length(x) i = rand(rng, 0:n-1) - r = (2.0*rand(rng)-1) * ((i>0) ? x[i] : A/y[1]) + r = (2*rand(rng)-1) * ((i>0) ? x[i] : A/y[1]) rabs = abs(r) if rabs < x[i+1] return r end s = (i>0) ? (y[i]+rand(rng)*(y[i+1]-y[i])) : rand(rng)*y[1] - if s < 2.0*ChernoffComputations._pdf(rabs) + if s < 2*ChernoffComputations._pdf(rabs) return r end if i > 0 return rand(rng, d) end F0 = ChernoffComputations._cdf(A/y[1]) - tau = 2.0*rand(rng)-1 # ~ U(-1,1) + tau = 2*rand(rng)-1 # ~ U(-1,1) tauabs = abs(tau) return quantile(d, tauabs + (1-tauabs)*F0) * sign(tau) end diff --git a/test/univariate/continuous/chernoff.jl b/test/univariate/continuous/chernoff.jl index a5f28d02b..8e218590f 100644 --- a/test/univariate/continuous/chernoff.jl +++ b/test/univariate/continuous/chernoff.jl @@ -1,5 +1,16 @@ @testset "Chernoff tests" begin d = Chernoff() + @test extrema(d) == (-Inf, Inf) + @test mean(d) == 0 + @test var(d) ≈ 0.26355964132470455 + @test modes(d) == [0.0] + @test mode(d) == 0 + @test skewness(d) == 0 + @test kurtosis(d) ≈ -0.16172525511461888 + @test kurtosis(d, true) ≈ kurtosis(d) + @test kurtosis(d, false) ≈ kurtosis(d) + 3 + @test entropy(d) ≈ -0.7515605300273104 + @test rand(d) isa Float64 cdftest = [ 0.005 0.5037916689930134; @@ -148,11 +159,18 @@ 2.000 0.0001212 ] + quantiletest = [0.0 -Inf; 0.05 -0.8450811888163325; 0.1 -0.664235196540868; 0.15 -0.5398553654712455; 0.2 -0.439827666120028; 0.25 -0.35330803528117866; 0.3 -0.27515128478220097; 0.35 -0.20241833430160983; 0.4 -0.13319637683484448; 0.45 -0.06609664086333447; 0.5 0.0; 0.55 0.06609664081686001; 0.6 0.1331963767858364; 0.65 0.20241833424985095; 0.7 0.2751512847290145; 0.75 0.3533080352204289; 0.8 0.4398276660488657; 0.85 0.5398553653947077; 0.9 0.6642351964332935; 0.95 0.845081188635776; 1.0 Inf] + for i=1:size(cdftest, 1) @test isapprox(cdf(d, cdftest[i, 1]), cdftest[i, 2]) end for i=1:size(pdftest, 1) @test isapprox(pdf(d, pdftest[i, 1]), pdftest[i, 2] ; atol = 1e-6) + @test isapprox(logpdf(d, pdftest[i, 1]), log(pdftest[i, 2]) ; atol = 1e-4) + end + + for i=1:size(quantiletest, 1) + @test isapprox(quantile(d, quantiletest[i,1]), quantiletest[i, 2]) end end