From bf63dc4d583df4f67a6c91c4151cd6f36964f66f Mon Sep 17 00:00:00 2001 From: cgarling Date: Sun, 28 Sep 2025 12:39:38 -0400 Subject: [PATCH 1/6] try using Chernoff sampler --- src/univariate/continuous/chernoff.jl | 157 ++++++++++++++------------ 1 file changed, 84 insertions(+), 73 deletions(-) diff --git a/src/univariate/continuous/chernoff.jl b/src/univariate/continuous/chernoff.jl index 3d319db252..62266985bb 100644 --- a/src/univariate/continuous/chernoff.jl +++ b/src/univariate/continuous/chernoff.jl @@ -149,75 +149,10 @@ 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=[ - 0.0 -Inf; - 0.01 -1.171534341573129; - 0.025 -0.9981810946684274; - 0.05 -0.8450811886357725; - 0.1 -0.6642351964332931; - 0.2 -0.43982766604886553; - 0.25 -0.353308035220429; - 0.3 -0.2751512847290148; - 0.4 -0.13319637678583637; - 0.5 0.0; - 0.6 0.13319637678583637; - 0.7 0.2751512847290147; - 0.75 0.353308035220429; - 0.8 0.4398276660488655; - 0.9 0.6642351964332931; - 0.95 0.8450811886357724; - 0.975 0.9981810946684272; - 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 +186,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 +220,98 @@ 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) + +@quantile_newton Chernoff + +# function quantile(d::Chernoff, tau::Real) +# # some commonly used quantiles were precomputed +# precomputedquants=[ +# 0.0 -Inf; +# 0.01 -1.171534341573129; +# 0.025 -0.9981810946684274; +# 0.05 -0.8450811886357725; +# 0.1 -0.6642351964332931; +# 0.2 -0.43982766604886553; +# 0.25 -0.353308035220429; +# 0.3 -0.2751512847290148; +# 0.4 -0.13319637678583637; +# 0.5 0.0; +# 0.6 0.13319637678583637; +# 0.7 0.2751512847290147; +# 0.75 0.353308035220429; +# 0.8 0.4398276660488655; +# 0.9 0.6642351964332931; +# 0.95 0.8450811886357724; +# 0.975 0.9981810946684272; +# 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 + +struct ChernoffSampler <: Sampleable{Univariate, Continuous} +end +ChernoffSampler(::Chernoff) = ChernoffSampler() +function rand(rng::AbstractRNG, d::ChernoffSampler) + 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) + return quantile(Chernoff(), tauabs + (1-tauabs)*F0) * sign(tau) end +sampler(d::Chernoff) = ChernoffSampler(d) +rand(rng::AbstractRNG, d::Chernoff) = rand(rng, sampler(d)) \ No newline at end of file From 801cbfc3d11ce0d470f6d137d9076e59b0608a73 Mon Sep 17 00:00:00 2001 From: cgarling Date: Sun, 28 Sep 2025 12:40:16 -0400 Subject: [PATCH 2/6] Simplify, remove sampler --- src/univariate/continuous/chernoff.jl | 7 +------ 1 file changed, 1 insertion(+), 6 deletions(-) diff --git a/src/univariate/continuous/chernoff.jl b/src/univariate/continuous/chernoff.jl index 62266985bb..61403aa5eb 100644 --- a/src/univariate/continuous/chernoff.jl +++ b/src/univariate/continuous/chernoff.jl @@ -287,10 +287,7 @@ entropy(d::Chernoff) = -0.7515605300273104 ### Random number generation -struct ChernoffSampler <: Sampleable{Univariate, Continuous} -end -ChernoffSampler(::Chernoff) = ChernoffSampler() -function rand(rng::AbstractRNG, d::ChernoffSampler) +function rand(rng::AbstractRNG, d::Chernoff) A = ChernoffComputations.randA x = ChernoffComputations.randx y = ChernoffComputations.randy @@ -313,5 +310,3 @@ function rand(rng::AbstractRNG, d::ChernoffSampler) tauabs = abs(tau) return quantile(Chernoff(), tauabs + (1-tauabs)*F0) * sign(tau) end -sampler(d::Chernoff) = ChernoffSampler(d) -rand(rng::AbstractRNG, d::Chernoff) = rand(rng, sampler(d)) \ No newline at end of file From b8368123cc8cc2951faa14f4ce8178d5a5f9449b Mon Sep 17 00:00:00 2001 From: cgarling Date: Sun, 28 Sep 2025 12:40:25 -0400 Subject: [PATCH 3/6] Add more Chernoff tests --- test/univariate/continuous/chernoff.jl | 12 ++++++++++++ 1 file changed, 12 insertions(+) diff --git a/test/univariate/continuous/chernoff.jl b/test/univariate/continuous/chernoff.jl index a5f28d02bc..760c273e63 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; @@ -154,5 +165,6 @@ 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 end From 88ed54670d930e621f7d01efda5ca1644e4b0b3b Mon Sep 17 00:00:00 2001 From: cgarling Date: Sun, 28 Sep 2025 13:06:39 -0400 Subject: [PATCH 4/6] Fix quantile and add tests --- src/univariate/continuous/chernoff.jl | 81 ++++++++++++-------------- test/univariate/continuous/chernoff.jl | 6 ++ 2 files changed, 44 insertions(+), 43 deletions(-) diff --git a/src/univariate/continuous/chernoff.jl b/src/univariate/continuous/chernoff.jl index 61403aa5eb..37248bbc79 100644 --- a/src/univariate/continuous/chernoff.jl +++ b/src/univariate/continuous/chernoff.jl @@ -150,6 +150,29 @@ module ChernoffComputations _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] + # some commonly used quantiles were precomputed + const precomputedquants=[ + 0.0 -Inf; + 0.01 -1.171534341573129; + 0.025 -0.9981810946684274; + 0.05 -0.8450811886357725; + 0.1 -0.6642351964332931; + 0.2 -0.43982766604886553; + 0.25 -0.353308035220429; + 0.3 -0.2751512847290148; + 0.4 -0.13319637678583637; + 0.5 0.0; + 0.6 0.13319637678583637; + 0.7 0.2751512847290147; + 0.75 0.353308035220429; + 0.8 0.4398276660488655; + 0.9 0.6642351964332931; + 0.95 0.8450811886357724; + 0.975 0.9981810946684272; + 0.99 1.17153434157313; + 1.0 Inf + ] + # constants needed for the Ziggurat random sampling algorithm const randA = 0.03248227216266608 const randx = [ @@ -226,51 +249,23 @@ 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) -@quantile_newton Chernoff - -# function quantile(d::Chernoff, tau::Real) -# # some commonly used quantiles were precomputed -# precomputedquants=[ -# 0.0 -Inf; -# 0.01 -1.171534341573129; -# 0.025 -0.9981810946684274; -# 0.05 -0.8450811886357725; -# 0.1 -0.6642351964332931; -# 0.2 -0.43982766604886553; -# 0.25 -0.353308035220429; -# 0.3 -0.2751512847290148; -# 0.4 -0.13319637678583637; -# 0.5 0.0; -# 0.6 0.13319637678583637; -# 0.7 0.2751512847290147; -# 0.75 0.353308035220429; -# 0.8 0.4398276660488655; -# 0.9 0.6642351964332931; -# 0.95 0.8450811886357724; -# 0.975 0.9981810946684272; -# 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 +# @quantile_newton Chernoff -# # 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) +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 -# 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 + # 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 diff --git a/test/univariate/continuous/chernoff.jl b/test/univariate/continuous/chernoff.jl index 760c273e63..8e218590f5 100644 --- a/test/univariate/continuous/chernoff.jl +++ b/test/univariate/continuous/chernoff.jl @@ -159,6 +159,8 @@ 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 @@ -167,4 +169,8 @@ @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 From 9fa8a17b9e1231e453eddc7d9c5a5a919280a48c Mon Sep 17 00:00:00 2001 From: cgarling Date: Sun, 28 Sep 2025 13:10:58 -0400 Subject: [PATCH 5/6] remove leftover code --- src/univariate/continuous/chernoff.jl | 2 -- 1 file changed, 2 deletions(-) diff --git a/src/univariate/continuous/chernoff.jl b/src/univariate/continuous/chernoff.jl index 37248bbc79..77915f6fce 100644 --- a/src/univariate/continuous/chernoff.jl +++ b/src/univariate/continuous/chernoff.jl @@ -249,8 +249,6 @@ 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) -# @quantile_newton Chernoff - function quantile(d::Chernoff, tau::Real) precomputedquants = ChernoffComputations.precomputedquants (0 <= tau && tau <= 1) || throw(DomainError(tau, "illegal value of tau")) From 8934077261123cd35b6aefa8cd1292643e5a52a0 Mon Sep 17 00:00:00 2001 From: cgarling Date: Sun, 28 Sep 2025 13:14:50 -0400 Subject: [PATCH 6/6] revert change in rand --- src/univariate/continuous/chernoff.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/univariate/continuous/chernoff.jl b/src/univariate/continuous/chernoff.jl index 77915f6fce..2d89c3546e 100644 --- a/src/univariate/continuous/chernoff.jl +++ b/src/univariate/continuous/chernoff.jl @@ -301,5 +301,5 @@ function rand(rng::AbstractRNG, d::Chernoff) F0 = ChernoffComputations._cdf(A/y[1]) tau = 2*rand(rng)-1 # ~ U(-1,1) tauabs = abs(tau) - return quantile(Chernoff(), tauabs + (1-tauabs)*F0) * sign(tau) + return quantile(d, tauabs + (1-tauabs)*F0) * sign(tau) end