Skip to content

Chernoff distribution fix #2000

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
wants to merge 6 commits into
base: master
Choose a base branch
from
Open

Chernoff distribution fix #2000

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
6 commits
Select commit Hold shift + click to select a range
bf63dc4
try using Chernoff sampler
cgarling Sep 28, 2025
801cbfc
Simplify, remove sampler
cgarling Sep 28, 2025
b836812
Add more Chernoff tests
cgarling Sep 28, 2025
88ed546
Fix quantile and add tests
cgarling Sep 28, 2025
9fa8a17
remove leftover code
cgarling Sep 28, 2025
8934077
revert change in rand
cgarling Sep 28, 2025
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
99 changes: 49 additions & 50 deletions src/univariate/continuous/chernoff.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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;
Expand All @@ -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
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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
18 changes: 18 additions & 0 deletions test/univariate/continuous/chernoff.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -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;
Expand Down Expand Up @@ -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
Loading