Merge pull request JuliaStats#386 from chrodan/master

Fix for Issue JuliaStats#384 (pdf of Binomial)

src/univariate/discrete/binomial.jl

@@ -4,8 +4,8 @@ immutable Binomial <: DiscreteUnivariateDistribution
     p::Float64
 
     function Binomial(n::Int, p::Float64)
-        n >= 0 || error("n must be non-negative.")
-        0.0 <= p <= 1.0 || error("p must be in [0, 1]")
+        n >= 0 || error("n must be non-negative but is $n.")
+        0.0 <= p <= 1.0 || error("p must be in [0, 1] but is $p")
         new(n, p)
     end
 
@@ -76,9 +76,37 @@ immutable RecursiveBinomProbEvaluator <: RecursiveProbabilityEvaluator
 end
 
 RecursiveBinomProbEvaluator(d::Binomial) = RecursiveBinomProbEvaluator(d.n, d.p / (1.0 - d.p))
-nextpdf(s::RecursiveBinomProbEvaluator, p::Float64, x::Integer) = ((s.n - x + 1) / x) * s.coef * p
-_pdf!(r::AbstractArray, d::Binomial, rgn::UnitRange) = _pdf!(r, d, rgn, RecursiveBinomProbEvaluator(d))
-
+nextpdf(s::RecursiveBinomProbEvaluator, pv::Float64, x::Integer) = ((s.n - x + 1) / x) * s.coef * pv
+
+function _pdf!(r::AbstractArray, d::Binomial, X::UnitRange)
+    vl,vr, vfirst, vlast = _pdf_fill_outside!(r, d, X)
+    if succprob(d) <= 0.5
+        # fill normal
+        rpe = RecursiveBinomProbEvaluator(d::Binomial)
+
+        # fill central part: with non-zero pdf
+        if vl <= vr
+            fm1 = vfirst - 1
+            r[vl - fm1] = pv = pdf(d, vl)
+            for v = (vl+1):vr
+                r[v - fm1] = pv = nextpdf(rpe, pv, v)
+            end
+        end
+    else
+        # fill reversed to avoid 1/0 for d.p==1.
+        rpe = RecursiveBinomProbEvaluator(d.n, (1.0 - d.p) / d.p)
+
+        # fill central part: with non-zero pdf
+        if vl <= vr
+            fm1 = vfirst - 1
+            r[vr - fm1] = pv = pdf(d, vr)
+            for v = (vr-1):-1:vl
+                r[v - fm1] = pv = nextpdf(rpe, pv, d.n-v)
+            end
+        end
+    end
+    return r
+end
 
 function mgf(d::Binomial, t::Real)
     n, p = params(d)

src/univariates.jl

@@ -316,7 +316,8 @@ for fun in [:pdf, :logpdf,
     end
 end
 
-function _pdf!(r::AbstractArray, d::DiscreteUnivariateDistribution, X::UnitRange)
+
+function _pdf_fill_outside!(r::AbstractArray, d::DiscreteUnivariateDistribution, X::UnitRange)
     vl = vfirst = first(X)
     vr = vlast = last(X)
     n = vlast - vfirst + 1
@@ -352,35 +353,25 @@ function _pdf!(r::AbstractArray, d::DiscreteUnivariateDistribution, X::UnitRange
             r[i] = 0.0
         end
     end
+    return vl, vr, vfirst, vlast
+end
+
+function _pdf!(r::AbstractArray, d::DiscreteUnivariateDistribution, X::UnitRange)
+    vl,vr, vfirst, vlast = _pdf_fill_outside!(r, d, X)
+
+    # fill central part: with non-zero pdf
+    fm1 = vfirst - 1
+    for v = vl:vr
+        r[v - fm1] = pdf(d, v)
+    end
     return r
 end
 
 
 abstract RecursiveProbabilityEvaluator
 
 function _pdf!(r::AbstractArray, d::DiscreteUnivariateDistribution, X::UnitRange, rpe::RecursiveProbabilityEvaluator)
-    vl = vfirst = first(X)
-    vr = vlast = last(X)
-    n = vlast - vfirst + 1
-    if islowerbounded(d)
-        lb = minimum(d)
-        if vl < lb
-            vl = lb
-        end
-    end
-    if isupperbounded(d)
-        ub = maximum(d)
-        if vr > ub
-            vr = ub
-        end
-    end
-
-    # fill left part
-    if vl > vfirst
-        for i = 1:(vl - vfirst)
-            r[i] = 0.0
-        end
-    end
+    vl,vr, vfirst, vlast = _pdf_fill_outside!(r, d, X)
 
     # fill central part: with non-zero pdf
     if vl <= vr
@@ -391,12 +382,6 @@ function _pdf!(r::AbstractArray, d::DiscreteUnivariateDistribution, X::UnitRange
         end
     end
 
-    # fill right part
-    if vr < vlast
-        for i = (vr-vfirst+2):n
-            r[i] = 0.0
-        end
-    end
     return r
 end
 

test/binomial.jl

@@ -0,0 +1,23 @@
+using Distributions
+using Base.Test
+
+# Test the consistency between the recursive and nonrecursive computation of the pdf 
+# of the Binomial distribution
+srand(1234)
+for (p, n) in [(0.6, 10), (0.8, 6), (0.5, 40), (0.04, 20), (1., 100), (0., 10), (0.999999, 1000), (1e-7, 1000)]
+
+    d = Binomial(n, p)
+
+    a = pdf(d, 0:n)
+    for t=0:n
+        @test_approx_eq pdf(d, t) a[1+t]
+    end
+
+    li = rand(0:n, 2)
+    rng = minimum(li):maximum(li)
+    b = pdf(d, rng)
+    for t in rng
+        @test_approx_eq pdf(d, t) b[t - first(rng) + 1]
+    end
+
+end

test/runtests.jl

@@ -7,6 +7,7 @@ tests = [
     "continuous",
     "fit",
     "multinomial",
+    "binomial",
     "poissonbinomial",
     "dirichlet",
     "mvnormal",