Add cdf and logcdf for Multinomial distribution

Project.toml

@@ -6,6 +6,7 @@ version = "0.25.125"
 [deps]
 AliasTables = "66dad0bd-aa9a-41b7-9441-69ab47430ed8"
 ChainRulesCore = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4"
+Combinatorics = "861a8166-3701-5b0c-9a16-15d98fcdc6aa"
 DensityInterface = "b429d917-457f-4dbc-8f4c-0cc954292b1d"
 FillArrays = "1a297f60-69ca-5386-bcde-b61e274b549b"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
@@ -36,6 +37,7 @@ Aqua = "0.8.9"
 Calculus = "0.5"
 ChainRulesCore = "1"
 ChainRulesTestUtils = "1"
+Combinatorics = "1.1.0"
 DensityInterface = "0.4"
 Distributed = "<0.0.1, 1"
 FillArrays = "0.9, 0.10, 0.11, 0.12, 0.13, 1"

src/multivariate/multinomial.jl

@@ -1,3 +1,5 @@
+using Combinatorics
+
 """
 The [Multinomial distribution](http://en.wikipedia.org/wiki/Multinomial_distribution)
 generalizes the *binomial distribution*. Consider n independent draws from a Categorical
@@ -19,6 +21,7 @@ Multinomial(n, k)   # Multinomial distribution for n trials with equal probabili
                     # over 1:k
 ```
 """
+
 struct Multinomial{T<:Real, TV<:AbstractVector{T}} <: DiscreteMultivariateDistribution
     n::Int
     p::TV
@@ -162,6 +165,53 @@ function _logpdf(d::Multinomial, x::AbstractVector{T}) where T<:Real
     return s
 end
 
+"""
+    cdf(d::Multinomial, x::AbstractVector{<:Real})
+
+Compute the cumulative distribution function (CDF) of the Multinomial distribution `d` evaluated at `x`. 
+
+This calculates the probability of observing counts less than or equal to the corresponding 
+values in `x` for all categories simultaneously.
+"""
+function cdf(d::Multinomial, x::AbstractVector{<:Real})
+    n = d.n
+    k = length(d.p)
+
+    any(x .< 0) && return 0.0
+
+    sum(x) < n && return 0.0
+
+    pool = Int[]
+    for i in 1:k
+        max_capacity = min(floor(Int, x[i]), n)
+        append!(pool, fill(i, max_capacity))
+    end
+
+    total_prob = 0.0
+
+    y = zeros(Int, k)
+
+    for combo in multiset_combinations(pool, n)
+        fill!(y, 0)
+        for category in combo
+            y[category] += 1
+        end
+
+        total_prob += pdf(d, y)
+    end
+
+    return total_prob
+end
+
+"""
+    logcdf(d::Multinomial, x::AbstractVector{<:Real})
+
+Compute the log of the cumulative distribution function of the Multinomial distribution `d` evaluated at `x`.
+"""
+function logcdf(d::Multinomial, x::AbstractVector{<:Real})
+    return log(cdf(d, x))
+end
+
 # Sampling
 
 # if only a single sample is requested, no alias table is created