Merge branch 'master' into rician

Project.toml

@@ -1,9 +1,10 @@
 name = "Distributions"
 uuid = "31c24e10-a181-5473-b8eb-7969acd0382f"
 authors = ["JuliaStats"]
-version = "0.25.14"
+version = "0.25.16"
 
 [deps]
+ChainRulesCore = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4"
 FillArrays = "1a297f60-69ca-5386-bcde-b61e274b549b"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 PDMats = "90014a1f-27ba-587c-ab20-58faa44d9150"
@@ -17,6 +18,7 @@ StatsBase = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91"
 StatsFuns = "4c63d2b9-4356-54db-8cca-17b64c39e42c"
 
 [compat]
+ChainRulesCore = "1"
 FillArrays = "0.9, 0.10, 0.11, 0.12"
 PDMats = "0.10, 0.11"
 QuadGK = "2"
@@ -27,6 +29,7 @@ julia = "1"
 
 [extras]
 Calculus = "49dc2e85-a5d0-5ad3-a950-438e2897f1b9"
+ChainRulesTestUtils = "cdddcdb0-9152-4a09-a978-84456f9df70a"
 Distributed = "8ba89e20-285c-5b6f-9357-94700520ee1b"
 FiniteDifferences = "26cc04aa-876d-5657-8c51-4c34ba976000"
 ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
@@ -36,4 +39,4 @@ StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [targets]
-test = ["StableRNGs", "Calculus", "Distributed", "FiniteDifferences", "ForwardDiff", "JSON", "StaticArrays", "Test"]
+test = ["StableRNGs", "Calculus", "ChainRulesTestUtils", "Distributed", "FiniteDifferences", "ForwardDiff", "JSON", "StaticArrays", "Test"]

src/Distributions.jl

@@ -25,6 +25,8 @@ import PDMats: dim, PDMat, invquad
 
 using SpecialFunctions
 
+import ChainRulesCore
+
 export
     # re-export Statistics
     mean, median, quantile, std, var, cov, cor,

src/mixtures/mixturemodel.jl

@@ -441,6 +441,20 @@ componentwise_logpdf(d::MultivariateMixture, x::AbstractVector) = componentwise_
 componentwise_logpdf(d::MultivariateMixture, x::AbstractMatrix) = componentwise_logpdf!(Matrix{eltype(x)}(undef, size(x,2), ncomponents(d)), d, x)
 
 
