Merge pull request #305 from JuliaStats/dh/probs

Add probs methods to discrete distributions

src/testutils.jl

@@ -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.
     #

src/univariate/bernoulli.jl

@@ -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) 

src/univariate/binomial.jl

@@ -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

src/univariate/categorical.jl

@@ -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
 

src/univariate/discreteuniform.jl

@@ -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

src/univariate/geometric.jl

@@ -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)

src/univariate/hypergeometric.jl

@@ -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)
 

src/univariate/negativebinomial.jl

@@ -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

src/univariate/poisson.jl

@@ -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

test/discrete_ref.csv

@@ -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

test/discrete_ref.txt

@@ -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)