Ce serveur Gitlab sera éteint le 30 juin 2020, pensez à migrer vos projets vers les serveurs gitlab-research.centralesupelec.fr et gitlab-student.centralesupelec.fr !

Commit 0792898c authored by Dahua Lin's avatar Dahua Lin

Merge pull request #305 from JuliaStats/dh/probs

Add probs methods to discrete distributions
parents 67d5eae4 1efccbee
......@@ -81,6 +81,15 @@ function test_samples(s::Sampleable{Univariate, Discrete}, # the sampleable
p0 = pdf(distr, rmin:rmax) # reference probability masses
@assert length(p0) == m
# check the consistency between probs and pdf
if isa(s, Distribution)
@test_approx_eq probs(s, rmin:rmax) p0
if isbounded(s)
@assert isfinite(vmin) && isfinite(vmax)
@test_approx_eq probs(s) probs(s, vmin:vmax)
end
end
# determine confidence intervals for counts:
# with probability q, the count will be out of this interval.
#
......
......@@ -18,6 +18,26 @@ end
Bernoulli() = Bernoulli(0.5)
probs(d::Bernoulli) = [d.p0, d.p1]
function probs(d::Bernoulli, rgn::UnitRange)
n = length(rgn)
if n == 0
return Float64[]
elseif n == 1
f = rgn[1]
return f == zero(f) ? [d.p0] :
f == one(f) ? [d.p1] :
throw(BoundsError())
elseif n == 2
f = rgn[1]
return f == zero(f) ? [d.p0, d.p1] :
throw(BoundsError())
else
throw(BoundsError())
end
end
cdf(d::Bernoulli, q::Real) = q >= zero(q) ? (q >= one(q) ? 1.0 : d.p0) : 0.
function entropy(d::Bernoulli)
......
......@@ -16,6 +16,31 @@ maximum(d::Binomial) = d.size
@_jl_dist_2p Binomial binom
function _probs(d::Binomial, f::Int, l::Int)
n = d.size
p = d.prob
b = f - 1
r = Array(Float64, l - b)
r[1] = v = pdf(d, f)
if l > f
c = p / (1.0 - p)
for k = f+1:l
v *= ((n - k + 1) / k * c)
r[k-b] = v
end
end
return r
end
probs(d::Binomial) = _probs(d, 0, d.size)
function probs(d::Binomial, rgn::UnitRange)
f, l = rgn[1], rgn[end]
0 <= f <= l <= d.size || throw(BoundsError())
_probs(d, f, l)
end
function entropy(d::Binomial; approx::Bool=false)
n = d.size
p1 = d.prob
......
......@@ -14,6 +14,7 @@ end
ncategories(d::Categorical) = d.K
probs(d::Categorical) = d.prob
probs(d::Categorical, rgn::UnitRange) = d.prob[rgn]
### handling support
......
......@@ -11,6 +11,17 @@ end
DiscreteUniform(b::Integer) = DiscreteUniform(0, b)
DiscreteUniform() = DiscreteUniform(0, 1)
function probs(d::DiscreteUniform)
n = d.b - d.a + 1
fill(1.0 / n, n)
end
function probs(d::DiscreteUniform, rgn::UnitRange)
f, l = rgn[1], rgn[end]
d.a <= f <= l <= d.b || throw(BoundsError())
fill(1.0 / (d.b - d.a + 1), l - f + 1)
end
function cdf(d::DiscreteUniform, k::Real)
k < d.a ? 0. : (k > d.b ? 1. : (ifloor(k) - d.a + 1.0) / (d.b - d.a + 1.0))
end
......
......@@ -8,6 +8,20 @@ end
Geometric() = Geometric(0.5) # Flips of a fair coin
function probs(d::Geometric, rgn::UnitRange)
f, l = rgn[1], rgn[end]
0 <= f <= l || throw(BoundsError())
p = d.prob
r = Array(Float64, l - f + 1)
pfail = 1.0 - p
r[1] = v = (pfail^f) * p
b = f - 1
for i = f+1:l
r[i - b] = (v *= pfail)
end
return r
end
## Support
insupport(::Geometric, x::Real) = isinteger(x) && x >= zero(x)
insupport(::Type{Geometric}, x::Real) = isinteger(x) && x >= zero(x)
......
......@@ -26,6 +26,31 @@ support(d::Hypergeometric) = minimum(d):maximum(d)
islowerbounded(d::Hypergeometric) = true
isupperbounded(d::Hypergeometric) = true
function _probs(d::Hypergeometric, f::Int, l::Int)
ns = d.ns
nf = d.nf
n = d.n
r = Array(Float64, l - f + 1)
r[1] = v = pdf(d, f)
b = f - 1
if l > f
for x = f+1:l
c = ((ns - x + 1) / x) * ((n - x + 1) / (nf - n + x))
r[x-b] = (v *= c)
end
end
return r
end
probs(d::Hypergeometric) = _probs(d, minimum(d), maximum(d))
function probs(d::Hypergeometric, rgn::UnitRange)
f, l = rgn[1], rgn[end]
minimum(d) <= f <= l <= maximum(d) || throw(BoundsError())
_probs(d, f, l)
end
# properties
mean(d::Hypergeometric) = d.n * d.ns / (d.ns + d.nf)
......
......@@ -29,6 +29,21 @@ maximum(::Union(NegativeBinomial, Type{NegativeBinomial})) = Inf
insupport(::NegativeBinomial, x::Real) = isinteger(x) && zero(x) <= x
insupport(::Type{NegativeBinomial}, x::Real) = isinteger(x) && zero(x) <= x
function probs(d::NegativeBinomial, rgn::UnitRange)
r = d.r
p0 = 1.0 - d.prob
f, l = rgn[1], rgn[end]
0 <= f <= l || throw(BoundsError())
res = Array(Float64, l - f + 1)
res[1] = v = pdf(d, f)
b = f - 1
for x = f+1:l
c = (x + r - 1) * p0 / x
res[x-b] = (v *= c)
end
return res
end
function mgf(d::NegativeBinomial, t::Real)
r, p = d.r, d.prob
return ((1.0 - p) * exp(t))^r / (1.0 - p * exp(t))^r
......
......@@ -9,21 +9,6 @@ end
@_jl_dist_1p Poisson pois
function entropy(d::Poisson)
λ = d.lambda
if λ < 50.0
s = 0.0
for k in 1:100
s += λ^k * lgamma(k + 1.0) / gamma(k + 1.0)
end
return λ * (1.0 - log(λ)) + exp(-λ) * s
else
return 0.5 * log(2 * pi * e * λ) -
(1 / (12 * λ)) -
(1 / (24 * λ * λ)) -
(19 / (360 * λ * λ * λ))
end
end
isupperbounded(::Union(Poisson, Type{Poisson})) = false
islowerbounded(::Union(Poisson, Type{Poisson})) = true
......@@ -35,12 +20,52 @@ maximum(::Union(Poisson, Type{Poisson})) = Inf
insupport(::Poisson, x::Real) = isinteger(x) && zero(x) <= x
insupport(::Type{Poisson}, x::Real) = isinteger(x) && zero(x) <= x
kurtosis(d::Poisson) = 1.0 / d.lambda
function probs(d::Poisson, rgn::UnitRange)
λ = d.lambda
f, l = rgn[1], rgn[end]
0 <= f <= l || throw(BoundsError())
r = Array(Float64, l - f + 1)
v = r[1] = pdf(d, f)
if l > f
b = f - 1
for x = f+1:l
c = λ / x
r[x - b] = (v *= c)
end
end
return r
end
mean(d::Poisson) = d.lambda
median(d::Poisson) = quantile(d, 0.5)
mode(d::Poisson) = ifloor(d.lambda)
modes(d::Poisson) = isinteger(d.lambda) ? [int(d.lambda)-1,int(d.lambda)] : [ifloor(d.lambda)]
var(d::Poisson) = d.lambda
skewness(d::Poisson) = 1.0 / sqrt(d.lambda)
kurtosis(d::Poisson) = 1.0 / d.lambda
function entropy(d::Poisson)
λ = d.lambda
if λ < 50.0
s = 0.0
for k in 1:100
s += λ^k * lgamma(k + 1.0) / gamma(k + 1.0)
end
return λ * (1.0 - log(λ)) + exp(-λ) * s
else
return 0.5 * log(2 * pi * e * λ) -
(1 / (12 * λ)) -
(1 / (24 * λ * λ)) -
(19 / (360 * λ * λ * λ))
end
end
function mgf(d::Poisson, t::Real)
l = d.lambda
return exp(l * (exp(t) - 1.0))
......@@ -51,13 +76,6 @@ function cf(d::Poisson, t::Real)
return exp(l * (exp(im * t) - 1.0))
end
mode(d::Poisson) = ifloor(d.lambda)
modes(d::Poisson) = isinteger(d.lambda) ? [int(d.lambda)-1,int(d.lambda)] : [ifloor(d.lambda)]
skewness(d::Poisson) = 1.0 / sqrt(d.lambda)
var(d::Poisson) = d.lambda
# model fitting
immutable PoissonStats <: SufficientStats
......
......@@ -2,6 +2,8 @@
"Bernoulli(0.9)", 9.0000000000000002e-01, 8.9999999999999983e-02, 3.2508297339144820e-01, 1, 1, 1, -1.0536051565782628e-01, -1.0536051565782628e-01, -1.0536051565782628e-01
"Bernoulli(0.1)", 1.0000000000000001e-01, 9.0000000000000011e-02, 3.2508297339144820e-01, 0, 0, 0, -1.0536051565782631e-01, -1.0536051565782631e-01, -1.0536051565782631e-01
"Binomial(1, 0.5)", 5.0000000000000000e-01, 2.5000000000000000e-01, 6.9314718055994529e-01, 0, 0, 1, -6.9314718055994529e-01, -6.9314718055994529e-01, -6.9314718055994529e-01
"Binomial(5, 0.4)", 2.0000000000000000e+00, 1.2000000000000000e+00, 1.4979981829038540e+00, 1, 2, 3, -1.3501553145040179e+00, -1.0624732420522367e+00, -1.4679383501604009e+00
"Binomial(6, 0.8)", 4.8000000000000007e+00, 9.5999999999999996e-01, 1.3425774508710622e+00, 4, 5, 6, -1.4033998290228298e+00, -9.3339619977709365e-01, -1.3388613078852583e+00
"Binomial(100, 0.1)", 1.0000000000000000e+01, 9.0000000000000000e+00, 2.5112331161825550e+00, 8, 10, 12, -2.1643632354294020e+00, -2.0259739768661866e+00, -2.3147790140625677e+00
"Binomial(100, 0.9)", 9.0000000000000000e+01, 8.9999999999999982e+00, 2.5112331161825514e+00, 88, 90, 92, -2.3147790140625695e+00, -2.0259739768661902e+00, -2.1643632354294020e+00
"DiscreteUniform(0, 4)", 2.0000000000000000e+00, 2.0000000000000000e+00, 1.6094379124341003e+00, 1, 2, 3, -1.6094379124341003e+00, -1.6094379124341003e+00, -1.6094379124341003e+00
......
......@@ -2,6 +2,8 @@ Bernoulli(0.5)
Bernoulli(0.9)
Bernoulli(0.1)
Binomial(1, 0.5)
Binomial(5, 0.4)
Binomial(6, 0.8)
Binomial(100, 0.1)
Binomial(100, 0.9)
DiscreteUniform(0, 4)
......
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment