Merge pull request #1198 from bdeonovic/master

Faster, more accurate Poisson Binomial PMF

src/univariate/discrete/poissonbinomial.jl

@@ -29,7 +29,7 @@ struct PoissonBinomial{T<:Real} <: DiscreteUnivariateDistribution
     pmf::Vector{T}
 
     function PoissonBinomial{T}(p::AbstractArray) where {T <: Real}
-        pb = poissonbinomial_pdf_fft(p)
+        pb = poissonbinomial_pdf(p)
         @assert isprobvec(pb)
         new{T}(p, pb)
     end
@@ -113,6 +113,28 @@ end
 pdf(d::PoissonBinomial, k::Real) = insupport(d, k) ? d.pmf[k+1] : zero(eltype(d.pmf))
 logpdf(d::PoissonBinomial, k::Real) = log(pdf(d, k))
 
+# Computes the pdf of a poisson-binomial random variable using
+# simple, fast recursive formula
+#
+#      Marlin A. Thomas & Audrey E. Taub (1982) 
+#      Calculating binomial probabilities when the trial probabilities are unequal, 
+#      Journal of Statistical Computation and Simulation, 14:2, 125-131, DOI: 10.1080/00949658208810534 
+#
+function poissonbinomial_pdf(p::AbstractArray{T,1}) where {T <: Real}
+  n = length(p)
+  S = zeros(T, n+1)
+  S[1] = 1-p[1]
+  S[2] = p[1]
+  @inbounds for col in 2:n
+    for r in 1:col
+        row = col - r + 1 
+        S[row+1] = (1-p[col])*S[row+1] + p[col] * S[row]
+    end
+    S[1] *= 1-p[col]
+  end
+  return S
+end
+
 # Computes the pdf of a poisson-binomial random variable using
 # fast fourier transform
 #

test/poissonbinomial.jl

@@ -1,6 +1,37 @@
 using Distributions
 using Test
 
+function naive_esf(x::AbstractVector{T}) where T <: Real
+    n = length(x)
+    S = zeros(T, n+1)
+    states = hcat(reverse.(digits.(0:2^n-1,base=2,pad=n))...)'
+
+    r_states = vec(mapslices(sum, states, dims=2))
+
+    for r in 0:n
+        idx = findall(r_states .== r)
+        S[r+1] = sum(mapslices(x->prod(x[x .!= 0]), states[idx, :] .* x', dims=2))
+    end
+    return S
+end
+
+function naive_pb(p::AbstractVector{T}) where T <: Real
+    x = p ./ (1 .- p)
+    naive_esf(x) * prod(1 ./ (1 .+ x))
+end
+
+# test PoissonBinomial PMF algorithms
+x = [3.5118, .6219, .2905, .8450, 1.8648]
+n = length(x)
+p = x ./ (1 .+ x)
+
+naive_sol = naive_pb(p)
+
+@test Distributions.poissonbinomial_pdf_fft(p) ≈ naive_sol
+@test Distributions.poissonbinomial_pdf(p) ≈ naive_sol
+
+@test Distributions.poissonbinomial_pdf_fft(p) ≈ Distributions.poissonbinomial_pdf(p)
+
 # Test the special base where PoissonBinomial distribution reduces
 # to Binomial distribution
 for (p, n) in [(0.8, 6), (0.5, 10), (0.04, 20)]