+function quantile(d::UnivariateMixture{Continuous}, p::Real)
+    ps = probs(d)
+    min_q, max_q = extrema(quantile(component(d, i), p) for (i, pi) in enumerate(ps) if pi > 0)
+    quantile_bisect(d, p, min_q, max_q)
+end
+
+# we also implement `median` since `median` is implemented more efficiently than
+# `quantile(d, 1//2)` for some distributions
+function median(d::UnivariateMixture{Continuous})
+    ps = probs(d)
+    min_q, max_q = extrema(median(component(d, i)) for (i, pi) in enumerate(ps) if pi > 0)
+    quantile_bisect(d, 1//2, min_q, max_q)
+end
+
 ## Sampling
 
 struct MixtureSampler{VF,VS,Sampler} <: Sampleable{VF,VS}

src/quantilealgs.jl

@@ -1,8 +1,12 @@
 # Various algorithms for computing quantile
 
-function quantile_bisect(d::ContinuousUnivariateDistribution, p::Real,
-                         lx::Real, rx::Real, tol::Real)
-
+function quantile_bisect(
+    d::ContinuousUnivariateDistribution, p::Real, lx::Real, rx::Real,
+    # base tolerance on types to support e.g. `Float32` (avoids an infinite loop)
+    # ≈ 3.7e-11 for Float64
+    # ≈ 2.4e-5 for Float32
+    tol::Real=(eps(Base.promote_typeof(float(lx), float(rx))))^(2 / 3),
+)
     # find quantile using bisect algorithm
     cl = cdf(d, lx)
     cr = cdf(d, rx)
@@ -22,7 +26,7 @@ function quantile_bisect(d::ContinuousUnivariateDistribution, p::Real,
 end
 
 quantile_bisect(d::ContinuousUnivariateDistribution, p::Real) =
-    quantile_bisect(d, p, minimum(d), maximum(d), 1.0e-12)
+    quantile_bisect(d, p, minimum(d), maximum(d))
 
 # if starting at mode, Newton is convergent for any unimodal continuous distribution, see:
 #   Göknur Giner, Gordon K. Smyth (2014)

src/univariate/continuous/normal.jl

@@ -198,7 +198,7 @@ end
 # quantile
 function quantile(d::Normal, p::Real)
     # Promote to ensure that we don't compute erfcinv in low precision and then promote
-    _p, _μ, _σ = promote(float(p), d.μ, d.σ)
+    _p, _μ, _σ = map(float, promote(p, d.μ, d.σ))
     q = xval(d, -erfcinv(2*_p) * sqrt2)
     if isnan(_p)
         return oftype(q, _p)
@@ -218,7 +218,7 @@ end
 # cquantile
 function cquantile(d::Normal, p::Real)
     # Promote to ensure that we don't compute erfcinv in low precision and then promote
-    _p, _μ, _σ = promote(float(p), d.μ, d.σ)
+    _p, _μ, _σ = map(float, promote(p, d.μ, d.σ))
     q = xval(d, erfcinv(2*_p) * sqrt2)
     if isnan(_p)
         return oftype(q, _p)

src/univariate/discrete/discretenonparametric.jl

@@ -178,71 +178,36 @@ insupport(d::DiscreteNonParametric, x::Real) =
 
 mean(d::DiscreteNonParametric) = dot(probs(d), support(d))
 
-function var(d::DiscreteNonParametric{T}) where T
-    m = mean(d)
+function var(d::DiscreteNonParametric)
     x = support(d)
     p = probs(d)
-    k = length(x)
-    σ² = zero(T)
-    for i in 1:k
-        @inbounds σ² += abs2(x[i] - m) * p[i]
-    end
-    σ²
+    return var(x, Weights(p, one(eltype(p))); corrected=false)
 end
 
-function skewness(d::DiscreteNonParametric{T}) where T
-    m = mean(d)
+function skewness(d::DiscreteNonParametric)
     x = support(d)
     p = probs(d)
-    k = length(x)
-    μ₃ = zero(T)
-    σ² = zero(T)
-    @inbounds for i in 1:k
-        d = x[i] - m
-        d²w = abs2(d) * p[i]
-        μ₃ += d * d²w
-        σ² += d²w
-    end
-    μ₃ / (σ² * sqrt(σ²))
+    return skewness(x, Weights(p, one(eltype(p))))
 end
 
-function kurtosis(d::DiscreteNonParametric{T}) where T
-    m = mean(d)
+function kurtosis(d::DiscreteNonParametric)
     x = support(d)
     p = probs(d)
-    k = length(x)
-    μ₄ = zero(T)
-    σ² = zero(T)
-    @inbounds for i in 1:k
-        d² = abs2(x[i] - m)
-        d²w = d² * p[i]
-        μ₄ += d² * d²w
-        σ² += d²w
-    end
-    μ₄ / abs2(σ²) - 3
+    return kurtosis(x, Weights(p, one(eltype(p))))
 end
 
 entropy(d::DiscreteNonParametric) = entropy(probs(d))
 entropy(d::DiscreteNonParametric, b::Real) = entropy(probs(d), b)
 
-mode(d::DiscreteNonParametric) = support(d)[argmax(probs(d))]
-function modes(d::DiscreteNonParametric{T,P}) where {T,P}
+function mode(d::DiscreteNonParametric)
     x = support(d)
     p = probs(d)
-    k = length(x)
-    mds = T[]
-    max_p = zero(P)
-    @inbounds for i in 1:k
-        pi = p[i]
-        xi = x[i]
-        if pi > max_p
-            max_p = pi
-            mds = [xi]
-        elseif pi == max_p
-            push!(mds, xi)
-        end
-    end
-    mds
+    return mode(x, Weights(p, one(eltype(p))))
+end
+function modes(d::DiscreteNonParametric)
+    x = support(d)
+    p = probs(d)
+    return modes(x, Weights(p, one(eltype(p))))
 end
 
 function mgf(d::DiscreteNonParametric, t::Real)

src/univariate/discrete/poissonbinomial.jl

@@ -201,3 +201,70 @@ end
 #### Sampling
 
 sampler(d::PoissonBinomial) = PoissBinAliasSampler(d)
+
+## ChainRules definitions
+
+# Compute matrix of partial derivatives [∂P(X=j-1)/∂pᵢ]_{i=1,…,n; j=1,…,n+1}
+#
+# This implementation uses the same dynamic programming "trick" as for the computation of
+# the primals.
+#
+# Reference (for the primal):
+#
+#      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_partialderivatives(p::AbstractVector{<:Real})
+    n = length(p)
+    A = zeros(eltype(p), n, n + 1)
+    @inbounds for j in 1:n
+        A[j, end] = 1
+    end
+    @inbounds for (i, pi) in enumerate(p)
+        qi = 1 - pi
+        for k in (n - i + 1):n
+            kp1 = k + 1
+            for j in 1:(i - 1)
+                A[j, k] = pi * A[j, k] + qi * A[j, kp1]
+            end
+            for j in (i+1):n
+                A[j, k] = pi * A[j, k] + qi * A[j, kp1]
+            end
+        end
+        for j in 1:(i-1)
+            A[j, end] *= pi
+        end
+        for j in (i+1):n
+            A[j, end] *= pi
+        end
+    end
+    @inbounds for j in 1:n, i in 1:n
+        A[i, j] -= A[i, j+1]
+    end
+    return A
+end
+
+for f in (:poissonbinomial_pdf, :poissonbinomial_pdf_fft)
+    pullback = Symbol(f, :_pullback)
+    @eval begin
+        function ChainRulesCore.frule(
+            (_, Δp)::Tuple{<:Any,<:AbstractVector{<:Real}}, ::typeof($f), p::AbstractVector{<:Real}
+        )
+            y = $f(p)
+            A = poissonbinomial_pdf_partialderivatives(p)
+            return y, A' * Δp
+        end
+        function ChainRulesCore.rrule(::typeof($f), p::AbstractVector{<:Real})
+            y = $f(p)
+            A = poissonbinomial_pdf_partialderivatives(p)
+            function $pullback(Δy)
+                p̄ = ChainRulesCore.InplaceableThunk(
+                    Δ -> LinearAlgebra.mul!(Δ, A, Δy, true, true),
+                    ChainRulesCore.@thunk(A * Δy),
+                )
+                return ChainRulesCore.NoTangent(), p̄
+            end
+            return y, $pullback
+        end
+    end
+end

src/univariates.jl

@@ -216,7 +216,7 @@ std(d::UnivariateDistribution) = sqrt(var(d))
 
 Return the median value of distribution `d`.
 """
-median(d::UnivariateDistribution) = quantile(d, 0.5)
+median(d::UnivariateDistribution) = quantile(d, 1//2)
 
 """
     modes(d::UnivariateDistribution)

test/discretenonparametric.jl

@@ -163,4 +163,11 @@ d = DiscreteNonParametric([1, 2], [0, 1])
 
 d = DiscreteNonParametric([2, 1], [1, 0])
 @test iszero(count(isone(rand(d)) for _ in 1:100))
+
+    @testset "type stability" begin
+        d = DiscreteNonParametric([0, 1], [0.5, 0.5])
+        @inferred(var(d))
+        @inferred(kurtosis(d))
+        @inferred(skewness(d))
+    end
 end

test/mixture.jl

@@ -69,8 +69,17 @@ function test_mixture(g::UnivariateMixture, n::Int, ns::Int,
     @test @inferred(componentwise_pdf(g, X)) ≈ P0
     @test @inferred(componentwise_logpdf(g, X)) ≈ LP0
 
+    # quantile
+    αs = float(partype(g))[0.0; 0.49; 0.5; 0.51; 1.0]
+    for α in αs
+        @test cdf(g, @inferred(quantile(g, α))) ≈ α
+    end
+    @test @inferred(median(g)) ≈ quantile(g, 1//2)
+
     # sampling
-    if (T <: AbstractFloat)
+    # sampling does not work with `Float32` since `AliasTable` does not support `Float32`
+    # Ref: https://github.com/JuliaStats/StatsBase.jl/issues/158
+    if T <: AbstractFloat && eltype(probs(g)) === Float64
         if ismissing(rng)
             Xs = rand(g, ns)
         else
@@ -172,6 +181,17 @@ end
          "rand(rng, ...)" => MersenneTwister(123))
 
     @testset "Testing UnivariateMixture" begin
+        g_u = MixtureModel([Normal(), Normal()])
+        @test isa(g_u, MixtureModel{Univariate, Continuous, <:Normal})
+        @test ncomponents(g_u) == 2
+        test_mixture(g_u, 1000, 10^6, rng)
+        test_params(g_u)
+        @test minimum(g_u) == -Inf
+        @test maximum(g_u) == Inf
+        @test extrema(g_u) == (-Inf, Inf)
+        @test @inferred(median(g_u)) === 0.0
+        @test @inferred(quantile(g_u, 0.5f0)) === 0.0
+
         g_u = MixtureModel(Normal{Float64}, [(0.0, 1.0), (2.0, 1.0), (-4.0, 1.5)], [0.2, 0.5, 0.3])
         @test isa(g_u, MixtureModel{Univariate,Continuous,<:Normal})
         @test ncomponents(g_u) == 3
@@ -181,6 +201,17 @@ end
         @test maximum(g_u) == Inf
         @test extrema(g_u) == (-Inf, Inf)
 
+        g_u = MixtureModel(Normal{Float32}, [(0f0, 1f0), (0f0, 2f0)], [0.4f0, 0.6f0])
+        @test isa(g_u, MixtureModel{Univariate,Continuous,<:Normal})
+        @test ncomponents(g_u) == 2
+        test_mixture(g_u, 1000, 10^6, rng)
+        test_params(g_u)
+        @test minimum(g_u) == -Inf
+        @test maximum(g_u) == Inf
+        @test extrema(g_u) == (-Inf, Inf)
+        @test @inferred(median(g_u)) === 0f0
+        @test @inferred(quantile(g_u, 0.5f0)) === 0f0
+
         g_u = MixtureModel([TriangularDist(-1,2,0),TriangularDist(-.5,3,1),TriangularDist(-2,0,-1)])
         @test minimum(g_u) ≈ -2.0
         @test maximum(g_u) ≈ 3.0

test/normal.jl

@@ -150,6 +150,8 @@ end
     @test @inferred(quantile(Normal(1.0f0, 0.0f0), 0.5f0)) ===  1.0f0
     @test isnan_type(Float32, @inferred(quantile(Normal(1.0f0, 0.0f0), NaN32)))
     @test @inferred(quantile(Normal(1//1, 0//1), 1//2))    ===  1.0
+    @test @inferred(quantile(Normal(1f0, 0f0), 1//2))      ===  1f0
+    @test @inferred(quantile(Normal(1f0, 0.0), 1//2))      ===  1.0
 
     @test @inferred(cquantile(Normal(1.0, 0.0), 0.0f0))     ===  Inf
     @test @inferred(cquantile(Normal(1.0, 0.0f0), 1.0))     === -Inf
@@ -160,6 +162,8 @@ end
     @test @inferred(cquantile(Normal(1.0f0, 0.0f0), 0.5f0)) ===  1.0f0
     @test isnan_type(Float32, @inferred(cquantile(Normal(1.0f0, 0.0f0), NaN32)))
     @test @inferred(cquantile(Normal(1//1, 0//1), 1//2))    ===  1.0
+    @test @inferred(cquantile(Normal(1f0, 0f0), 1//2))      ===  1f0
+    @test @inferred(cquantile(Normal(1f0, 0.0), 1//2))      ===  1.0
 end
 
 @testset "Normal: Sampling with integer-valued parameters" begin