diff --git a/src/samplers.jl b/src/samplers.jl index fa667fff2b..f6df74b828 100644 --- a/src/samplers.jl +++ b/src/samplers.jl @@ -14,19 +14,20 @@ Samples from the sampler and returns the result. """ rand(::AbstractRNG, ::Sampleable) -for fname in ["aliastable.jl", - "binomial.jl", - "poissonbinomial.jl", - "poisson.jl", - "exponential.jl", - "gamma.jl", - "multinomial.jl", - "vonmises.jl", - "vonmisesfisher.jl", - "discretenonparametric.jl", - "categorical.jl", - "productnamedtuple.jl", - ] - +for fname in [ + "aliastable.jl" + "binomial.jl" + "poissonbinomial.jl" + "poisson.jl" + "exponential.jl" + "gamma.jl" + "multinomial.jl" + "vonmises.jl" + "vonmisesfisher.jl" + "discretenonparametric.jl" + "categorical.jl" + "productnamedtuple.jl" + "hypergeometric.jl" +] include(joinpath("samplers", fname)) end diff --git a/src/samplers/hypergeometric.jl b/src/samplers/hypergeometric.jl new file mode 100644 index 0000000000..a54ba16011 --- /dev/null +++ b/src/samplers/hypergeometric.jl @@ -0,0 +1,195 @@ +""" +### Hypergeometric Sampler + +***Sampling implementation reference:*** V. Kachitvichyanukul & B. Schmeiser + "Computer generation of hypergeometric random variates" + Journal of Statistical Computation and Simulation, 22(2):127-145 + doi:10.1080/00949658508810839 +""" +struct HypergeometricSampler{T} <: Sampleable{Univariate,Discrete} + + dist::Hypergeometric + ns_opt::Int + nf_opt::Int + n_opt::Int + N::Int + cache::T + + # some addapted names between this implementation and the original paper (see docstring) + # | paper | here | + # |-------|--------| + # | n1* | ns | successes in population ─┬─ stored in the distribution struct + # | n2* | nf | failures in population ──┤ + # | k* | n | sample size ─────────────╯ + # | n1 | ns_opt | ─┬─ optimized parameters for the algorithms + # | n2 | nf_opt | ─┤ + # | k | n_opt | ─╯ + # | n | N | total population size + # | x | y | random variate being generated (specifically on HIN) + + function HypergeometricSampler(dist::Hypergeometric) + + # Step 0: Set-up constants + N = dist.ns + dist.nf + ns_opt, nf_opt = if dist.ns > dist.nf + dist.nf, dist.ns + else + dist.ns, dist.nf + end + n_opt = if dist.n > N / 2 + N - dist.n + else + dist.n + end + M = floor(Int, (n_opt + 1) * (ns_opt + 1) / (N + 2)) + cache = if M - max(0, n_opt - nf_opt) < 10 + # for the HIN algorithm + if n_opt < nf_opt + p = exp( + logabsgamma(nf_opt + 1)[1] + + logabsgamma(N - n_opt + 1)[1] - + logabsgamma(N + 1)[1] - + logabsgamma(nf_opt - n_opt + 1)[1] + ) + y = 0 + HINCache(p, y) + else + p = exp( + logabsgamma(ns_opt + 1)[1] + + logabsgamma(n_opt + 1)[1] - + logabsgamma(n_opt - nf_opt + 1)[1] - + logabsgamma(N + 1)[1] + ) + y = n_opt - nf_opt + HINCache(p, y) + end + else + # for the H2PE algorithm + A = logabsgamma(M + 1)[1] + + logabsgamma(ns_opt - M + 1)[1] + + logabsgamma(n_opt - M + 1)[1] + + logabsgamma(nf_opt - n_opt + M + 1)[1] + D = 1.5 * sqrt((N - n_opt) * n_opt * ns_opt * nf_opt / ((N - 1) * N * N)) + 0.5 + xL = M - D + 0.5 + xR = M + D + 0.5 + kL = exp( + A - + logabsgamma(xL + 1)[1] - + logabsgamma(ns_opt - xL + 1)[1] - + logabsgamma(n_opt - xL + 1)[1] - + logabsgamma(nf_opt - n_opt + xL + 1)[1] + ) + kR = exp( + A - + logabsgamma(xR)[1] - + logabsgamma(ns_opt - xR + 2)[1] - + logabsgamma(n_opt - xR + 2)[1] - + logabsgamma(nf_opt - n_opt + xR)[1] + ) + λL = -log(xL * (nf_opt - n_opt + xL) / ((ns_opt - xL + 1) * (n_opt - xL + 1))) + λR = -log((ns_opt - xR + 1) * (n_opt - xR + 1) / (xR * (nf_opt - n_opt + xR))) + p1 = 2D + p2 = p1 + kL / λL + p3 = p2 + kR / λR + H2PECache(A, xL, xR, λL, λR, p1, p2, p3) + end + new{typeof(cache)}(dist, ns_opt, nf_opt, n_opt, N, cache) + end + +end + +struct H2PECache + A::Float64 + xL::Float64 + xR::Float64 + λL::Float64 + λR::Float64 + p1::Float64 + p2::Float64 + p3::Float64 +end + +struct HINCache + p::Float64 + y::Int +end + +function Random.rand(rng::AbstractRNG, spl::HypergeometricSampler{H2PECache}) + # Steps 1, 2, 3, 4: base variate generation + cache = spl.cache + y::Int = 0 + while true + u = rand(rng) * cache.p3 + v = rand(rng) + if u <= cache.p1 + # Region 1: central region + y = floor(Int, cache.xL + u) + elseif u <= cache.p2 + # Region 2: left exponential tail + y = floor(Int, cache.xL + log(v) / cache.λL) + if y < max(0, spl.n_opt - spl.nf_opt) + continue + end + v = v * (u - cache.p1) * cache.λL + else + # Region 3: right exponential tail + y = floor(Int, cache.xR - log(v) / cache.λR) + if y > min(spl.ns_opt, spl.n_opt) + continue + end + v = v * (u - cache.p2) * cache.λR + end + # Step 4: Acceptance/Regection Comparison + # note: seems like optimizations found in the paper don't improve the + # performance compared to this direct implementation with the lgamma function, + # so the check is done directly for simplicity. + if log(v) > cache.A - + logabsgamma(y + 1)[1] - + logabsgamma(spl.ns_opt - y + 1)[1] - + logabsgamma(spl.n_opt - y + 1)[1] - + logabsgamma(spl.nf_opt - spl.n_opt + y + 1)[1] + continue + else + break + end + end + + # Step 5: return appropriate random variate + return _hg_correct_variate(spl, y) + +end + +function Random.rand(rng::AbstractRNG, spl::HypergeometricSampler{HINCache}) + cache = spl.cache + u = rand(rng) + p = cache.p + y = cache.y + while u > p + u -= p + p *= (spl.ns_opt - y) * + (spl.n_opt - y) / + (y + 1) / + (spl.nf_opt - spl.n_opt + 1 + y) + y += 1 + end + return _hg_correct_variate(spl, y) +end + +"***Internal:*** Addapts the sampled variable from the optimized distribution for output" +function _hg_correct_variate(spl::HypergeometricSampler, y::Integer) + dist = spl.dist + if dist.n < spl.N / 2 + if dist.ns <= dist.nf + y + else + dist.n - y + end + else + if dist.ns <= dist.nf + dist.ns - y + else + dist.n - dist.nf + y + end + end +end + diff --git a/src/univariate/discrete/hypergeometric.jl b/src/univariate/discrete/hypergeometric.jl index a13127f2f6..bdf33dc842 100644 --- a/src/univariate/discrete/hypergeometric.jl +++ b/src/univariate/discrete/hypergeometric.jl @@ -1,7 +1,8 @@ """ Hypergeometric(s, f, n) -A *Hypergeometric distribution* describes the number of successes in `n` draws without replacement from a finite population containing `s` successes and `f` failures. +A *Hypergeometric distribution* describes the number of successes in `n` draws without +replacement from a finite population containing `s` successes and `f` failures. ```math P(X = k) = {{{s \\choose k} {f \\choose {n-k}}}\\over {s+f \\choose n}}, \\quad \\text{for } k = \\max(0, n - f), \\ldots, \\min(n, s). @@ -17,12 +18,11 @@ params(d) # Get the parameters, i.e. (s, f, n) External links * [Hypergeometric distribution on Wikipedia](http://en.wikipedia.org/wiki/Hypergeometric_distribution) - """ struct Hypergeometric <: DiscreteUnivariateDistribution - ns::Int # number of successes in population - nf::Int # number of failures in population - n::Int # sample size + ns::Int # number of successes in population + nf::Int # number of failures in population + n::Int # sample size function Hypergeometric(ns::Real, nf::Real, n::Real; check_args::Bool=true) @check_args( @@ -58,22 +58,27 @@ mode(d::Hypergeometric) = floor(Int, (d.n + 1) * (d.ns + 1) / (d.ns + d.nf + 2)) function modes(d::Hypergeometric) if (d.ns == d.nf) && mod(d.n, 2) == 1 - [(d.n-1)/2, (d.n+1)/2] + [(d.n - 1) / 2, (d.n + 1) / 2] else [mode(d)] end end -skewness(d::Hypergeometric) = (d.nf-d.ns)*sqrt(d.ns+d.nf-1)*(d.ns+d.nf-2*d.n)/sqrt(d.n*d.ns*d.nf*(d.ns+d.nf-d.n))/(d.ns+d.nf-2) +skewness(d::Hypergeometric) = (d.nf - d.ns) * + sqrt(d.ns + d.nf - 1) * + (d.ns + d.nf - 2 * d.n) / + sqrt(d.n * d.ns * d.nf * (d.ns + d.nf - d.n)) / + (d.ns + d.nf - 2) function kurtosis(d::Hypergeometric) ns = float(d.ns) nf = float(d.nf) n = float(d.n) N = ns + nf - a = (N-1) * N^2 * (N * (N+1) - 6*ns * (N-ns) - 6*n*(N-n)) + 6*n*ns*(nf)*(N-n)*(5*N-6) - b = (n*ns*(N-ns) * (N-n)*(N-2)*(N-3)) - a/b + a = (N - 1) * N^2 * (N * (N + 1) - 6 * ns * (N - ns) - 6 * n * (N - n)) + + 6 * n * ns * (nf) * (N - n) * (5 * N - 6) + b = (n * ns * (N - ns) * (N - n) * (N - 2) * (N - 3)) + a / b end entropy(d::Hypergeometric) = entropy(map(Base.Fix1(pdf, d), support(d))) @@ -84,11 +89,4 @@ entropy(d::Hypergeometric) = entropy(map(Base.Fix1(pdf, d), support(d))) ## sampling -# TODO: remove RFunctions dependency. Implement: -# V. Kachitvichyanukul & B. Schmeiser -# "Computer generation of hypergeometric random variates" -# Journal of Statistical Computation and Simulation, 22(2):127-145 -# doi:10.1080/00949658508810839 -@rand_rdist(Hypergeometric) -rand(d::Hypergeometric) = - convert(Int, StatsFuns.RFunctions.hyperrand(d.ns, d.nf, d.n)) +sampler(dist::Hypergeometric) = HypergeometricSampler(dist) diff --git a/test/ref/discrete_test.lst b/test/ref/discrete_test.lst index f130624af1..01d997cc3d 100644 --- a/test/ref/discrete_test.lst +++ b/test/ref/discrete_test.lst @@ -36,6 +36,7 @@ Hypergeometric(3, 2, 0) Hypergeometric(3, 2, 5) Hypergeometric(4, 5, 6) Hypergeometric(60, 80, 100) +Hypergeometric(200, 200, 10) NegativeBinomial() NegativeBinomial(6) diff --git a/test/ref/discrete_test.ref.json b/test/ref/discrete_test.ref.json index 05a465be71..7b3e5cda76 100644 --- a/test/ref/discrete_test.ref.json +++ b/test/ref/discrete_test.ref.json @@ -19,7 +19,7 @@ { "x": 1, "pdf": 0.5, "logpdf": -0.693147180559945, "cdf": 1 } ], "quans": [ - { "q": 0.10, "x": 0 }, + { "q": 0.10, "x": -0 }, { "q": 0.25, "x": 0 }, { "q": 0.50, "x": 0 }, { "q": 0.75, "x": 1 }, @@ -46,8 +46,8 @@ { "x": 1, "pdf": 0.25, "logpdf": -1.38629436111989, "cdf": 1 } ], "quans": [ - { "q": 0.10, "x": 0 }, - { "q": 0.25, "x": 0 }, + { "q": 0.10, "x": -0 }, + { "q": 0.25, "x": -0 }, { "q": 0.50, "x": 0 }, { "q": 0.75, "x": 0 }, { "q": 0.90, "x": 1 } @@ -173,10 +173,10 @@ }, "points": [ { "x": 0, "pdf": 0.343519287158209, "logpdf": -1.06851201996383, "cdf": 0.343519287158209 }, - { "x": 1, "pdf": 0.0742744404666397, "logpdf": -2.59998839092822, "cdf": 0.417793727624848 }, - { "x": 3, "pdf": 0.039339841571296, "logpdf": -3.23551749319517, "cdf": 0.50574956659249 }, - { "x": 6, "pdf": 0.0345020592656908, "logpdf": -3.36673626786553, "cdf": 0.609338682032028 }, - { "x": 9, "pdf": 0.0657942390366793, "logpdf": -2.72122299711787, "cdf": 0.757877200345019 }, + { "x": 1, "pdf": 0.0742744404666398, "logpdf": -2.59998839092822, "cdf": 0.417793727624849 }, + { "x": 3, "pdf": 0.0393398415712961, "logpdf": -3.23551749319517, "cdf": 0.505749566592491 }, + { "x": 6, "pdf": 0.0345020592656909, "logpdf": -3.36673626786553, "cdf": 0.609338682032029 }, + { "x": 9, "pdf": 0.0657942390366794, "logpdf": -2.72122299711787, "cdf": 0.757877200345021 }, { "x": 10, "pdf": 0.24212279965498, "logpdf": -1.41831024493703, "cdf": 1 } ], "quans": [ @@ -200,16 +200,16 @@ "kurtosis": -0.811058355437666 }, "points": [ - { "x": 0, "pdf": 0.0518518518518519, "logpdf": -2.95936462938312, "cdf": 0.0518518518518519 }, + { "x": 0, "pdf": 0.0518518518518519, "logpdf": -2.95936462938312, "cdf": 0.0518518518518518 }, { "x": 1, "pdf": 0.0901771336553946, "logpdf": -2.40597939119833, "cdf": 0.142028985507246 }, { "x": 2, "pdf": 0.115942028985507, "logpdf": -2.15466496291742, "cdf": 0.257971014492754 }, - { "x": 3, "pdf": 0.1301805237732, "logpdf": -2.03883314739231, "cdf": 0.388151538265954 }, - { "x": 4, "pdf": 0.134009362707706, "logpdf": -2.00984561051905, "cdf": 0.522160900973661 }, - { "x": 5, "pdf": 0.128648988199398, "logpdf": -2.05066760503931, "cdf": 0.650809889173059 }, - { "x": 6, "pdf": 0.115454220178947, "logpdf": -2.15888118967954, "cdf": 0.766264109352006 }, - { "x": 7, "pdf": 0.0959619492396443, "logpdf": -2.34380352817355, "cdf": 0.86222605859165 }, - { "x": 8, "pdf": 0.0719714619297336, "logpdf": -2.63148560062533, "cdf": 0.934197520521384 }, - { "x": 10, "pdf": 0.0201063131740208, "logpdf": -3.90672142497276, "cdf": 0.999999999999997 } + { "x": 3, "pdf": 0.1301805237732, "logpdf": -2.03883314739231, "cdf": 0.388151538265955 }, + { "x": 4, "pdf": 0.134009362707706, "logpdf": -2.00984561051905, "cdf": 0.522160900973662 }, + { "x": 5, "pdf": 0.128648988199398, "logpdf": -2.05066760503931, "cdf": 0.650809889173061 }, + { "x": 6, "pdf": 0.115454220178947, "logpdf": -2.15888118967954, "cdf": 0.766264109352009 }, + { "x": 7, "pdf": 0.0959619492396443, "logpdf": -2.34380352817355, "cdf": 0.862226058591654 }, + { "x": 8, "pdf": 0.0719714619297336, "logpdf": -2.63148560062533, "cdf": 0.934197520521388 }, + { "x": 10, "pdf": 0.0201063131740208, "logpdf": -3.90672142497276, "cdf": 1 } ], "quans": [ { "q": 0.10, "x": 1 }, @@ -232,13 +232,13 @@ "kurtosis": -0.212506273471436 }, "points": [ - { "x": 0, "pdf": 0.000192751819501661, "logpdf": -8.55410710586806, "cdf": 0.000192751819501661 }, - { "x": 4, "pdf": 0.113798908976598, "logpdf": -2.17332234453623, "cdf": 0.177330155587457 }, - { "x": 5, "pdf": 0.194216804653393, "logpdf": -1.63878019415293, "cdf": 0.371546960240851 }, - { "x": 6, "pdf": 0.239092657243757, "logpdf": -1.43090411496952, "cdf": 0.610639617484608 }, - { "x": 7, "pdf": 0.209702529941374, "logpdf": -1.56206527657206, "cdf": 0.820342147425982 }, - { "x": 8, "pdf": 0.125447049161356, "logpdf": -2.0758715284762, "cdf": 0.945789196587338 }, - { "x": 10, "pdf": 0.00797557662147816, "logpdf": -4.83137132929139, "cdf": 0.999999999999993 } + { "x": 0, "pdf": 0.000192751819501661, "logpdf": -8.55410710586806, "cdf": 0.000192751819501639 }, + { "x": 4, "pdf": 0.113798908976598, "logpdf": -2.17332234453623, "cdf": 0.177330155587428 }, + { "x": 5, "pdf": 0.194216804653393, "logpdf": -1.63878019415293, "cdf": 0.371546960240797 }, + { "x": 6, "pdf": 0.239092657243757, "logpdf": -1.43090411496952, "cdf": 0.610639617484534 }, + { "x": 7, "pdf": 0.209702529941374, "logpdf": -1.56206527657206, "cdf": 0.820342147425872 }, + { "x": 8, "pdf": 0.125447049161356, "logpdf": -2.0758715284762, "cdf": 0.945789196587212 }, + { "x": 10, "pdf": 0.00797557662147816, "logpdf": -4.83137132929139, "cdf": 1 } ], "quans": [ { "q": 0.10, "x": 4 }, @@ -268,7 +268,7 @@ { "x": 1, "pdf": 0.5, "logpdf": -0.693147180559945, "cdf": 1 } ], "quans": [ - { "q": 0.10, "x": 0 }, + { "q": 0.10, "x": -0 }, { "q": 0.25, "x": 0 }, { "q": 0.50, "x": 0 }, { "q": 0.75, "x": 1 }, @@ -351,7 +351,7 @@ "kurtosis": 0.0416666666666667 }, "points": [ - { "x": 0, "pdf": 6.4e-05, "logpdf": -9.6566274746046, "cdf": 6.39999999999999e-05 }, + { "x": 0, "pdf": 6.39999999999999e-05, "logpdf": -9.6566274746046, "cdf": 6.39999999999999e-05 }, { "x": 1, "pdf": 0.001536, "logpdf": -6.47857364425666, "cdf": 0.0016 }, { "x": 2, "pdf": 0.01536, "logpdf": -4.17598855126261, "cdf": 0.01696 }, { "x": 3, "pdf": 0.0819199999999999, "logpdf": -2.50201211769094, "cdf": 0.09888 }, @@ -390,9 +390,9 @@ { "x": 9, "pdf": 0.130416277079164, "logpdf": -2.03702381305276, "cdf": 0.451290165442003 }, { "x": 10, "pdf": 0.131865346824488, "logpdf": -2.02597397686618, "cdf": 0.583155512266491 }, { "x": 11, "pdf": 0.119877588022262, "logpdf": -2.1212841566705, "cdf": 0.703033100288754 }, - { "x": 12, "pdf": 0.0987880123516788, "logpdf": -2.31477901406258, "cdf": 0.801821112640432 }, - { "x": 14, "pdf": 0.0513038273344495, "logpdf": -2.96998992267902, "cdf": 0.927427034735119 }, - { "x": 100, "pdf": 1.00000000000002e-100, "logpdf": -230.258509299405, "cdf": 1 } + { "x": 12, "pdf": 0.0987880123516788, "logpdf": -2.31477901406258, "cdf": 0.801821112640433 }, + { "x": 14, "pdf": 0.0513038273344494, "logpdf": -2.96998992267902, "cdf": 0.927427034735119 }, + { "x": 100, "pdf": 1.00000000000001e-100, "logpdf": -230.258509299405, "cdf": 1 } ], "quans": [ { "q": 0.10, "x": 6 }, @@ -418,7 +418,7 @@ "kurtosis": 0.0511111111111111 }, "points": [ - { "x": 0, "pdf": 9.99999999999989e-101, "logpdf": -230.258509299405, "cdf": 9.99999999999998e-101 }, + { "x": 0, "pdf": 9.99999999999978e-101, "logpdf": -230.258509299405, "cdf": 9.99999999999998e-101 }, { "x": 86, "pdf": 0.0513038273344493, "logpdf": -2.96998992267903, "cdf": 0.12387679259933 }, { "x": 88, "pdf": 0.0987880123516788, "logpdf": -2.31477901406258, "cdf": 0.296966899711246 }, { "x": 89, "pdf": 0.119877588022262, "logpdf": -2.1212841566705, "cdf": 0.416844487733508 }, @@ -732,7 +732,7 @@ { "x": 4, "pdf": 0.06561, "logpdf": -2.72402715562535, "cdf": 0.40951 }, { "x": 6, "pdf": 0.0531441, "logpdf": -2.934748186941, "cdf": 0.5217031 }, { "x": 8, "pdf": 0.043046721, "logpdf": -3.14546921825666, "cdf": 0.612579511 }, - { "x": 11, "pdf": 0.031381059609, "logpdf": -3.46155076523013, "cdf": 0.717570463519 }, + { "x": 11, "pdf": 0.031381059609, "logpdf": -3.46155076523014, "cdf": 0.717570463519 }, { "x": 15, "pdf": 0.0205891132094649, "logpdf": -3.88299282786144, "cdf": 0.814697981114816 }, { "x": 21, "pdf": 0.0109418989131512, "logpdf": -4.5151559218084, "cdf": 0.901522909781639 }, { "x": 43, "pdf": 0.00107752636643058, "logpdf": -6.83308726628058, "cdf": 0.990302262702125 } @@ -931,15 +931,15 @@ "kurtosis": -0.02787878186819 }, "points": [ - { "x": 20, "pdf": 2.37106739517318e-20, "logpdf": -45.1883616284133, "cdf": 2.37106739517318e-20 }, - { "x": 39, "pdf": 0.0523495571361668, "logpdf": -2.94981180137345, "cdf": 0.102455999861417 }, + { "x": 20, "pdf": 2.37106739517317e-20, "logpdf": -45.1883616284133, "cdf": 2.37106739517317e-20 }, + { "x": 39, "pdf": 0.0523495571361666, "logpdf": -2.94981180137346, "cdf": 0.102455999861417 }, { "x": 41, "pdf": 0.116828889706323, "logpdf": -2.14704489579628, "cdf": 0.303109617932027 }, { "x": 42, "pdf": 0.141737213593927, "logpdf": -1.95378054436581, "cdf": 0.444846831525954 }, { "x": 43, "pdf": 0.14961946510825, "logpdf": -1.89966010754594, "cdf": 0.594466296634204 }, { "x": 44, "pdf": 0.137292861448764, "logpdf": -1.98563895992138, "cdf": 0.731759158082968 }, { "x": 45, "pdf": 0.109346136762749, "logpdf": -2.21323686158497, "cdf": 0.841105294845717 }, { "x": 46, "pdf": 0.0754268919977156, "logpdf": -2.58459140976086, "cdf": 0.916532186843432 }, - { "x": 60, "pdf": 6.08101843998881e-13, "logpdf": -28.1284340203951, "cdf": 1 } + { "x": 60, "pdf": 6.08101843998883e-13, "logpdf": -28.1284340203951, "cdf": 1 } ], "quans": [ { "q": 0.10, "x": 39 }, @@ -949,6 +949,34 @@ { "q": 0.90, "x": 46 } ] }, +{ + "expr": "Hypergeometric(200, 200, 10)", + "dtype": "Hypergeometric", + "minimum": 0, + "maximum": 10, + "properties": { + "mean": 5, + "var": 2.44360902255639, + "skewness": 0, + "kurtosis": -0.191048246463864 + }, + "points": [ + { "x": 0, "pdf": 0.000870258877244979, "logpdf": -7.04671983052616, "cdf": 0.000870258877244979 }, + { "x": 3, "pdf": 0.116274924318445, "logpdf": -2.15179785475634, "cdf": 0.168759849917107 }, + { "x": 4, "pdf": 0.206627732777236, "logpdf": -1.57683649714626, "cdf": 0.375387582694344 }, + { "x": 5, "pdf": 0.249224834611313, "logpdf": -1.38939983968554, "cdf": 0.624612417305657 }, + { "x": 6, "pdf": 0.206627732777236, "logpdf": -1.57683649714626, "cdf": 0.831240150082893 }, + { "x": 7, "pdf": 0.116274924318445, "logpdf": -2.15179785475634, "cdf": 0.947515074401337 }, + { "x": 10, "pdf": 0.000870258877244979, "logpdf": -7.04671983052616, "cdf": 1 } + ], + "quans": [ + { "q": 0.10, "x": 3 }, + { "q": 0.25, "x": 4 }, + { "q": 0.50, "x": 5 }, + { "q": 0.75, "x": 6 }, + { "q": 0.90, "x": 7 } + ] +}, { "expr": "NegativeBinomial()", "dtype": "NegativeBinomial", @@ -1002,7 +1030,7 @@ { "x": 7, "pdf": 0.0966796875, "logpdf": -2.33635195546486, "cdf": 0.70947265625 }, { "x": 9, "pdf": 0.06109619140625, "logpdf": -2.79530574852401, "cdf": 0.84912109375 }, { "x": 11, "pdf": 0.0333251953125, "logpdf": -3.40144155209433, "cdf": 0.928268432617188 }, - { "x": 16, "pdf": 0.00485157966613769, "logpdf": -5.32845092270461, "cdf": 0.991549730300903 } + { "x": 16, "pdf": 0.0048515796661377, "logpdf": -5.32845092270461, "cdf": 0.991549730300903 } ], "quans": [ { "q": 0.10, "x": 2 }, @@ -1062,7 +1090,7 @@ { "x": 3, "pdf": 0.1741824, "logpdf": -1.74765225296301, "cdf": 0.5940864 }, { "x": 4, "pdf": 0.13934592, "logpdf": -1.97079580427721, "cdf": 0.73343232 }, { "x": 5, "pdf": 0.1003290624, "logpdf": -2.29929987124925, "cdf": 0.8337613824 }, - { "x": 6, "pdf": 0.0668860416, "logpdf": -2.70476497935742, "cdf": 0.900647424 }, + { "x": 6, "pdf": 0.0668860416, "logpdf": -2.70476497935741, "cdf": 0.900647424 }, { "x": 10, "pdf": 0.008161880702976, "logpdf": -4.80828065825628, "cdf": 0.990652339224576 } ], "quans": [ @@ -1118,7 +1146,7 @@ { "x": 0, "pdf": 0.367879441171442, "logpdf": -1, "cdf": 0.367879441171442 }, { "x": 1, "pdf": 0.367879441171442, "logpdf": -1, "cdf": 0.735758882342885 }, { "x": 2, "pdf": 0.183939720585721, "logpdf": -1.69314718055995, "cdf": 0.919698602928606 }, - { "x": 3, "pdf": 0.0613132401952404, "logpdf": -2.79175946922806, "cdf": 0.981011843123846 }, + { "x": 3, "pdf": 0.0613132401952404, "logpdf": -2.79175946922805, "cdf": 0.981011843123846 }, { "x": 4, "pdf": 0.0153283100488101, "logpdf": -4.17805383034795, "cdf": 0.996340153172656 } ], "quans": [ @@ -1228,8 +1256,8 @@ { "x": 10, "pdf": 0.125110035721133, "logpdf": -2.07856164313506, "cdf": 0.583039750192985 }, { "x": 11, "pdf": 0.113736396110121, "logpdf": -2.17387182293938, "cdf": 0.696776146303107 }, { "x": 12, "pdf": 0.0947803300917677, "logpdf": -2.35619337973334, "cdf": 0.791556476394874 }, - { "x": 13, "pdf": 0.0729079462244366, "logpdf": -2.61855764420083, "cdf": 0.864464422619311 }, - { "x": 14, "pdf": 0.0520771044460263, "logpdf": -2.95502988082204, "cdf": 0.916541527065337 }, + { "x": 13, "pdf": 0.0729079462244367, "logpdf": -2.61855764420083, "cdf": 0.864464422619311 }, + { "x": 14, "pdf": 0.0520771044460262, "logpdf": -2.95502988082204, "cdf": 0.916541527065337 }, { "x": 18, "pdf": 0.00709110899319529, "logpdf": -4.94891353414023, "cdf": 0.992813495396146 } ], "quans": [ @@ -1259,11 +1287,11 @@ { "x": 75, "pdf": 0.0392264266855184, "logpdf": -3.23840460918207, "cdf": 0.312384231977445 }, { "x": 78, "pdf": 0.0439997074482216, "logpdf": -3.12357229399004, "cdf": 0.440574629646123 }, { "x": 80, "pdf": 0.0445566657703509, "logpdf": -3.11099351178318, "cdf": 0.529687961186824 }, - { "x": 82, "pdf": 0.0429332521725754, "logpdf": -3.1481086443721, "cdf": 0.61662779683629 }, + { "x": 82, "pdf": 0.0429332521725755, "logpdf": -3.1481086443721, "cdf": 0.61662779683629 }, { "x": 85, "pdf": 0.0370926143436918, "logpdf": -3.29433740348069, "cdf": 0.734512761797335 }, - { "x": 87, "pdf": 0.0317285126703592, "logpdf": -3.45053954904101, "cdf": 0.80074603199671 }, + { "x": 87, "pdf": 0.0317285126703592, "logpdf": -3.45053954904102, "cdf": 0.80074603199671 }, { "x": 92, "pdf": 0.0176179446534293, "logpdf": -4.03883731377811, "cdf": 0.916442472359937 }, - { "x": 102, "pdf": 0.00244734344970641, "logpdf": -6.01275214895309, "cdf": 0.992398851472148 } + { "x": 102, "pdf": 0.00244734344970639, "logpdf": -6.0127521489531, "cdf": 0.992398851472148 } ], "quans": [ { "q": 0.10, "x": 69 }, @@ -1344,7 +1372,7 @@ }, "points": [ { "x": -6, "pdf": 0.0154316649323744, "logpdf": -4.17133371613553, "cdf": 0.0244462081254431 }, - { "x": -4, "pdf": 0.068410944672505, "logpdf": -2.6822224573219, "cdf": 0.127731053628502 }, + { "x": -4, "pdf": 0.068410944672505, "logpdf": -2.6822224573219, "cdf": 0.127731053628501 }, { "x": -3, "pdf": 0.114463860117042, "logpdf": -2.16749613797573, "cdf": 0.242194913745544 }, { "x": -2, "pdf": 0.160071156565379, "logpdf": -1.83213683407712, "cdf": 0.402266070310923 }, { "x": -1, "pdf": 0.183023344454947, "logpdf": -1.6981415689489, "cdf": 0.58528941476587 }, diff --git a/test/runtests.jl b/test/runtests.jl index f63a3d5e8e..9b05725b0d 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -100,6 +100,7 @@ const tests = [ "univariate/continuous/triangular", "statsapi", "univariate/continuous/inversegaussian", + "univariate/discrete/hypergeometric" ### missing files compared to /src: # "common", diff --git a/test/truncate.jl b/test/truncate.jl index 9c2a286d09..5525995106 100644 --- a/test/truncate.jl +++ b/test/truncate.jl @@ -78,12 +78,15 @@ function verify_and_test(d::UnivariateDistribution, dct::Dict, n_tsamples::Int) end @test cdf(d, x) ≈ cf atol=sqrt(eps()) # NOTE: some distributions use pdf() in StatsFuns.jl which have no generic support yet - if !any(T -> d isa T, [Distributions.Truncated{Distributions.NoncentralChisq{Float64},Distributions.Continuous, Float64}, - Distributions.Truncated{Distributions.NoncentralF{Float64},Distributions.Continuous, Float64}, - Distributions.Truncated{Distributions.NoncentralT{Float64},Distributions.Continuous, Float64}, - Distributions.Truncated{Distributions.StudentizedRange{Float64},Distributions.Continuous, Float64}, - Distributions.Truncated{Distributions.Rician{Float64},Distributions.Continuous, Float64}]) - @test isapprox(logpdf(d, Dual(float(x))), lp, atol=sqrt(eps())) + if !(d isa Truncated{T} where {T<:Union{ + NoncentralChisq, + NoncentralF, + NoncentralT, + StudentizedRange, + Rician, + Hypergeometric, + }}) + @test logpdf(d, Dual(float(x))) ≈ lp atol = sqrt(eps()) end # NOTE: this test is disabled as StatsFuns.jl doesn't have generic support for cdf() # @test isapprox(cdf(d, Dual(x)) , cf, atol=sqrt(eps())) diff --git a/test/univariate/discrete/hypergeometric.jl b/test/univariate/discrete/hypergeometric.jl new file mode 100644 index 0000000000..9a72b6c58b --- /dev/null +++ b/test/univariate/discrete/hypergeometric.jl @@ -0,0 +1,9 @@ +using Distributions +using Test + +@testset "Hypergeometric unimodal or bimodal" begin + # bimodal if ns == nf and n is odd, else unimodal + for ns ∈ 1:10, nf ∈ 1:10, n ∈ 0:ns+nf + @test length(modes(Hypergeometric(ns, nf, n))) == (ns == nf && isodd(n) ? 2 : 1) + end +end \ No newline at end of file