From 87aebc29b2b9608801b70aae09fbc1d2dad56e3f Mon Sep 17 00:00:00 2001 From: "dependabot[bot]" <49699333+dependabot[bot]@users.noreply.github.com> Date: Mon, 25 Sep 2023 09:53:40 +0200 Subject: [PATCH 01/26] Bump actions/checkout from 3 to 4 (#1764) Bumps [actions/checkout](https://github.com/actions/checkout) from 3 to 4. - [Release notes](https://github.com/actions/checkout/releases) - [Changelog](https://github.com/actions/checkout/blob/main/CHANGELOG.md) - [Commits](https://github.com/actions/checkout/compare/v3...v4) --- updated-dependencies: - dependency-name: actions/checkout dependency-type: direct:production update-type: version-update:semver-major ... Signed-off-by: dependabot[bot] Co-authored-by: dependabot[bot] <49699333+dependabot[bot]@users.noreply.github.com> --- .github/workflows/CI.yml | 4 ++-- .github/workflows/DocPreviewCleanup.yml | 2 +- .github/workflows/IntegrationTest.yml | 4 ++-- 3 files changed, 5 insertions(+), 5 deletions(-) diff --git a/.github/workflows/CI.yml b/.github/workflows/CI.yml index b7ac9ffd17..ec2f9a00c9 100644 --- a/.github/workflows/CI.yml +++ b/.github/workflows/CI.yml @@ -34,7 +34,7 @@ jobs: arch: - x64 steps: - - uses: actions/checkout@v3 + - uses: actions/checkout@v4 - uses: julia-actions/setup-julia@v1 with: version: ${{ matrix.version }} @@ -64,7 +64,7 @@ jobs: name: Documentation runs-on: ubuntu-latest steps: - - uses: actions/checkout@v3 + - uses: actions/checkout@v4 - uses: julia-actions/setup-julia@v1 with: version: '1' diff --git a/.github/workflows/DocPreviewCleanup.yml b/.github/workflows/DocPreviewCleanup.yml index 844c4c9a51..94b9f70ade 100644 --- a/.github/workflows/DocPreviewCleanup.yml +++ b/.github/workflows/DocPreviewCleanup.yml @@ -9,7 +9,7 @@ jobs: runs-on: ubuntu-latest steps: - name: Checkout gh-pages branch - uses: actions/checkout@v3 + uses: actions/checkout@v4 with: ref: gh-pages - name: Delete preview and history + push changes diff --git a/.github/workflows/IntegrationTest.yml b/.github/workflows/IntegrationTest.yml index 29af889646..3aec5d9002 100644 --- a/.github/workflows/IntegrationTest.yml +++ b/.github/workflows/IntegrationTest.yml @@ -31,14 +31,14 @@ jobs: #- {user: TuringLang, repo: DistributionsAD.jl, group: ForwardDiff} takes > 1 hour steps: - - uses: actions/checkout@v3 + - uses: actions/checkout@v4 - uses: julia-actions/setup-julia@v1 with: version: 1 arch: x64 - uses: julia-actions/julia-buildpkg@latest - name: Clone Downstream - uses: actions/checkout@v3 + uses: actions/checkout@v4 with: repository: ${{ matrix.package.user }}/${{ matrix.package.repo }} path: downstream From 902c98e1438b140a99266a4129ce0575d0366ef9 Mon Sep 17 00:00:00 2001 From: Chris Elrod Date: Tue, 26 Sep 2023 18:13:13 -0400 Subject: [PATCH 02/26] Fix undef ref in lkj_chol sampling (#1782) * Fix undef ref in lkj_chol sampling * Add test for dim=1 lkjcholesky * Update test/cholesky/lkjcholesky.jl Co-authored-by: David Widmann * Update test/cholesky/lkjcholesky.jl Co-authored-by: Seth Axen * Update test/cholesky/lkjcholesky.jl Co-authored-by: David Widmann --------- Co-authored-by: David Widmann Co-authored-by: Seth Axen --- src/cholesky/lkjcholesky.jl | 2 +- test/cholesky/lkjcholesky.jl | 1 + 2 files changed, 2 insertions(+), 1 deletion(-) diff --git a/src/cholesky/lkjcholesky.jl b/src/cholesky/lkjcholesky.jl index 3540e31e29..7556620f60 100644 --- a/src/cholesky/lkjcholesky.jl +++ b/src/cholesky/lkjcholesky.jl @@ -238,7 +238,7 @@ function _lkj_cholesky_onion_tri!( # equivalent steps in algorithm in reference are marked. @assert size(A) == (d, d) A[1, 1] = 1 - d > 1 || return R + d > 1 || return A β = η + (d - 2)//2 # 1. Initialization w0 = 2 * rand(rng, Beta(β, β)) - 1 diff --git a/test/cholesky/lkjcholesky.jl b/test/cholesky/lkjcholesky.jl index a0e8283436..b9afb59beb 100644 --- a/test/cholesky/lkjcholesky.jl +++ b/test/cholesky/lkjcholesky.jl @@ -181,6 +181,7 @@ using FiniteDifferences nkstests = 4 # use for appropriate Bonferroni correction for KS test @testset "rand" begin + @test rand(LKJCholesky(1, 0.5)).factors == ones(1, 1) @testset for p in (2, 4, 10), η in (0.5, 1, 3), uplo in ('L', 'U') d = LKJCholesky(p, η, uplo) test_draw(d, rand(rng, d)) From b21e515dd4684d2d0f2c706c573ed6eb98ce537f Mon Sep 17 00:00:00 2001 From: David Widmann Date: Wed, 27 Sep 2023 00:13:33 +0200 Subject: [PATCH 03/26] Update Project.toml --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index 278fd45246..bd2c0fe893 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "Distributions" uuid = "31c24e10-a181-5473-b8eb-7969acd0382f" authors = ["JuliaStats"] -version = "0.25.100" +version = "0.25.101" [deps] ChainRulesCore = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4" From cd5d4cc7ca98c165bfff07b23a14afd66ff2e9db Mon Sep 17 00:00:00 2001 From: Seth Axen Date: Thu, 28 Sep 2023 13:43:09 +0200 Subject: [PATCH 04/26] Add MvLogitNormal (#1774) * Create MvLogitNormal * Add MvLogitNormal to docs * Simplify constructors * Fix conversions * Rearrange code * Fix computation of -Inf * Add meanform and canonform * Add back type constructor * Add MvLogitNormal tests * Update and test show method * Fix testset name * Fix for older Julia versions * Restrict testing of `show` method to newer versions * Add kldivergence tests * Improve documentation * Remove constructor with type and AbstractMvNormal params * Update show method * Update docstring * Remove reference to Dirichlet * Apply suggestions from code review Co-authored-by: David Widmann --------- Co-authored-by: David Widmann --- docs/src/multivariate.md | 1 + src/Distributions.jl | 1 + src/multivariate/mvlogitnormal.jl | 140 +++++++++++++++++++++++++ src/multivariates.jl | 1 + test/multivariate/mvlogitnormal.jl | 158 +++++++++++++++++++++++++++++ test/runtests.jl | 1 + 6 files changed, 302 insertions(+) create mode 100644 src/multivariate/mvlogitnormal.jl create mode 100644 test/multivariate/mvlogitnormal.jl diff --git a/docs/src/multivariate.md b/docs/src/multivariate.md index 2ae5ef59bd..c4e7c17646 100644 --- a/docs/src/multivariate.md +++ b/docs/src/multivariate.md @@ -55,6 +55,7 @@ Multinomial Distributions.AbstractMvNormal MvNormal MvNormalCanon +MvLogitNormal MvLogNormal Dirichlet Product diff --git a/src/Distributions.jl b/src/Distributions.jl index 1e2d580c55..7e4f6bd62b 100644 --- a/src/Distributions.jl +++ b/src/Distributions.jl @@ -122,6 +122,7 @@ export Logistic, LogNormal, LogUniform, + MvLogitNormal, LogitNormal, MatrixBeta, MatrixFDist, diff --git a/src/multivariate/mvlogitnormal.jl b/src/multivariate/mvlogitnormal.jl new file mode 100644 index 0000000000..0d60ddf654 --- /dev/null +++ b/src/multivariate/mvlogitnormal.jl @@ -0,0 +1,140 @@ +""" + MvLogitNormal{<:AbstractMvNormal} + +The [multivariate logit-normal distribution](https://en.wikipedia.org/wiki/Logit-normal_distribution#Multivariate_generalization) +is a multivariate generalization of [`LogitNormal`](@ref) capable of handling correlations +between variables. + +If ``\\mathbf{y} \\sim \\mathrm{MvNormal}(\\boldsymbol{\\mu}, \\boldsymbol{\\Sigma})`` is a +length ``d-1`` vector, then +```math +\\mathbf{x} = \\operatorname{softmax}\\left(\\begin{bmatrix}\\mathbf{y} \\\\ 0 \\end{bmatrix}\\right) \\sim \\mathrm{MvLogitNormal}(\\boldsymbol{\\mu}, \\boldsymbol{\\Sigma}) +``` +is a length ``d`` probability vector. + +```julia +MvLogitNormal(μ, Σ) # MvLogitNormal with y ~ MvNormal(μ, Σ) +MvLogitNormal(MvNormal(μ, Σ)) # same as above +MvLogitNormal(MvNormalCanon(μ, J)) # MvLogitNormal with y ~ MvNormalCanon(μ, J) +``` + +# Fields + +- `normal::AbstractMvNormal`: contains the ``d-1``-dimensional distribution of ``y`` +""" +struct MvLogitNormal{D<:AbstractMvNormal} <: ContinuousMultivariateDistribution + normal::D + MvLogitNormal{D}(normal::D) where {D<:AbstractMvNormal} = new{D}(normal) +end +MvLogitNormal(d::AbstractMvNormal) = MvLogitNormal{typeof(d)}(d) +MvLogitNormal(args...) = MvLogitNormal(MvNormal(args...)) + +function Base.show(io::IO, d::MvLogitNormal; indent::String=" ") + print(io, distrname(d)) + println(io, "(") + normstr = strip(sprint(show, d.normal; context=IOContext(io))) + normstr = replace(normstr, "\n" => "\n$indent") + print(io, indent) + println(io, normstr) + println(io, ")") +end + +# Conversions + +function convert(::Type{MvLogitNormal{D}}, d::MvLogitNormal) where {D} + return MvLogitNormal(convert(D, d.normal)) +end +Base.convert(::Type{MvLogitNormal{D}}, d::MvLogitNormal{D}) where {D} = d + +meanform(d::MvLogitNormal{<:MvNormalCanon}) = MvLogitNormal(meanform(d.normal)) +canonform(d::MvLogitNormal{<:MvNormal}) = MvLogitNormal(canonform(d.normal)) + +# Properties + +length(d::MvLogitNormal) = length(d.normal) + 1 +Base.eltype(::Type{<:MvLogitNormal{D}}) where {D} = eltype(D) +Base.eltype(d::MvLogitNormal) = eltype(d.normal) +params(d::MvLogitNormal) = params(d.normal) +@inline partype(d::MvLogitNormal) = partype(d.normal) + +location(d::MvLogitNormal) = mean(d.normal) +minimum(d::MvLogitNormal) = fill(zero(eltype(d)), length(d)) +maximum(d::MvLogitNormal) = fill(oneunit(eltype(d)), length(d)) + +function insupport(d::MvLogitNormal, x::AbstractVector{<:Real}) + return length(d) == length(x) && all(≥(0), x) && sum(x) ≈ 1 +end + +# Evaluation + +function _logpdf(d::MvLogitNormal, x::AbstractVector{<:Real}) + if !insupport(d, x) + return oftype(logpdf(d.normal, _inv_softmax1(abs.(x))), -Inf) + else + return logpdf(d.normal, _inv_softmax1(x)) - sum(log, x) + end +end + +function gradlogpdf(d::MvLogitNormal, x::AbstractVector{<:Real}) + y = _inv_softmax1(x) + ∂y = gradlogpdf(d.normal, y) + ∂x = (vcat(∂y, -sum(∂y)) .- 1) ./ x + return ∂x +end + +# Statistics + +kldivergence(p::MvLogitNormal, q::MvLogitNormal) = kldivergence(p.normal, q.normal) + +# Sampling + +function _rand!(rng::AbstractRNG, d::MvLogitNormal, x::AbstractVecOrMat{<:Real}) + y = @views _drop1(x) + rand!(rng, d.normal, y) + _softmax1!(x, y) + return x +end + +# Fitting + +function fit_mle(::Type{MvLogitNormal{D}}, x::AbstractMatrix{<:Real}; kwargs...) where {D} + y = similar(x, size(x, 1) - 1, size(x, 2)) + map(_inv_softmax1!, eachcol(y), eachcol(x)) + normal = fit_mle(D, y; kwargs...) + return MvLogitNormal(normal) +end +function fit_mle(::Type{MvLogitNormal}, x::AbstractMatrix{<:Real}; kwargs...) + return fit_mle(MvLogitNormal{MvNormal}, x; kwargs...) +end + +# Utility + +function _softmax1!(x::AbstractVector, y::AbstractVector) + u = max(0, maximum(y)) + _drop1(x) .= exp.(y .- u) + x[end] = exp(-u) + LinearAlgebra.normalize!(x, 1) + return x +end +function _softmax1!(x::AbstractMatrix, y::AbstractMatrix) + map(_softmax1!, eachcol(x), eachcol(y)) + return x +end + +_drop1(x::AbstractVector) = @views x[firstindex(x, 1):(end - 1)] +_drop1(x::AbstractMatrix) = @views x[firstindex(x, 1):(end - 1), :] + +_last1(x::AbstractVector) = x[end] +_last1(x::AbstractMatrix) = @views x[end, :] + +function _inv_softmax1!(y::AbstractVecOrMat, x::AbstractVecOrMat) + x₋ = _drop1(x) + xd = _last1(x) + @. y = log(x₋) - log(xd) + return y +end +function _inv_softmax1(x::AbstractVecOrMat) + y = similar(_drop1(x)) + _inv_softmax1!(y, x) + return y +end diff --git a/src/multivariates.jl b/src/multivariates.jl index 7a6f926a65..56d91233cf 100644 --- a/src/multivariates.jl +++ b/src/multivariates.jl @@ -115,6 +115,7 @@ for fname in ["dirichlet.jl", "jointorderstatistics.jl", "mvnormal.jl", "mvnormalcanon.jl", + "mvlogitnormal.jl", "mvlognormal.jl", "mvtdist.jl", "product.jl", # deprecated diff --git a/test/multivariate/mvlogitnormal.jl b/test/multivariate/mvlogitnormal.jl new file mode 100644 index 0000000000..cb53c9b377 --- /dev/null +++ b/test/multivariate/mvlogitnormal.jl @@ -0,0 +1,158 @@ +# Tests on Multivariate Logit-Normal distributions +using Distributions +using ForwardDiff +using LinearAlgebra +using Random +using Test + +####### Core testing procedure + +function test_mvlogitnormal(d::MvLogitNormal; nsamples::Int=10^6) + @test d.normal isa AbstractMvNormal + dnorm = d.normal + + @testset "properties" begin + @test length(d) == length(dnorm) + 1 + @test params(d) == params(dnorm) + @test partype(d) == partype(dnorm) + @test eltype(d) == eltype(dnorm) + @test eltype(typeof(d)) == eltype(typeof(dnorm)) + @test location(d) == mean(dnorm) + @test minimum(d) == fill(0, length(d)) + @test maximum(d) == fill(1, length(d)) + @test insupport(d, normalize(rand(length(d)), 1)) + @test !insupport(d, normalize(rand(length(d) + 1), 1)) + @test !insupport(d, rand(length(d))) + x = rand(length(d) - 1) + x = vcat(x, -sum(x)) + @test !insupport(d, x) + end + + @testset "conversions" begin + @test convert(typeof(d), d) === d + T = partype(d) <: Float64 ? Float32 : Float64 + if dnorm isa MvNormal + @test convert(MvLogitNormal{MvNormal{T}}, d).normal == + convert(MvNormal{T}, dnorm) + @test partype(convert(MvLogitNormal{MvNormal{T}}, d)) <: T + @test canonform(d) isa MvLogitNormal{<:MvNormalCanon} + @test canonform(d).normal == canonform(dnorm) + elseif dnorm isa MvNormalCanon + @test convert(MvLogitNormal{MvNormalCanon{T}}, d).normal == + convert(MvNormalCanon{T}, dnorm) + @test partype(convert(MvLogitNormal{MvNormalCanon{T}}, d)) <: T + @test meanform(d) isa MvLogitNormal{<:MvNormal} + @test meanform(d).normal == meanform(dnorm) + end + end + + @testset "sampling" begin + X = rand(d, nsamples) + Y = @views log.(X[1:(end - 1), :]) .- log.(X[end, :]') + Ymean = vec(mean(Y; dims=2)) + Ycov = cov(Y; dims=2) + for i in 1:(length(d) - 1) + @test isapprox( + Ymean[i], mean(dnorm)[i], atol=sqrt(var(dnorm)[i] / nsamples) * 8 + ) + end + for i in 1:(length(d) - 1), j in 1:(length(d) - 1) + @test isapprox( + Ycov[i, j], + cov(dnorm)[i, j], + atol=sqrt(prod(var(dnorm)[[i, j]]) / nsamples) * 20, + ) + end + end + + @testset "fitting" begin + X = rand(d, nsamples) + dfit = fit_mle(MvLogitNormal, X) + dfit_norm = dfit.normal + for i in 1:(length(d) - 1) + @test isapprox( + mean(dfit_norm)[i], mean(dnorm)[i], atol=sqrt(var(dnorm)[i] / nsamples) * 8 + ) + end + for i in 1:(length(d) - 1), j in 1:(length(d) - 1) + @test isapprox( + cov(dfit_norm)[i, j], + cov(dnorm)[i, j], + atol=sqrt(prod(var(dnorm)[[i, j]]) / nsamples) * 20, + ) + end + @test fit_mle(MvLogitNormal{IsoNormal}, X) isa MvLogitNormal{<:IsoNormal} + end + + @testset "evaluation" begin + X = rand(d, nsamples) + for i in 1:min(100, nsamples) + @test @inferred(logpdf(d, X[:, i])) ≈ log(pdf(d, X[:, i])) + if dnorm isa MvNormal + @test @inferred(gradlogpdf(d, X[:, i])) ≈ + ForwardDiff.gradient(x -> logpdf(d, x), X[:, i]) + end + end + @test logpdf(d, X) ≈ log.(pdf(d, X)) + @test isequal(logpdf(d, zeros(length(d))), -Inf) + @test isequal(logpdf(d, ones(length(d))), -Inf) + @test isequal(pdf(d, zeros(length(d))), 0) + @test isequal(pdf(d, ones(length(d))), 0) + end +end + +@testset "Results MvLogitNormal consistent with univariate LogitNormal" begin + μ = randn() + σ = rand() + d = MvLogitNormal([μ], fill(σ^2, 1, 1)) + duni = LogitNormal(μ, σ) + @test location(d) ≈ [location(duni)] + x = normalize(rand(2), 1) + @test logpdf(d, x) ≈ logpdf(duni, x[1]) + @test pdf(d, x) ≈ pdf(duni, x[1]) + @test (Random.seed!(9274); rand(d)[1]) ≈ (Random.seed!(9274); rand(duni)) +end + +###### General Testing + +@testset "MvLogitNormal tests" begin + mvnorm_params = [ + (randn(5), I * rand()), + (randn(4), Diagonal(rand(4))), + (Diagonal(rand(6)),), + (randn(5), exp(Symmetric(randn(5, 5)))), + (exp(Symmetric(randn(5, 5))),), + ] + @testset "wraps MvNormal" begin + @testset "$(typeof(prms))" for prms in mvnorm_params + d = MvLogitNormal(prms...) + @test d == MvLogitNormal(MvNormal(prms...)) + test_mvlogitnormal(d; nsamples=10^4) + end + end + @testset "wraps MvNormalCanon" begin + @testset "$(typeof(prms))" for prms in mvnorm_params + d = MvLogitNormal(MvNormalCanon(prms...)) + test_mvlogitnormal(d; nsamples=10^4) + end + end + + @testset "kldivergence" begin + d1 = MvLogitNormal(randn(5), exp(Symmetric(randn(5, 5)))) + d2 = MvLogitNormal(randn(5), exp(Symmetric(randn(5, 5)))) + @test kldivergence(d1, d2) ≈ kldivergence(d1.normal, d2.normal) + end + + VERSION ≥ v"1.8" && @testset "show" begin + d = MvLogitNormal([1.0, 2.0, 3.0], Diagonal([4.0, 5.0, 6.0])) + @test sprint(show, d) === """ + MvLogitNormal{DiagNormal}( + DiagNormal( + dim: 3 + μ: [1.0, 2.0, 3.0] + Σ: [4.0 0.0 0.0; 0.0 5.0 0.0; 0.0 0.0 6.0] + ) + ) + """ + end +end diff --git a/test/runtests.jl b/test/runtests.jl index cb724278b4..ce3f16b799 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -26,6 +26,7 @@ const tests = [ "univariate/continuous/uniform", "univariate/continuous/lognormal", "multivariate/mvnormal", + "multivariate/mvlogitnormal", "multivariate/mvlognormal", "types", # extra file compared to /src "utils", From e407fa5fd098e50df51801c6d062946eac7a7d0f Mon Sep 17 00:00:00 2001 From: David Widmann Date: Thu, 28 Sep 2023 13:43:46 +0200 Subject: [PATCH 05/26] Update Project.toml --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index bd2c0fe893..92dc0b54d9 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "Distributions" uuid = "31c24e10-a181-5473-b8eb-7969acd0382f" authors = ["JuliaStats"] -version = "0.25.101" +version = "0.25.102" [deps] ChainRulesCore = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4" From e666d74abb99ea517e6c90d18c3e73b882b7e62b Mon Sep 17 00:00:00 2001 From: David Widmann Date: Fri, 3 Nov 2023 15:06:37 +0100 Subject: [PATCH 06/26] Move test utilities to an extension (#1791) * Move test utilities to an extension * Fix signature and docstring * Also qualify AbstractRNG * Fix Julia < 1.9 * Fix for 1.3? * Simplify the TestUtils stub --- Project.toml | 4 +- ext/DistributionsTestExt.jl | 101 +++++++++++++++++++++++++++++++++++ src/Distributions.jl | 7 +-- src/test_utils.jl | 103 +++++------------------------------- 4 files changed, 121 insertions(+), 94 deletions(-) create mode 100644 ext/DistributionsTestExt.jl diff --git a/Project.toml b/Project.toml index 92dc0b54d9..cf6ee221ca 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "Distributions" uuid = "31c24e10-a181-5473-b8eb-7969acd0382f" authors = ["JuliaStats"] -version = "0.25.102" +version = "0.25.103" [deps] ChainRulesCore = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4" @@ -22,10 +22,12 @@ Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" [weakdeps] ChainRulesCore = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4" DensityInterface = "b429d917-457f-4dbc-8f4c-0cc954292b1d" +Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" [extensions] DistributionsChainRulesCoreExt = "ChainRulesCore" DistributionsDensityInterfaceExt = "DensityInterface" +DistributionsTestExt = "Test" [compat] ChainRulesCore = "1" diff --git a/ext/DistributionsTestExt.jl b/ext/DistributionsTestExt.jl new file mode 100644 index 0000000000..2665d72fae --- /dev/null +++ b/ext/DistributionsTestExt.jl @@ -0,0 +1,101 @@ +module DistributionsTestExt + +using Distributions +using Distributions.LinearAlgebra +using Distributions.Random +using Test + +__rand(::Nothing, args...) = rand(args...) +__rand(rng::AbstractRNG, args...) = rand(rng, args...) + +__rand!(::Nothing, args...) = rand!(args...) +__rand!(rng::AbstractRNG, args...) = rand!(rng, args...) + +""" + test_mvnormal( + g::AbstractMvNormal, + n_tsamples::Int=10^6, + rng::Union{Random.AbstractRNG, Nothing}=nothing, + ) + +Test that `AbstractMvNormal` implements the expected API. + +!!! Note + On Julia >= 1.9, you have to load the `Test` standard library to be able to use + this function. +""" +function Distributions.TestUtils.test_mvnormal( + g::AbstractMvNormal, n_tsamples::Int=10^6, rng::Union{AbstractRNG, Nothing}=nothing +) + d = length(g) + μ = mean(g) + Σ = cov(g) + @test length(μ) == d + @test size(Σ) == (d, d) + @test var(g) ≈ diag(Σ) + @test entropy(g) ≈ 0.5 * logdet(2π * ℯ * Σ) + ldcov = logdetcov(g) + @test ldcov ≈ logdet(Σ) + vs = diag(Σ) + @test g == typeof(g)(params(g)...) + @test g == deepcopy(g) + @test minimum(g) == fill(-Inf, d) + @test maximum(g) == fill(Inf, d) + @test extrema(g) == (minimum(g), maximum(g)) + @test isless(extrema(g)...) + + # test sampling for AbstractMatrix (here, a SubArray): + subX = view(__rand(rng, d, 2d), :, 1:d) + @test isa(__rand!(rng, g, subX), SubArray) + + # sampling + @test isa(__rand(rng, g), Vector{Float64}) + X = __rand(rng, g, n_tsamples) + emp_mu = vec(mean(X, dims=2)) + Z = X .- emp_mu + emp_cov = (Z * Z') * inv(n_tsamples) + + mean_atols = 8 .* sqrt.(vs ./ n_tsamples) + cov_atols = 10 .* sqrt.(vs .* vs') ./ sqrt.(n_tsamples) + for i = 1:d + @test isapprox(emp_mu[i], μ[i], atol=mean_atols[i]) + end + for i = 1:d, j = 1:d + @test isapprox(emp_cov[i,j], Σ[i,j], atol=cov_atols[i,j]) + end + + X = rand(MersenneTwister(14), g, n_tsamples) + Y = rand(MersenneTwister(14), g, n_tsamples) + @test X == Y + emp_mu = vec(mean(X, dims=2)) + Z = X .- emp_mu + emp_cov = (Z * Z') * inv(n_tsamples) + for i = 1:d + @test isapprox(emp_mu[i] , μ[i] , atol=mean_atols[i]) + end + for i = 1:d, j = 1:d + @test isapprox(emp_cov[i,j], Σ[i,j], atol=cov_atols[i,j]) + end + + + # evaluation of sqmahal & logpdf + U = X .- μ + sqm = vec(sum(U .* (Σ \ U), dims=1)) + for i = 1:min(100, n_tsamples) + @test sqmahal(g, X[:,i]) ≈ sqm[i] + end + @test sqmahal(g, X) ≈ sqm + + lp = -0.5 .* sqm .- 0.5 * (d * log(2.0 * pi) + ldcov) + for i = 1:min(100, n_tsamples) + @test logpdf(g, X[:,i]) ≈ lp[i] + end + @test logpdf(g, X) ≈ lp + + # log likelihood + @test loglikelihood(g, X) ≈ sum(i -> Distributions._logpdf(g, X[:,i]), 1:n_tsamples) + @test loglikelihood(g, X[:, 1]) ≈ logpdf(g, X[:, 1]) + @test loglikelihood(g, [X[:, i] for i in axes(X, 2)]) ≈ loglikelihood(g, X) +end + +end # module \ No newline at end of file diff --git a/src/Distributions.jl b/src/Distributions.jl index 7e4f6bd62b..ceb6063c70 100644 --- a/src/Distributions.jl +++ b/src/Distributions.jl @@ -316,15 +316,16 @@ include("mixtures/unigmm.jl") # Interface for StatsAPI include("statsapi.jl") +# Testing utilities for other packages which implement distributions. +include("test_utils.jl") + # Extensions: Implementation of DensityInterface and ChainRulesCore API if !isdefined(Base, :get_extension) include("../ext/DistributionsChainRulesCoreExt/DistributionsChainRulesCoreExt.jl") include("../ext/DistributionsDensityInterfaceExt.jl") + include("../ext/DistributionsTestExt.jl") end -# Testing utilities for other packages which implement distributions. -include("test_utils.jl") - include("deprecates.jl") """ diff --git a/src/test_utils.jl b/src/test_utils.jl index 41e0d588ae..6f74f80ec3 100644 --- a/src/test_utils.jl +++ b/src/test_utils.jl @@ -1,96 +1,19 @@ module TestUtils -using Distributions -using LinearAlgebra -using Random -using Test - - -__rand(::Nothing, args...) = rand(args...) -__rand(rng::AbstractRNG, args...) = rand(rng, args...) - -__rand!(::Nothing, args...) = rand!(args...) -__rand!(rng::AbstractRNG, args...) = rand!(rng, args...) - -""" - test_mvnormal( - g::AbstractMvNormal, n_tsamples::Int=10^6, rng::AbstractRNG=Random.default_rng() - ) - -Test that `AbstractMvNormal` implements the expected API. -""" -function test_mvnormal( - g::AbstractMvNormal, n_tsamples::Int=10^6, rng::Union{AbstractRNG, Nothing}=nothing -) - d = length(g) - μ = mean(g) - Σ = cov(g) - @test length(μ) == d - @test size(Σ) == (d, d) - @test var(g) ≈ diag(Σ) - @test entropy(g) ≈ 0.5 * logdet(2π * ℯ * Σ) - ldcov = logdetcov(g) - @test ldcov ≈ logdet(Σ) - vs = diag(Σ) - @test g == typeof(g)(params(g)...) - @test g == deepcopy(g) - @test minimum(g) == fill(-Inf, d) - @test maximum(g) == fill(Inf, d) - @test extrema(g) == (minimum(g), maximum(g)) - @test isless(extrema(g)...) - - # test sampling for AbstractMatrix (here, a SubArray): - subX = view(__rand(rng, d, 2d), :, 1:d) - @test isa(__rand!(rng, g, subX), SubArray) - - # sampling - @test isa(__rand(rng, g), Vector{Float64}) - X = __rand(rng, g, n_tsamples) - emp_mu = vec(mean(X, dims=2)) - Z = X .- emp_mu - emp_cov = (Z * Z') * inv(n_tsamples) - - mean_atols = 8 .* sqrt.(vs ./ n_tsamples) - cov_atols = 10 .* sqrt.(vs .* vs') ./ sqrt.(n_tsamples) - for i = 1:d - @test isapprox(emp_mu[i], μ[i], atol=mean_atols[i]) - end - for i = 1:d, j = 1:d - @test isapprox(emp_cov[i,j], Σ[i,j], atol=cov_atols[i,j]) +import ..Distributions + +function test_mvnormal end + +if isdefined(Base, :get_extension) && isdefined(Base.Experimental, :register_error_hint) + function __init__() + # Better error message if users forget to load Test + Base.Experimental.register_error_hint(MethodError) do io, exc, _, _ + if exc.f === test_mvnormal && + (Base.get_extension(Distributions, :DistributionsTestExt) === nothing) + print(io, "\nDid you forget to load Test?") + end + end end - - X = rand(MersenneTwister(14), g, n_tsamples) - Y = rand(MersenneTwister(14), g, n_tsamples) - @test X == Y - emp_mu = vec(mean(X, dims=2)) - Z = X .- emp_mu - emp_cov = (Z * Z') * inv(n_tsamples) - for i = 1:d - @test isapprox(emp_mu[i] , μ[i] , atol=mean_atols[i]) - end - for i = 1:d, j = 1:d - @test isapprox(emp_cov[i,j], Σ[i,j], atol=cov_atols[i,j]) - end - - - # evaluation of sqmahal & logpdf - U = X .- μ - sqm = vec(sum(U .* (Σ \ U), dims=1)) - for i = 1:min(100, n_tsamples) - @test sqmahal(g, X[:,i]) ≈ sqm[i] - end - @test sqmahal(g, X) ≈ sqm - - lp = -0.5 .* sqm .- 0.5 * (d * log(2.0 * pi) + ldcov) - for i = 1:min(100, n_tsamples) - @test logpdf(g, X[:,i]) ≈ lp[i] - end - @test logpdf(g, X) ≈ lp - - # log likelihood - @test loglikelihood(g, X) ≈ sum(i -> Distributions._logpdf(g, X[:,i]), 1:n_tsamples) - @test loglikelihood(g, X[:, 1]) ≈ logpdf(g, X[:, 1]) - @test loglikelihood(g, [X[:, i] for i in axes(X, 2)]) ≈ loglikelihood(g, X) end end From 277f0ab8d68c70101274970a80b6af10156b1d5a Mon Sep 17 00:00:00 2001 From: "github-actions[bot]" <41898282+github-actions[bot]@users.noreply.github.com> Date: Fri, 3 Nov 2023 15:15:58 +0100 Subject: [PATCH 07/26] CompatHelper: add new compat entry for Statistics at version 1, (keep existing compat) (#1793) Co-authored-by: CompatHelper Julia --- Project.toml | 1 + 1 file changed, 1 insertion(+) diff --git a/Project.toml b/Project.toml index cf6ee221ca..4babf7a402 100644 --- a/Project.toml +++ b/Project.toml @@ -36,6 +36,7 @@ FillArrays = "0.9, 0.10, 0.11, 0.12, 0.13, 1" PDMats = "0.10, 0.11" QuadGK = "2" SpecialFunctions = "1.2, 2" +Statistics = "1" StatsAPI = "1.6" StatsBase = "0.32, 0.33, 0.34" StatsFuns = "0.9.15, 1" From 7e232ca86c56c4b2e1e7a44b16aa23c26668af3d Mon Sep 17 00:00:00 2001 From: David Widmann Date: Fri, 3 Nov 2023 20:09:09 +0100 Subject: [PATCH 08/26] Add compat entries for stdlibs (#1794) * Add compat entries for stdlibs * Update Project.toml * Update Project.toml * Update Project.toml * Update Project.toml --- Project.toml | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/Project.toml b/Project.toml index 4babf7a402..68b96caa8f 100644 --- a/Project.toml +++ b/Project.toml @@ -33,13 +33,17 @@ DistributionsTestExt = "Test" ChainRulesCore = "1" DensityInterface = "0.4" FillArrays = "0.9, 0.10, 0.11, 0.12, 0.13, 1" +LinearAlgebra = "<0.0.1, 1" PDMats = "0.10, 0.11" +Printf = "<0.0.1, 1" QuadGK = "2" +Random = "<0.0.1, 1" SpecialFunctions = "1.2, 2" Statistics = "1" StatsAPI = "1.6" StatsBase = "0.32, 0.33, 0.34" StatsFuns = "0.9.15, 1" +Test = "<0.0.1, 1" julia = "1.3" [extras] From 40f40bcdc4ea7d476a2c5fa664957df658e3df73 Mon Sep 17 00:00:00 2001 From: David Widmann Date: Fri, 8 Dec 2023 12:04:41 +0100 Subject: [PATCH 09/26] Fix method ambiguities & unbound parameters + add Aqua tests (#1804) * Fix method ambiguities & unbound parameters + add Aqua tests * Fix compat entries for Julia 1.3 * Update aqua.jl --- Project.toml | 14 +++- README.md | 1 + src/censored.jl | 1 + src/common.jl | 107 ++++++++++++------------- src/deprecates.jl | 6 +- src/mixtures/mixturemodel.jl | 16 ++-- src/multivariate/mvtdist.jl | 2 +- src/product.jl | 14 ++-- src/univariate/discrete/categorical.jl | 2 +- src/univariates.jl | 6 +- test/aqua.jl | 21 +++++ test/runtests.jl | 6 +- 12 files changed, 112 insertions(+), 84 deletions(-) create mode 100644 test/aqua.jl diff --git a/Project.toml b/Project.toml index 68b96caa8f..df9f7004cc 100644 --- a/Project.toml +++ b/Project.toml @@ -30,15 +30,26 @@ DistributionsDensityInterfaceExt = "DensityInterface" DistributionsTestExt = "Test" [compat] +Aqua = "0.8" +Calculus = "0.5" ChainRulesCore = "1" +ChainRulesTestUtils = "1" DensityInterface = "0.4" +Distributed = "<0.0.1, 1" FillArrays = "0.9, 0.10, 0.11, 0.12, 0.13, 1" +FiniteDifferences = "0.12" +ForwardDiff = "0.10" +JSON = "0.21" LinearAlgebra = "<0.0.1, 1" +OffsetArrays = "1" PDMats = "0.10, 0.11" Printf = "<0.0.1, 1" QuadGK = "2" Random = "<0.0.1, 1" +SparseArrays = "<0.0.1, 1" SpecialFunctions = "1.2, 2" +StableRNGs = "1" +StaticArrays = "1" Statistics = "1" StatsAPI = "1.6" StatsBase = "0.32, 0.33, 0.34" @@ -47,6 +58,7 @@ Test = "<0.0.1, 1" julia = "1.3" [extras] +Aqua = "4c88cf16-eb10-579e-8560-4a9242c79595" Calculus = "49dc2e85-a5d0-5ad3-a950-438e2897f1b9" ChainRulesCore = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4" ChainRulesTestUtils = "cdddcdb0-9152-4a09-a978-84456f9df70a" @@ -62,4 +74,4 @@ StaticArrays = "90137ffa-7385-5640-81b9-e52037218182" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" [targets] -test = ["StableRNGs", "Calculus", "ChainRulesCore", "ChainRulesTestUtils", "DensityInterface", "Distributed", "FiniteDifferences", "ForwardDiff", "JSON", "SparseArrays", "StaticArrays", "Test", "OffsetArrays"] +test = ["Aqua", "StableRNGs", "Calculus", "ChainRulesCore", "ChainRulesTestUtils", "DensityInterface", "Distributed", "FiniteDifferences", "ForwardDiff", "JSON", "SparseArrays", "StaticArrays", "Test", "OffsetArrays"] diff --git a/README.md b/README.md index 80a39ecb72..0071cf3617 100644 --- a/README.md +++ b/README.md @@ -4,6 +4,7 @@ Distributions.jl [![Build Status](https://github.com/JuliaStats/Distributions.jl/workflows/CI/badge.svg)](https://github.com/JuliaStats/Distributions.jl/actions) [![](https://zenodo.org/badge/DOI/10.5281/zenodo.2647458.svg)](https://zenodo.org/record/2647458) [![Coverage Status](https://coveralls.io/repos/JuliaStats/Distributions.jl/badge.svg?branch=master)](https://coveralls.io/r/JuliaStats/Distributions.jl?branch=master) +[![Aqua QA](https://raw.githubusercontent.com/JuliaTesting/Aqua.jl/master/badge.svg)](https://github.com/JuliaTesting/Aqua.jl) [![](https://img.shields.io/badge/docs-latest-blue.svg)](https://JuliaStats.github.io/Distributions.jl/latest/) [![](https://img.shields.io/badge/docs-stable-blue.svg)](https://JuliaStats.github.io/Distributions.jl/stable/) diff --git a/src/censored.jl b/src/censored.jl index 2d5e1d67fa..5b8cfac5ee 100644 --- a/src/censored.jl +++ b/src/censored.jl @@ -431,6 +431,7 @@ _in_open_interval(x::Real, l::Real, ::Nothing) = x > l _clamp(x, l, u) = clamp(x, l, u) _clamp(x, ::Nothing, u) = min(x, u) _clamp(x, l, ::Nothing) = max(x, l) +_clamp(x, ::Nothing, u::Nothing) = x _to_truncated(d::Censored) = truncated(d.uncensored, d.lower, d.upper) diff --git a/src/common.jl b/src/common.jl index 703e126930..31a218e197 100644 --- a/src/common.jl +++ b/src/common.jl @@ -212,13 +212,25 @@ usually it is sufficient to implement `logpdf`. See also: [`logpdf`](@ref). """ @inline function pdf( - d::Distribution{ArrayLikeVariate{N}}, x::AbstractArray{<:Real,N} -) where {N} - @boundscheck begin - size(x) == size(d) || - throw(DimensionMismatch("inconsistent array dimensions")) + d::Distribution{ArrayLikeVariate{N}}, x::AbstractArray{<:Real,M} +) where {N,M} + if M == N + @boundscheck begin + size(x) == size(d) || + throw(DimensionMismatch("inconsistent array dimensions")) + end + return _pdf(d, x) + else + @boundscheck begin + M > N || + throw(DimensionMismatch( + "number of dimensions of the variates ($M) must be greater than or equal to the dimension of the distribution ($N)" + )) + ntuple(i -> size(x, i), Val(N)) == size(d) || + throw(DimensionMismatch("inconsistent array dimensions")) + end + return @inbounds map(Base.Fix1(pdf, d), eachvariate(x, variate_form(typeof(d)))) end - return _pdf(d, x) end function _pdf(d::Distribution{ArrayLikeVariate{N}}, x::AbstractArray{<:Real,N}) where {N} @@ -241,13 +253,25 @@ size of `x`. See also: [`pdf`](@ref). """ @inline function logpdf( - d::Distribution{ArrayLikeVariate{N}}, x::AbstractArray{<:Real,N} -) where {N} - @boundscheck begin - size(x) == size(d) || - throw(DimensionMismatch("inconsistent array dimensions")) + d::Distribution{ArrayLikeVariate{N}}, x::AbstractArray{<:Real,M} +) where {N,M} + if M == N + @boundscheck begin + size(x) == size(d) || + throw(DimensionMismatch("inconsistent array dimensions")) + end + return _logpdf(d, x) + else + @boundscheck begin + M > N || + throw(DimensionMismatch( + "number of dimensions of the variates ($M) must be greater than or equal to the dimension of the distribution ($N)" + )) + ntuple(i -> size(x, i), Val(N)) == size(d) || + throw(DimensionMismatch("inconsistent array dimensions")) + end + return @inbounds map(Base.Fix1(logpdf, d), eachvariate(x, variate_form(typeof(d)))) end - return _logpdf(d, x) end # `_logpdf` should be implemented and has no default definition @@ -272,20 +296,6 @@ Base.@propagate_inbounds function pdf( return map(Base.Fix1(pdf, d), x) end -@inline function pdf( - d::Distribution{ArrayLikeVariate{N}}, x::AbstractArray{<:Real,M}, -) where {N,M} - @boundscheck begin - M > N || - throw(DimensionMismatch( - "number of dimensions of `x` ($M) must be greater than number of dimensions of `d` ($N)" - )) - ntuple(i -> size(x, i), Val(N)) == size(d) || - throw(DimensionMismatch("inconsistent array dimensions")) - end - return @inbounds map(Base.Fix1(pdf, d), eachvariate(x, variate_form(typeof(d)))) -end - """ logpdf(d::Distribution{ArrayLikeVariate{N}}, x) where {N} @@ -305,20 +315,6 @@ Base.@propagate_inbounds function logpdf( return map(Base.Fix1(logpdf, d), x) end -@inline function logpdf( - d::Distribution{ArrayLikeVariate{N}}, x::AbstractArray{<:Real,M}, -) where {N,M} - @boundscheck begin - M > N || - throw(DimensionMismatch( - "number of dimensions of `x` ($M) must be greater than number of dimensions of `d` ($N)" - )) - ntuple(i -> size(x, i), Val(N)) == size(d) || - throw(DimensionMismatch("inconsistent array dimensions")) - end - return @inbounds map(Base.Fix1(logpdf, d), eachvariate(x, variate_form(typeof(d)))) -end - """ pdf!(out, d::Distribution{ArrayLikeVariate{N}}, x) where {N} @@ -365,7 +361,7 @@ end @boundscheck begin M > N || throw(DimensionMismatch( - "number of dimensions of `x` ($M) must be greater than number of dimensions of `d` ($N)" + "number of dimensions of the variates ($M) must be greater than the dimension of the distribution ($N)" )) ntuple(i -> size(x, i), Val(N)) == size(d) || throw(DimensionMismatch("inconsistent array dimensions")) @@ -414,7 +410,7 @@ See also: [`pdf!`](@ref). @boundscheck begin M > N || throw(DimensionMismatch( - "number of dimensions of `x` ($M) must be greater than number of dimensions of `d` ($N)" + "number of dimensions of the variates ($M) must be greater than the dimension of the distribution ($N)" )) ntuple(i -> size(x, i), Val(N)) == size(d) || throw(DimensionMismatch("inconsistent array dimensions")) @@ -445,23 +441,22 @@ be - an array of dimension `N + 1` with `size(x)[1:N] == size(d)`, or - an array of arrays `xi` of dimension `N` with `size(xi) == size(d)`. """ -Base.@propagate_inbounds function loglikelihood( - d::Distribution{ArrayLikeVariate{N}}, x::AbstractArray{<:Real,N}, -) where {N} - return logpdf(d, x) -end -@inline function loglikelihood( +Base.@propagate_inbounds @inline function loglikelihood( d::Distribution{ArrayLikeVariate{N}}, x::AbstractArray{<:Real,M}, ) where {N,M} - @boundscheck begin - M > N || - throw(DimensionMismatch( - "number of dimensions of `x` ($M) must be greater than number of dimensions of `d` ($N)" - )) - ntuple(i -> size(x, i), Val(N)) == size(d) || - throw(DimensionMismatch("inconsistent array dimensions")) + if M == N + return logpdf(d, x) + else + @boundscheck begin + M > N || + throw(DimensionMismatch( + "number of dimensions of the variates ($M) must be greater than or equal to the dimension of the distribution ($N)" + )) + ntuple(i -> size(x, i), Val(N)) == size(d) || + throw(DimensionMismatch("inconsistent array dimensions")) + end + return @inbounds sum(Base.Fix1(logpdf, d), eachvariate(x, ArrayLikeVariate{N})) end - return @inbounds sum(Base.Fix1(logpdf, d), eachvariate(x, ArrayLikeVariate{N})) end Base.@propagate_inbounds function loglikelihood( d::Distribution{ArrayLikeVariate{N}}, x::AbstractArray{<:AbstractArray{<:Real,N}}, diff --git a/src/deprecates.jl b/src/deprecates.jl index e779d6dbb0..cacaa12b59 100644 --- a/src/deprecates.jl +++ b/src/deprecates.jl @@ -34,9 +34,9 @@ for fun in [:pdf, :logpdf, fun! = Symbol(fun, '!') @eval begin - @deprecate ($_fun!)(r::AbstractArray, d::UnivariateDistribution, X::AbstractArray) r .= ($fun).(d, X) false - @deprecate ($fun!)(r::AbstractArray, d::UnivariateDistribution, X::AbstractArray) r .= ($fun).(d, X) false - @deprecate ($fun)(d::UnivariateDistribution, X::AbstractArray) ($fun).(d, X) + @deprecate ($_fun!)(r::AbstractArray{<:Real}, d::UnivariateDistribution, X::AbstractArray{<:Real}) r .= ($fun).(d, X) false + @deprecate ($fun!)(r::AbstractArray{<:Real}, d::UnivariateDistribution, X::AbstractArray{<:Real}) r .= ($fun).(d, X) false + @deprecate ($fun)(d::UnivariateDistribution, X::AbstractArray{<:Real}) ($fun).(d, X) end end diff --git a/src/mixtures/mixturemodel.jl b/src/mixtures/mixturemodel.jl index 7f5f6782e3..fc05012d0d 100644 --- a/src/mixtures/mixturemodel.jl +++ b/src/mixtures/mixturemodel.jl @@ -362,14 +362,14 @@ end pdf(d::UnivariateMixture, x::Real) = _mixpdf1(d, x) logpdf(d::UnivariateMixture, x::Real) = _mixlogpdf1(d, x) -_pdf!(r::AbstractArray, d::UnivariateMixture{Discrete}, x::UnitRange) = _mixpdf!(r, d, x) -_pdf!(r::AbstractArray, d::UnivariateMixture, x::AbstractArray) = _mixpdf!(r, d, x) -_logpdf!(r::AbstractArray, d::UnivariateMixture, x::AbstractArray) = _mixlogpdf!(r, d, x) - -_pdf(d::MultivariateMixture, x::AbstractVector) = _mixpdf1(d, x) -_logpdf(d::MultivariateMixture, x::AbstractVector) = _mixlogpdf1(d, x) -_pdf!(r::AbstractArray, d::MultivariateMixture, x::AbstractMatrix) = _mixpdf!(r, d, x) -_logpdf!(r::AbstractArray, d::MultivariateMixture, x::AbstractMatrix) = _mixlogpdf!(r, d, x) +_pdf!(r::AbstractArray{<:Real}, d::UnivariateMixture{Discrete}, x::UnitRange) = _mixpdf!(r, d, x) +_pdf!(r::AbstractArray{<:Real}, d::UnivariateMixture, x::AbstractArray{<:Real}) = _mixpdf!(r, d, x) +_logpdf!(r::AbstractArray{<:Real}, d::UnivariateMixture, x::AbstractArray{<:Real}) = _mixlogpdf!(r, d, x) + +_pdf(d::MultivariateMixture, x::AbstractVector{<:Real}) = _mixpdf1(d, x) +_logpdf(d::MultivariateMixture, x::AbstractVector{<:Real}) = _mixlogpdf1(d, x) +_pdf!(r::AbstractArray{<:Real}, d::MultivariateMixture, x::AbstractMatrix{<:Real}) = _mixpdf!(r, d, x) +_logpdf!(r::AbstractArray{<:Real}, d::MultivariateMixture, x::AbstractMatrix{<:Real}) = _mixlogpdf!(r, d, x) ## component-wise pdf and logpdf diff --git a/src/multivariate/mvtdist.jl b/src/multivariate/mvtdist.jl index dbd13ce60f..9076c364b5 100644 --- a/src/multivariate/mvtdist.jl +++ b/src/multivariate/mvtdist.jl @@ -139,7 +139,7 @@ function _logpdf(d::AbstractMvTDist, x::AbstractVector{T}) where T<:Real v - shdfhdim * log1p(sqmahal(d, x) / d.df) end -function _logpdf!(r::AbstractArray, d::AbstractMvTDist, x::AbstractMatrix) +function _logpdf!(r::AbstractArray{<:Real}, d::AbstractMvTDist, x::AbstractMatrix{<:Real}) sqmahal!(r, d, x) shdfhdim, v = mvtdist_consts(d) for i = 1:size(x, 2) diff --git a/src/product.jl b/src/product.jl index 7a4904ae7a..d6f7aaa6bf 100644 --- a/src/product.jl +++ b/src/product.jl @@ -12,7 +12,9 @@ struct ProductDistribution{N,M,D,S<:ValueSupport,T} <: Distribution{ArrayLikeVar size::Dims{N} function ProductDistribution{N,M,D}(dists::D) where {N,M,D} - isempty(dists) && error("product distribution must consist of at least one distribution") + if isempty(dists) + throw(ArgumentError("a product distribution must consist of at least one distribution")) + end return new{N,M,D,_product_valuesupport(dists),_product_eltype(dists)}( dists, _product_size(dists), @@ -20,11 +22,11 @@ struct ProductDistribution{N,M,D,S<:ValueSupport,T} <: Distribution{ArrayLikeVar end end -function ProductDistribution(dists::AbstractArray{<:Distribution{ArrayLikeVariate{M}},N}) where {M,N} +function ProductDistribution(dists::AbstractArray{<:Distribution{<:ArrayLikeVariate{M}},N}) where {M,N} return ProductDistribution{M + N,M,typeof(dists)}(dists) end -function ProductDistribution(dists::NTuple{N,Distribution{ArrayLikeVariate{M}}}) where {M,N} +function ProductDistribution(dists::Tuple{Distribution{<:ArrayLikeVariate{M}},Vararg{Distribution{<:ArrayLikeVariate{M}}}}) where {M} return ProductDistribution{M + 1,M,typeof(dists)}(dists) end @@ -54,10 +56,10 @@ function _product_size(dists::AbstractArray{<:Distribution{<:ArrayLikeVariate{M} size_dists = size(dists) return ntuple(i -> i <= M ? size_d[i] : size_dists[i-M], Val(M + N)) end -function _product_size(dists::NTuple{N,Distribution{<:ArrayLikeVariate{M}}}) where {M,N} +function _product_size(dists::Tuple{Distribution{<:ArrayLikeVariate{M}},Vararg{Distribution{<:ArrayLikeVariate{M}}, N}}) where {M,N} size_d = size(first(dists)) all(size(d) == size_d for d in dists) || error("all distributions must be of the same size") - return ntuple(i -> i <= M ? size_d[i] : N, Val(M + 1)) + return ntuple(i -> i <= M ? size_d[i] : N + 1, Val(M + 1)) end ## aliases @@ -167,7 +169,7 @@ function _rand!( end # `_logpdf` for arrays of distributions -# we have to fix a method ambiguity +# we have to fix some method ambiguities _logpdf(d::ProductDistribution{N}, x::AbstractArray{<:Real,N}) where {N} = __logpdf(d, x) _logpdf(d::ProductDistribution{2}, x::AbstractMatrix{<:Real}) = __logpdf(d, x) function __logpdf( diff --git a/src/univariate/discrete/categorical.jl b/src/univariate/discrete/categorical.jl index 9572acb17e..1ac338f6f0 100644 --- a/src/univariate/discrete/categorical.jl +++ b/src/univariate/discrete/categorical.jl @@ -89,7 +89,7 @@ function pdf(d::Categorical, x::Real) return insupport(d, x) ? ps[round(Int, x)] : zero(eltype(ps)) end -function _pdf!(r::AbstractArray, d::Categorical{T}, rgn::UnitRange) where {T<:Real} +function _pdf!(r::AbstractArray{<:Real}, d::Categorical{T}, rgn::UnitRange) where {T<:Real} vfirst = round(Int, first(rgn)) vlast = round(Int, last(rgn)) vl = max(vfirst, 1) diff --git a/src/univariates.jl b/src/univariates.jl index 0e2ae05e32..ca83d727e5 100644 --- a/src/univariates.jl +++ b/src/univariates.jl @@ -311,7 +311,7 @@ See also: [`logpdf`](@ref). pdf(d::UnivariateDistribution, x::Real) = exp(logpdf(d, x)) # extract value from array of zero dimension -_pdf(d::UnivariateDistribution, x::AbstractArray{<:Real,0}) = pdf(d, first(x)) +pdf(d::UnivariateDistribution, x::AbstractArray{<:Real,0}) = pdf(d, first(x)) """ logpdf(d::UnivariateDistribution, x::Real) @@ -323,7 +323,7 @@ See also: [`pdf`](@ref). logpdf(d::UnivariateDistribution, x::Real) # extract value from array of zero dimension -_logpdf(d::UnivariateDistribution, x::AbstractArray{<:Real,0}) = logpdf(d, first(x)) +logpdf(d::UnivariateDistribution, x::AbstractArray{<:Real,0}) = logpdf(d, first(x)) # loglikelihood for `Real` Base.@propagate_inbounds loglikelihood(d::UnivariateDistribution, x::Real) = logpdf(d, x) @@ -452,7 +452,7 @@ function _pdf_fill_outside!(r::AbstractArray, d::DiscreteUnivariateDistribution, return vl, vr, vfirst, vlast end -function _pdf!(r::AbstractArray, d::DiscreteUnivariateDistribution, X::UnitRange) +function _pdf!(r::AbstractArray{<:Real}, d::DiscreteUnivariateDistribution, X::UnitRange) vl,vr, vfirst, vlast = _pdf_fill_outside!(r, d, X) # fill central part: with non-zero pdf diff --git a/test/aqua.jl b/test/aqua.jl new file mode 100644 index 0000000000..4dd8c07064 --- /dev/null +++ b/test/aqua.jl @@ -0,0 +1,21 @@ +using Distributions +using Test + +import Aqua + +@testset "Aqua" begin + # Test ambiguities separately without Base and Core + # Ref: https://github.com/JuliaTesting/Aqua.jl/issues/77 + Aqua.test_all( + Distributions; + ambiguities = false, + # On older Julia versions, installed dependencies are quite old + # Thus unbound type parameters show up that are fixed in newer versions + unbound_args = VERSION >= v"1.6", + ) + # Tests are not reliable on older Julia versions and + # show ambiguities in loaded packages + if VERSION >= v"1.6" + Aqua.test_ambiguities(Distributions) + end +end diff --git a/test/runtests.jl b/test/runtests.jl index ce3f16b799..583132c536 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -11,6 +11,7 @@ import JSON import ForwardDiff const tests = [ + "aqua", "univariate/continuous/loguniform", "univariate/continuous/arcsine", "univariate/discrete/dirac", @@ -174,8 +175,3 @@ include("testutils.jl") include("$t.jl") end end - -# print method ambiguities -println("Potentially stale exports: ") -display(Test.detect_ambiguities(Distributions)) -println() From 8e231ed369efd80ea6063c353929e82ee10fb36e Mon Sep 17 00:00:00 2001 From: David Widmann Date: Fri, 8 Dec 2023 12:05:04 +0100 Subject: [PATCH 10/26] Update Project.toml --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index df9f7004cc..4e71ee9135 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "Distributions" uuid = "31c24e10-a181-5473-b8eb-7969acd0382f" authors = ["JuliaStats"] -version = "0.25.103" +version = "0.25.104" [deps] ChainRulesCore = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4" From b2d765246e0f35473fe354dc578a3f7edf11958d Mon Sep 17 00:00:00 2001 From: "github-actions[bot]" <41898282+github-actions[bot]@users.noreply.github.com> Date: Thu, 14 Dec 2023 09:12:05 +0100 Subject: [PATCH 11/26] CompatHelper: bump compat for GR to 0.73 for package docs, (keep existing compat) (#1811) Co-authored-by: CompatHelper Julia --- docs/Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/docs/Project.toml b/docs/Project.toml index 540725c6d6..58c4b62081 100644 --- a/docs/Project.toml +++ b/docs/Project.toml @@ -4,4 +4,4 @@ GR = "28b8d3ca-fb5f-59d9-8090-bfdbd6d07a71" [compat] Documenter = "0.26, 0.27" -GR = "0.72.1" +GR = "0.72.1, 0.73" From 9e72f1ff39eb547cc420a278b5d96323a6de5ad5 Mon Sep 17 00:00:00 2001 From: Seth Axen Date: Wed, 3 Jan 2024 11:53:07 +0100 Subject: [PATCH 12/26] Clarify definitions of inverse CDF functions (#1814) * Clarify definitions of inverse CDF functions * Such that to for which * Wrap "generalized" in parentheses --- src/univariates.jl | 18 ++++++++++++++---- 1 file changed, 14 insertions(+), 4 deletions(-) diff --git a/src/univariates.jl b/src/univariates.jl index ca83d727e5..b60e5a2949 100644 --- a/src/univariates.jl +++ b/src/univariates.jl @@ -182,7 +182,8 @@ std(d::UnivariateDistribution) = sqrt(var(d)) """ median(d::UnivariateDistribution) -Return the median value of distribution `d`. The median is the smallest `x` such that `cdf(d, x) ≥ 1/2`. +Return the median value of distribution `d`. The median is the smallest `x` in the support +of `d` for which `cdf(d, x) ≥ 1/2`. Corresponding to this definition as 1/2-quantile, a fallback is provided calling the `quantile` function. """ median(d::UnivariateDistribution) = quantile(d, 1//2) @@ -381,7 +382,10 @@ logccdf(d::UnivariateDistribution, x::Real) = log(ccdf(d, x)) """ quantile(d::UnivariateDistribution, q::Real) -Evaluate the inverse cumulative distribution function at `q`. +Evaluate the (generalized) inverse cumulative distribution function at `q`. + +For a given `0 ≤ q ≤ 1`, `quantile(d, q)` is the smallest value `x` in the support of `d` +for which `cdf(d, x) ≥ q`. See also: [`cquantile`](@ref), [`invlogcdf`](@ref), and [`invlogccdf`](@ref). """ @@ -397,14 +401,20 @@ cquantile(d::UnivariateDistribution, p::Real) = quantile(d, 1.0 - p) """ invlogcdf(d::UnivariateDistribution, lp::Real) -The inverse function of logcdf. +The (generalized) inverse function of [`logcdf`](@ref). + +For a given `lp ≤ 0`, `invlogcdf(d, lp)` is the smallest value `x` in the support of `d` for +which `logcdf(d, x) ≥ lp`. """ invlogcdf(d::UnivariateDistribution, lp::Real) = quantile(d, exp(lp)) """ invlogccdf(d::UnivariateDistribution, lp::Real) -The inverse function of logccdf. +The (generalized) inverse function of [`logccdf`](@ref). + +For a given `lp ≤ 0`, `invlogccdf(d, lp)` is the smallest value `x` in the support of `d` +for which `logccdf(d, x) ≤ lp`. """ invlogccdf(d::UnivariateDistribution, lp::Real) = quantile(d, -expm1(lp)) From 28bf738139dd4dc5bca69671f699912aa978d5a8 Mon Sep 17 00:00:00 2001 From: David Widmann Date: Fri, 5 Jan 2024 11:46:00 +0100 Subject: [PATCH 13/26] Fix `eachvariate` with zero variates (#1819) --- Project.toml | 2 +- src/eachvariate.jl | 2 +- test/eachvariate.jl | 27 +++++++++++++++++++++++++++ 3 files changed, 29 insertions(+), 2 deletions(-) diff --git a/Project.toml b/Project.toml index 4e71ee9135..60ae8b2618 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "Distributions" uuid = "31c24e10-a181-5473-b8eb-7969acd0382f" authors = ["JuliaStats"] -version = "0.25.104" +version = "0.25.105" [deps] ChainRulesCore = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4" diff --git a/src/eachvariate.jl b/src/eachvariate.jl index 701be99faa..583600f36d 100644 --- a/src/eachvariate.jl +++ b/src/eachvariate.jl @@ -7,7 +7,7 @@ end function EachVariate{V}(x::AbstractArray{<:Real,M}) where {V,M} ax = ntuple(i -> axes(x, i + V), Val(M - V)) - T = typeof(view(x, ntuple(i -> i <= V ? Colon() : firstindex(x, i), Val(M))...)) + T = Base.promote_op(view, typeof(x), ntuple(i -> i <= V ? Colon : eltype(axes(x, i)), Val(M))...) return EachVariate{V,typeof(x),typeof(ax),T,M-V}(x, ax) end diff --git a/test/eachvariate.jl b/test/eachvariate.jl index f41a5207d2..60cd1f68d3 100644 --- a/test/eachvariate.jl +++ b/test/eachvariate.jl @@ -1,12 +1,24 @@ +using Distributions using ChainRulesTestUtils using ChainRulesTestUtils: FiniteDifferences +using Random +using Test + # Without this, `to_vec` will also include the `axes` field of `EachVariate`. function FiniteDifferences.to_vec(xs::Distributions.EachVariate{V}) where {V} vals, vals_from_vec = FiniteDifferences.to_vec(xs.parent) return vals, x -> Distributions.EachVariate{V}(vals_from_vec(x)) end +# MWE in #1817 +struct FooEachvariate <: Sampleable{Multivariate, Continuous} end +Base.length(::FooEachvariate) = 3 +Base.eltype(::FooEachvariate) = Float64 +function Distributions._rand!(rng::AbstractRNG, ::FooEachvariate, x::AbstractVector{<:Real}) + return rand!(rng, x) +end + @testset "eachvariate.jl" begin @testset "ChainRules" begin xs = randn(2, 3, 4, 5) @@ -14,4 +26,19 @@ end test_rrule(Distributions.EachVariate{2}, xs) test_rrule(Distributions.EachVariate{3}, xs) end + + @testset "No variates (#1817)" begin + @test size(Distributions.eachvariate(rand(0), Univariate)) == (0,) + @test size(Distributions.eachvariate(rand(3, 0, 1), Univariate)) == (3, 0, 1) + @test size(Distributions.eachvariate(rand(3, 2, 0), Univariate)) == (3, 2, 0) + + @test size(Distributions.eachvariate(rand(4, 0), Multivariate)) == (0,) + @test size(Distributions.eachvariate(rand(4, 0, 2), Multivariate)) == (0, 2) + @test size(Distributions.eachvariate(rand(4, 5, 0), Multivariate)) == (5, 0) + @test size(Distributions.eachvariate(rand(4, 5, 0, 2), Multivariate)) == (5, 0, 2) + + draws = @inferred(rand(FooEachvariate(), 0)) + @test draws isa Matrix{Float64} + @test size(draws) == (3, 0) + end end From d58adb6088a56f6790c218389525ea7ee28d2b7b Mon Sep 17 00:00:00 2001 From: David Widmann Date: Fri, 5 Jan 2024 12:20:54 +0100 Subject: [PATCH 14/26] Replace broadcasting over distributions with broadcasting with partially applied functions (#1818) * Replace broadcasting over distributions with broadcasting with partially applied functions * Delete src/multivariate/mvnormal copy.jl --- src/deprecates.jl | 8 ++-- src/mixtures/mixturemodel.jl | 8 ++-- src/multivariate/jointorderstatistics.jl | 2 +- src/qq.jl | 4 +- src/univariate/continuous/uniform.jl | 2 +- src/univariate/discrete/betabinomial.jl | 13 ++++-- src/univariate/discrete/hypergeometric.jl | 2 +- .../discrete/noncentralhypergeometric.jl | 9 ++-- test/censored.jl | 8 ++-- test/matrixvariates.jl | 2 +- test/mixture.jl | 2 +- test/multivariate/product.jl | 3 +- test/product.jl | 2 +- test/testutils.jl | 46 +++++++++---------- test/univariate/continuous/johnsonsu.jl | 8 ++-- test/univariate/continuous/logitnormal.jl | 2 +- test/univariate/continuous/rician.jl | 10 ++-- test/univariate/continuous/skewnormal.jl | 2 +- test/univariate/discrete/binomial.jl | 4 +- test/univariate/discrete/categorical.jl | 4 +- test/univariate/discrete/poissonbinomial.jl | 6 +-- test/univariate/discrete/soliton.jl | 4 +- test/univariate/locationscale.jl | 4 +- test/univariates.jl | 4 +- 24 files changed, 82 insertions(+), 77 deletions(-) diff --git a/src/deprecates.jl b/src/deprecates.jl index cacaa12b59..3cfc20a5c6 100644 --- a/src/deprecates.jl +++ b/src/deprecates.jl @@ -34,13 +34,13 @@ for fun in [:pdf, :logpdf, fun! = Symbol(fun, '!') @eval begin - @deprecate ($_fun!)(r::AbstractArray{<:Real}, d::UnivariateDistribution, X::AbstractArray{<:Real}) r .= ($fun).(d, X) false - @deprecate ($fun!)(r::AbstractArray{<:Real}, d::UnivariateDistribution, X::AbstractArray{<:Real}) r .= ($fun).(d, X) false - @deprecate ($fun)(d::UnivariateDistribution, X::AbstractArray{<:Real}) ($fun).(d, X) + @deprecate ($_fun!)(r::AbstractArray{<:Real}, d::UnivariateDistribution, X::AbstractArray{<:Real}) r .= Base.Fix1($fun, d).(X) false + @deprecate ($fun!)(r::AbstractArray{<:Real}, d::UnivariateDistribution, X::AbstractArray{<:Real}) r .= Base.Fix1($fun, d).(X) false + @deprecate ($fun)(d::UnivariateDistribution, X::AbstractArray{<:Real}) map(Base.Fix1($fun, d), X) end end -@deprecate pdf(d::DiscreteUnivariateDistribution) pdf.(Ref(d), support(d)) +@deprecate pdf(d::DiscreteUnivariateDistribution) map(Base.Fix1(pdf, d), support(d)) # Wishart constructors @deprecate Wishart(df::Real, S::AbstractPDMat, warn::Bool) Wishart(df, S) diff --git a/src/mixtures/mixturemodel.jl b/src/mixtures/mixturemodel.jl index fc05012d0d..c3d3b1f919 100644 --- a/src/mixtures/mixturemodel.jl +++ b/src/mixtures/mixturemodel.jl @@ -297,7 +297,7 @@ function _mixpdf!(r::AbstractArray, d::AbstractMixtureModel, x) pi = p[i] if pi > 0.0 if d isa UnivariateMixture - t .= pdf.(component(d, i), x) + t .= Base.Fix1(pdf, component(d, i)).(x) else pdf!(t, component(d, i), x) end @@ -326,7 +326,7 @@ function _mixlogpdf!(r::AbstractArray, d::AbstractMixtureModel, x) lp_i = view(Lp, :, i) # compute logpdf in batch and store if d isa UnivariateMixture - lp_i .= logpdf.(component(d, i), x) + lp_i .= Base.Fix1(logpdf, component(d, i)).(x) else logpdf!(lp_i, component(d, i), x) end @@ -398,7 +398,7 @@ function _cwise_pdf!(r::AbstractMatrix, d::AbstractMixtureModel, X) size(r) == (n, K) || error("The size of r is incorrect.") for i = 1:K if d isa UnivariateMixture - view(r,:,i) .= pdf.(Ref(component(d, i)), X) + view(r,:,i) .= Base.Fix1(pdf, component(d, i)).(X) else pdf!(view(r,:,i),component(d, i), X) end @@ -412,7 +412,7 @@ function _cwise_logpdf!(r::AbstractMatrix, d::AbstractMixtureModel, X) size(r) == (n, K) || error("The size of r is incorrect.") for i = 1:K if d isa UnivariateMixture - view(r,:,i) .= logpdf.(Ref(component(d, i)), X) + view(r,:,i) .= Base.Fix1(logpdf, component(d, i)).(X) else logpdf!(view(r,:,i), component(d, i), X) end diff --git a/src/multivariate/jointorderstatistics.jl b/src/multivariate/jointorderstatistics.jl index 71c3f93bf3..1fbed0d1b6 100644 --- a/src/multivariate/jointorderstatistics.jl +++ b/src/multivariate/jointorderstatistics.jl @@ -162,7 +162,7 @@ function _rand!(rng::AbstractRNG, d::JointOrderStatistics, x::AbstractVector{<:R else s += randexp(rng, T) end - x .= quantile.(d.dist, x ./ s) + x .= Base.Fix1(quantile, d.dist).(x ./ s) end return x end diff --git a/src/qq.jl b/src/qq.jl index fe80e31a68..6431ed2bf9 100644 --- a/src/qq.jl +++ b/src/qq.jl @@ -35,14 +35,14 @@ function qqbuild(x::AbstractVector, d::UnivariateDistribution) n = length(x) grid = ppoints(n) qx = quantile(x, grid) - qd = quantile.(Ref(d), grid) + qd = map(Base.Fix1(quantile, d), grid) return QQPair(qx, qd) end function qqbuild(d::UnivariateDistribution, x::AbstractVector) n = length(x) grid = ppoints(n) - qd = quantile.(Ref(d), grid) + qd = map(Base.Fix1(quantile, d), grid) qx = quantile(x, grid) return QQPair(qd, qx) end diff --git a/src/univariate/continuous/uniform.jl b/src/univariate/continuous/uniform.jl index 1f535159d0..d4beca7d79 100644 --- a/src/univariate/continuous/uniform.jl +++ b/src/univariate/continuous/uniform.jl @@ -154,7 +154,7 @@ Base.:*(c::Real, d::Uniform) = Uniform(minmax(c * d.a, c * d.b)...) rand(rng::AbstractRNG, d::Uniform) = d.a + (d.b - d.a) * rand(rng) _rand!(rng::AbstractRNG, d::Uniform, A::AbstractArray{<:Real}) = - A .= quantile.(d, rand!(rng, A)) + A .= Base.Fix1(quantile, d).(rand!(rng, A)) #### Fitting diff --git a/src/univariate/discrete/betabinomial.jl b/src/univariate/discrete/betabinomial.jl index 618b3f83ef..7f2a0ac0dc 100644 --- a/src/univariate/discrete/betabinomial.jl +++ b/src/univariate/discrete/betabinomial.jl @@ -103,12 +103,15 @@ for f in (:ccdf, :logcdf, :logccdf) end end -entropy(d::BetaBinomial) = entropy(Categorical(pdf.(Ref(d),support(d)))) -median(d::BetaBinomial) = median(Categorical(pdf.(Ref(d),support(d)))) - 1 -mode(d::BetaBinomial) = argmax(pdf.(Ref(d),support(d))) - 1 -modes(d::BetaBinomial) = modes(Categorical(pdf.(Ref(d),support(d)))) .- 1 +# Shifted categorical distribution corresponding to `BetaBinomial` +_categorical(d::BetaBinomial) = Categorical(map(Base.Fix1(pdf, d), support(d))) -quantile(d::BetaBinomial, p::Float64) = quantile(Categorical(pdf.(Ref(d), support(d))), p) - 1 +entropy(d::BetaBinomial) = entropy(_categorical(d)) +median(d::BetaBinomial) = median(_categorical(d)) - 1 +mode(d::BetaBinomial) = mode(_categorical(d)) - 1 +modes(d::BetaBinomial) = modes(_categorical(d)) .- 1 + +quantile(d::BetaBinomial, p::Float64) = quantile(_categorical(d), p) - 1 #### Sampling diff --git a/src/univariate/discrete/hypergeometric.jl b/src/univariate/discrete/hypergeometric.jl index 9f72c2fea4..de575abdcb 100644 --- a/src/univariate/discrete/hypergeometric.jl +++ b/src/univariate/discrete/hypergeometric.jl @@ -75,7 +75,7 @@ function kurtosis(d::Hypergeometric) a/b end -entropy(d::Hypergeometric) = entropy(pdf.(Ref(d), support(d))) +entropy(d::Hypergeometric) = entropy(map(Base.Fix1(pdf, d), support(d))) ### Evaluation & Sampling diff --git a/src/univariate/discrete/noncentralhypergeometric.jl b/src/univariate/discrete/noncentralhypergeometric.jl index 4c33d26215..52853cfc08 100644 --- a/src/univariate/discrete/noncentralhypergeometric.jl +++ b/src/univariate/discrete/noncentralhypergeometric.jl @@ -256,9 +256,12 @@ end Base.convert(::Type{WalleniusNoncentralHypergeometric{T}}, d::WalleniusNoncentralHypergeometric{T}) where {T<:Real} = d # Properties -mean(d::WalleniusNoncentralHypergeometric) = sum(support(d) .* pdf.(Ref(d), support(d))) -var(d::WalleniusNoncentralHypergeometric) = sum((support(d) .- mean(d)).^2 .* pdf.(Ref(d), support(d))) -mode(d::WalleniusNoncentralHypergeometric) = support(d)[argmax(pdf.(Ref(d), support(d)))] +function _discretenonparametric(d::WalleniusNoncentralHypergeometric) + return DiscreteNonParametric(support(d), map(Base.Fix1(pdf, d), support(d))) +end +mean(d::WalleniusNoncentralHypergeometric) = mean(_discretenonparametric(d)) +var(d::WalleniusNoncentralHypergeometric) = var(_discretenonparametric(d)) +mode(d::WalleniusNoncentralHypergeometric) = mode(_discretenonparametric(d)) entropy(d::WalleniusNoncentralHypergeometric) = 1 diff --git a/test/censored.jl b/test/censored.jl index eaad72cdc2..27f30cadfc 100644 --- a/test/censored.jl +++ b/test/censored.jl @@ -207,7 +207,7 @@ end end @test @inferred(median(d)) ≈ clamp(median(d0), l, u) @inferred quantile(d, 0.5) - @test quantile.(d, 0:0.01:1) ≈ clamp.(quantile.(d0, 0:0.01:1), l, u) + @test Base.Fix1(quantile, d).(0:0.01:1) ≈ clamp.(Base.Fix1(quantile, d0).(0:0.01:1), l, u) # special-case pdf/logpdf/loglikelihood since when replacing Dirac(μ) with # Normal(μ, 0), they are infinite if lower === nothing || !isfinite(lower) @@ -253,7 +253,7 @@ end @test f(d) ≈ f(dmix) end @test median(d) ≈ clamp(median(d0), l, u) - @test quantile.(d, 0:0.01:1) ≈ clamp.(quantile.(d0, 0:0.01:1), l, u) + @test Base.Fix1(quantile, d).(0:0.01:1) ≈ clamp.(Base.Fix1(quantile, d0).(0:0.01:1), l, u) # special-case pdf/logpdf/loglikelihood since when replacing Dirac(μ) with # Normal(μ, 0), they are infinite if lower === nothing @@ -311,7 +311,7 @@ end end @test @inferred(median(d)) ≈ clamp(median(d0), l, u) @inferred quantile(d, 0.5) - @test quantile.(d, 0:0.01:1) ≈ clamp.(quantile.(d0, 0:0.01:1), l, u) + @test Base.Fix1(quantile, d).(0:0.01:1) ≈ clamp.(Base.Fix1(quantile, d0).(0:0.01:1), l, u) # rand x = rand(d, 10_000) @test all(x -> insupport(d, x), x) @@ -346,7 +346,7 @@ end @test f(d, 5) ≈ f(dmix, 5) end @test median(d) ≈ clamp(median(d0), l, u) - @test quantile.(d, 0:0.01:0.99) ≈ clamp.(quantile.(d0, 0:0.01:0.99), l, u) + @test Base.Fix1(quantile, d).(0:0.01:0.99) ≈ clamp.(Base.Fix1(quantile, d0).(0:0.01:0.99), l, u) x = rand(d, 100) @test loglikelihood(d, x) ≈ loglikelihood(dmix, x) # rand diff --git a/test/matrixvariates.jl b/test/matrixvariates.jl index 3d0dddfff6..18e984a470 100644 --- a/test/matrixvariates.jl +++ b/test/matrixvariates.jl @@ -117,7 +117,7 @@ function test_convert(d::MatrixDistribution) @test d == deepcopy(d) for elty in (Float32, Float64, BigFloat) del1 = convert(distname{elty}, d) - del2 = convert(distname{elty}, getfield.(Ref(d), fieldnames(typeof(d)))...) + del2 = convert(distname{elty}, (Base.Fix1(getfield, d)).(fieldnames(typeof(d)))...) @test del1 isa distname{elty} @test del2 isa distname{elty} @test partype(del1) == elty diff --git a/test/mixture.jl b/test/mixture.jl index 9b84f69326..0b25a2346a 100644 --- a/test/mixture.jl +++ b/test/mixture.jl @@ -40,7 +40,7 @@ function test_mixture(g::UnivariateMixture, n::Int, ns::Int, for i = 1:n @test @inferred(cdf(g, X[i])) ≈ cf[i] end - @test cdf.(g, X) ≈ cf + @test Base.Fix1(cdf, g).(X) ≈ cf # evaluation P0 = zeros(T, n, K) diff --git a/test/multivariate/product.jl b/test/multivariate/product.jl index 840eb409f2..c452f96f76 100644 --- a/test/multivariate/product.jl +++ b/test/multivariate/product.jl @@ -67,8 +67,7 @@ end for a in ([0, 1], [-0.5, 0.5]) # Construct independent distributions and `Product` distribution from these. - support = fill(a, N) - ds = DiscreteNonParametric.(support, Ref([0.5, 0.5])) + ds = [DiscreteNonParametric(copy(a), [0.5, 0.5]) for _ in 1:N] x = rand.(ds) d_product = product_distribution(ds) @test d_product isa Product diff --git a/test/product.jl b/test/product.jl index 829e80d1b8..16f12328d3 100644 --- a/test/product.jl +++ b/test/product.jl @@ -93,7 +93,7 @@ end for a in ([0, 1], [-0.5, 0.5]) # Construct independent distributions and `ProductDistribution` from these. - ds1 = DiscreteNonParametric.(fill(a, N), Ref([0.5, 0.5])) + ds1 = [DiscreteNonParametric(copy(a), [0.5, 0.5]) for _ in 1:N] # Replace with # d_product1 = @inferred(product_distribution(ds1)) # when `Product` is removed diff --git a/test/testutils.jl b/test/testutils.jl index ad109e58e3..2859856a73 100644 --- a/test/testutils.jl +++ b/test/testutils.jl @@ -128,7 +128,7 @@ function test_samples(s::Sampleable{Univariate, Discrete}, # the sampleable rmin = floor(Int,quantile(distr, 0.00001))::Int rmax = floor(Int,quantile(distr, 0.99999))::Int m = rmax - rmin + 1 # length of the range - p0 = pdf.(Ref(distr), rmin:rmax) # reference probability masses + p0 = map(Base.Fix1(pdf, distr), rmin:rmax) # reference probability masses @assert length(p0) == m # determine confidence intervals for counts: @@ -233,7 +233,7 @@ function test_samples(s::Sampleable{Univariate, Continuous}, # the sampleable # clb = Vector{Int}(undef, nbins) cub = Vector{Int}(undef, nbins) - cdfs = cdf.(Ref(distr), edges) + cdfs = map(Base.Fix1(cdf, distr), edges) for i = 1:nbins pi = cdfs[i+1] - cdfs[i] @@ -385,19 +385,19 @@ function test_range_evaluation(d::DiscreteUnivariateDistribution) rmin = round(Int, islowerbounded(d) ? vmin : quantile(d, 0.001))::Int rmax = round(Int, isupperbounded(d) ? vmax : quantile(d, 0.999))::Int - p0 = pdf.(Ref(d), collect(rmin:rmax)) - @test pdf.(Ref(d), rmin:rmax) ≈ p0 + p0 = map(Base.Fix1(pdf, d), collect(rmin:rmax)) + @test map(Base.Fix1(pdf, d), rmin:rmax) ≈ p0 if rmin + 2 <= rmax - @test pdf.(Ref(d), rmin+1:rmax-1) ≈ p0[2:end-1] + @test map(Base.Fix1(pdf, d), rmin+1:rmax-1) ≈ p0[2:end-1] end if isbounded(d) - @test pdf.(Ref(d), support(d)) ≈ p0 - @test pdf.(Ref(d), rmin-2:rmax) ≈ vcat(0.0, 0.0, p0) - @test pdf.(Ref(d), rmin:rmax+3) ≈ vcat(p0, 0.0, 0.0, 0.0) - @test pdf.(Ref(d), rmin-2:rmax+3) ≈ vcat(0.0, 0.0, p0, 0.0, 0.0, 0.0) + @test map(Base.Fix1(pdf, d), support(d)) ≈ p0 + @test map(Base.Fix1(pdf, d), rmin-2:rmax) ≈ vcat(0.0, 0.0, p0) + @test map(Base.Fix1(pdf, d), rmin:rmax+3) ≈ vcat(p0, 0.0, 0.0, 0.0) + @test map(Base.Fix1(pdf, d), rmin-2:rmax+3) ≈ vcat(0.0, 0.0, p0, 0.0, 0.0, 0.0) elseif islowerbounded(d) - @test pdf.(Ref(d), rmin-2:rmax) ≈ vcat(0.0, 0.0, p0) + @test map(Base.Fix1(pdf, d), rmin-2:rmax) ≈ vcat(0.0, 0.0, p0) end end @@ -444,13 +444,13 @@ function test_evaluation(d::DiscreteUnivariateDistribution, vs::AbstractVector, end # consistency of scalar-based and vectorized evaluation - @test pdf.(Ref(d), vs) ≈ p - @test cdf.(Ref(d), vs) ≈ c - @test ccdf.(Ref(d), vs) ≈ cc + @test Base.Fix1(pdf, d).(vs) ≈ p + @test Base.Fix1(cdf, d).(vs) ≈ c + @test Base.Fix1(ccdf, d).(vs) ≈ cc - @test logpdf.(Ref(d), vs) ≈ lp - @test logcdf.(Ref(d), vs) ≈ lc - @test logccdf.(Ref(d), vs) ≈ lcc + @test Base.Fix1(logpdf, d).(vs) ≈ lp + @test Base.Fix1(logcdf, d).(vs) ≈ lc + @test Base.Fix1(logccdf, d).(vs) ≈ lcc end @@ -511,15 +511,15 @@ function test_evaluation(d::ContinuousUnivariateDistribution, vs::AbstractVector # consistency of scalar-based and vectorized evaluation if !isa(d, StudentizedRange) - @test pdf.(Ref(d), vs) ≈ p - @test logpdf.(Ref(d), vs) ≈ lp + @test Base.Fix1(pdf, d).(vs) ≈ p + @test Base.Fix1(logpdf, d).(vs) ≈ lp end - @test cdf.(Ref(d), vs) ≈ c - @test ccdf.(Ref(d), vs) ≈ cc + @test Base.Fix1(cdf, d).(vs) ≈ c + @test Base.Fix1(ccdf, d).(vs) ≈ cc - @test logcdf.(Ref(d), vs) ≈ lc - @test logccdf.(Ref(d), vs) ≈ lcc + @test Base.Fix1(logcdf, d).(vs) ≈ lc + @test Base.Fix1(logccdf, d).(vs) ≈ lcc end function test_nonfinite(distr::UnivariateDistribution) @@ -550,7 +550,7 @@ function test_stats(d::DiscreteUnivariateDistribution, vs::AbstractVector) # using definition (or an approximation) vf = Float64[v for v in vs] - p = pdf.(Ref(d), vf) + p = Base.Fix1(pdf, d).(vf) xmean = dot(p, vf) xvar = dot(p, abs2.(vf .- xmean)) xstd = sqrt(xvar) diff --git a/test/univariate/continuous/johnsonsu.jl b/test/univariate/continuous/johnsonsu.jl index 5f259167f8..716f1b1df8 100644 --- a/test/univariate/continuous/johnsonsu.jl +++ b/test/univariate/continuous/johnsonsu.jl @@ -10,10 +10,10 @@ @test rand(d1) isa Float64 @test median(d1) == quantile(d1, 0.5) - x = quantile.(d1, [0.25, 0.45, 0.60, 0.80, 0.90]) - @test all(cdf.(d1, x) .≈ [0.25, 0.45, 0.60, 0.80, 0.90]) - y = cquantile.(d1, [0.25, 0.45, 0.60, 0.80, 0.90]) - @test all(ccdf.(d1, y) .≈ [0.25, 0.45, 0.60, 0.80, 0.90]) + x = Base.Fix1(quantile, d1).([0.25, 0.45, 0.60, 0.80, 0.90]) + @test all(Base.Fix1(cdf, d1).(x) .≈ [0.25, 0.45, 0.60, 0.80, 0.90]) + y = Base.Fix1(cquantile, d1).([0.25, 0.45, 0.60, 0.80, 0.90]) + @test all(Base.Fix1(ccdf, d1).(y) .≈ [0.25, 0.45, 0.60, 0.80, 0.90]) @test mean(d1) ≈ 7.581281 @test var(d1) ≈ 19.1969485 diff --git a/test/univariate/continuous/logitnormal.jl b/test/univariate/continuous/logitnormal.jl index 20ccce3bd1..454376606c 100644 --- a/test/univariate/continuous/logitnormal.jl +++ b/test/univariate/continuous/logitnormal.jl @@ -43,7 +43,7 @@ function test_logitnormal(g::LogitNormal, n_tsamples::Int=10^6, for i = 1:min(100, n_tsamples) @test logpdf(g, X[i]) ≈ log(pdf(g, X[i])) end - @test logpdf.(g, X) ≈ log.(pdf.(g, X)) + @test Base.Fix1(logpdf, g).(X) ≈ log.(Base.Fix1(pdf, g).(X)) @test isequal(logpdf(g, 0),-Inf) @test isequal(logpdf(g, 1),-Inf) @test isequal(logpdf(g, -eps()),-Inf) diff --git a/test/univariate/continuous/rician.jl b/test/univariate/continuous/rician.jl index 5674373b36..a75397f891 100644 --- a/test/univariate/continuous/rician.jl +++ b/test/univariate/continuous/rician.jl @@ -14,14 +14,14 @@ @test var(d1) ≈ var(d2) @test mode(d1) ≈ mode(d2) @test median(d1) ≈ median(d2) - @test quantile.(d1, [0.25, 0.45, 0.60, 0.80, 0.90]) ≈ quantile.(d2, [0.25, 0.45, 0.60, 0.80, 0.90]) - @test pdf.(d1, 0.0:0.1:1.0) ≈ pdf.(d2, 0.0:0.1:1.0) - @test cdf.(d1, 0.0:0.1:1.0) ≈ cdf.(d2, 0.0:0.1:1.0) + @test Base.Fix1(quantile, d1).([0.25, 0.45, 0.60, 0.80, 0.90]) ≈ Base.Fix1(quantile, d2).([0.25, 0.45, 0.60, 0.80, 0.90]) + @test Base.Fix1(pdf, d1).(0.0:0.1:1.0) ≈ Base.Fix1(pdf, d2).(0.0:0.1:1.0) + @test Base.Fix1(cdf, d1).(0.0:0.1:1.0) ≈ Base.Fix1(cdf, d2).(0.0:0.1:1.0) d1 = Rician(10.0, 10.0) @test median(d1) == quantile(d1, 0.5) - x = quantile.(d1, [0.25, 0.45, 0.60, 0.80, 0.90]) - @test all(cdf.(d1, x) .≈ [0.25, 0.45, 0.60, 0.80, 0.90]) + x = Base.Fix1(quantile, d1).([0.25, 0.45, 0.60, 0.80, 0.90]) + @test all(Base.Fix1(cdf, d1).(x) .≈ [0.25, 0.45, 0.60, 0.80, 0.90]) x = rand(Rician(5.0, 5.0), 100000) d1 = fit(Rician, x) diff --git a/test/univariate/continuous/skewnormal.jl b/test/univariate/continuous/skewnormal.jl index 81d59ed961..5d1257162c 100644 --- a/test/univariate/continuous/skewnormal.jl +++ b/test/univariate/continuous/skewnormal.jl @@ -25,7 +25,7 @@ import Distributions: normpdf, normcdf, normlogpdf, normlogcdf d4 = Normal(0.5, 2.2) # @test pdf(d3, 3.3) == Distributions.pdf(d4, 3.3) - @test pdf.(d3, 1:3) == Distributions.pdf.(d4, 1:3) + @test Base.Fix1(pdf, d3).(1:3) == Base.Fix1(pdf, d4).(1:3) a = mean(d3), var(d3), std(d3) b = Distributions.mean(d4), Distributions.var(d4), Distributions.std(d4) @test a == b diff --git a/test/univariate/discrete/binomial.jl b/test/univariate/discrete/binomial.jl index 9fd3251f2b..0caa974216 100644 --- a/test/univariate/discrete/binomial.jl +++ b/test/univariate/discrete/binomial.jl @@ -9,14 +9,14 @@ Random.seed!(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) + a = Base.Fix1(pdf, d).(0:n) for t=0:n @test pdf(d, t) ≈ a[1+t] end li = rand(0:n, 2) rng = minimum(li):maximum(li) - b = pdf.(d, rng) + b = Base.Fix1(pdf, d).(rng) for t in rng @test pdf(d, t) ≈ b[t - first(rng) + 1] end diff --git a/test/univariate/discrete/categorical.jl b/test/univariate/discrete/categorical.jl index 45d2c84f7d..a835f7c13c 100644 --- a/test/univariate/discrete/categorical.jl +++ b/test/univariate/discrete/categorical.jl @@ -56,8 +56,8 @@ for p in Any[ @test iszero(ccdf(d, Inf)) @test isnan(ccdf(d, NaN)) - @test pdf.(d, support(d)) == p - @test pdf.(d, 1:k) == p + @test Base.Fix1(pdf, d).(support(d)) == p + @test Base.Fix1(pdf, d).(1:k) == p @test cf(d, 0) ≈ 1.0 @test cf(d, 1) ≈ p' * cis.(1:length(p)) diff --git a/test/univariate/discrete/poissonbinomial.jl b/test/univariate/discrete/poissonbinomial.jl index 8e53e6edd6..611a54c258 100644 --- a/test/univariate/discrete/poissonbinomial.jl +++ b/test/univariate/discrete/poissonbinomial.jl @@ -104,9 +104,9 @@ for (n₁, n₂, n₃, p₁, p₂, p₃) in [(10, 10, 10, 0.1, 0.5, 0.9), b2 = Binomial(n₂, p₂) b3 = Binomial(n₃, p₃) - pmf1 = pdf.(b1, support(b1)) - pmf2 = pdf.(b2, support(b2)) - pmf3 = pdf.(b3, support(b3)) + pmf1 = Base.Fix1(pdf, b1).(support(b1)) + pmf2 = Base.Fix1(pdf, b2).(support(b2)) + pmf3 = Base.Fix1(pdf, b3).(support(b3)) @test @inferred(mean(d)) ≈ (mean(b1) + mean(b2) + mean(b3)) @test @inferred(var(d)) ≈ (var(b1) + var(b2) + var(b3)) diff --git a/test/univariate/discrete/soliton.jl b/test/univariate/discrete/soliton.jl index e0e15e101e..af0aefddc1 100644 --- a/test/univariate/discrete/soliton.jl +++ b/test/univariate/discrete/soliton.jl @@ -5,7 +5,7 @@ using Distributions Ω = Soliton(K, M, δ, atol) @test pdf(Ω, M) > pdf(Ω, M-1) @test pdf(Ω, M) > pdf(Ω, M+1) - @test cumsum(pdf.(Ω, 1:K)) ≈ cdf.(Ω, 1:K) + @test cumsum(Base.Fix1(pdf, Ω).(1:K)) ≈ Base.Fix1(cdf, Ω).(1:K) @test cdf(Ω, 0) ≈ 0 @test cdf(Ω, K) ≈ 1 @test mean(Ω) ≈ 7.448826535558562 @@ -20,7 +20,7 @@ using Distributions K, M, δ, atol = 100, 60, 0.2, 1e-3 Ω = Soliton(K, M, δ, atol) ds = [d for d in 1:K if pdf(Ω, d) > 0] - @test all(pdf.(Ω, ds) .> atol) + @test all(Base.Fix1(pdf, Ω).(ds) .> atol) @test cdf(Ω, 0) ≈ 0 @test cdf(Ω, K) ≈ 1 @test minimum(Ω) == 1 diff --git a/test/univariate/locationscale.jl b/test/univariate/locationscale.jl index 761ce55bbf..c4876e1970 100644 --- a/test/univariate/locationscale.jl +++ b/test/univariate/locationscale.jl @@ -83,9 +83,9 @@ function test_location_scale( insupport(dtest, -x) == insupport(dref, -x) @test pdf(dtest, x) ≈ pdf(dref, x) - @test pdf.(dtest, xs) ≈ pdf.(dref, xs) + @test Base.Fix1(pdf, dtest).(xs) ≈ Base.Fix1(pdf, dref).(xs) @test logpdf(dtest, x) ≈ logpdf(dref, x) - @test logpdf.(dtest, xs) ≈ logpdf.(dref, xs) + @test Base.Fix1(logpdf, dtest).(xs) ≈ Base.Fix1(logpdf, dref).(xs) @test loglikelihood(dtest, x) ≈ loglikelihood(dref, x) @test loglikelihood(dtest, xs) ≈ loglikelihood(dref, xs) diff --git a/test/univariates.jl b/test/univariates.jl index e16c52f59b..669ddfe0e7 100644 --- a/test/univariates.jl +++ b/test/univariates.jl @@ -99,8 +99,8 @@ function verify_and_test(D::Union{Type,Function}, d::UnivariateDistribution, dct # pdf method is not implemented for StudentizedRange if !isa(d, StudentizedRange) - @test isapprox(pdf.(d, x), p; atol=1e-16, rtol=1e-8) - @test isapprox(logpdf.(d, x), lp; atol=isa(d, NoncentralHypergeometric) ? 1e-4 : 1e-12) + @test Base.Fix1(pdf, d).(x) ≈ p atol=1e-16 rtol=1e-8 + @test Base.Fix1(logpdf, d).(x) ≈ lp atol=isa(d, NoncentralHypergeometric) ? 1e-4 : 1e-12 end # cdf method is not implemented for NormalInverseGaussian From de9d5cfd8aa3f9e770faf68d8d689bcc425af60c Mon Sep 17 00:00:00 2001 From: David Widmann Date: Fri, 5 Jan 2024 12:21:42 +0100 Subject: [PATCH 15/26] Update Project.toml --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index 60ae8b2618..0633a6f0b0 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "Distributions" uuid = "31c24e10-a181-5473-b8eb-7969acd0382f" authors = ["JuliaStats"] -version = "0.25.105" +version = "0.25.106" [deps] ChainRulesCore = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4" From 07d0ed520c80e272e1985aa505fc9b3cc20cef84 Mon Sep 17 00:00:00 2001 From: Ruben Seyer Date: Mon, 8 Jan 2024 15:09:20 +0100 Subject: [PATCH 16/26] Fix Gumbel gradlogpdf --- src/univariate/continuous/gumbel.jl | 2 +- test/gradlogpdf.jl | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/src/univariate/continuous/gumbel.jl b/src/univariate/continuous/gumbel.jl index 92898ca7a2..8f062fa8d5 100644 --- a/src/univariate/continuous/gumbel.jl +++ b/src/univariate/continuous/gumbel.jl @@ -102,4 +102,4 @@ logcdf(d::Gumbel, x::Real) = -exp(-zval(d, x)) quantile(d::Gumbel, p::Real) = d.μ - d.θ * log(-log(p)) -gradlogpdf(d::Gumbel, x::Real) = - (1 + exp((d.μ - x) / d.θ)) / d.θ +gradlogpdf(d::Gumbel, x::Real) = (-1 + exp(-zval(d, x))) / d.θ diff --git a/test/gradlogpdf.jl b/test/gradlogpdf.jl index 5f54340ffd..f4216a67e6 100644 --- a/test/gradlogpdf.jl +++ b/test/gradlogpdf.jl @@ -9,7 +9,7 @@ using Test @test isapprox(gradlogpdf(Chisq(7.0), 12.0) , -0.29166666666666663, atol=1.0e-8) @test isapprox(gradlogpdf(Exponential(2.0), 7.0) , -0.5 , atol=1.0e-8) @test isapprox(gradlogpdf(Gamma(9.0, 0.5), 11.0) , -1.2727272727272727 , atol=1.0e-8) -@test isapprox(gradlogpdf(Gumbel(3.5, 1.0), 4.0) , -1.6065306597126334 , atol=1.0e-8) +@test isapprox(gradlogpdf(Gumbel(3.5, 1.0), 4.0) , -0.3934693402873666 , atol=1.0e-8) @test isapprox(gradlogpdf(Laplace(7.0), 34.0) , -1.0 , atol=1.0e-8) @test isapprox(gradlogpdf(Logistic(-6.0), 1.0) , -0.9981778976111987 , atol=1.0e-8) @test isapprox(gradlogpdf(LogNormal(5.5), 2.0) , 1.9034264097200273 , atol=1.0e-8) From 8f68d065590a4ddba3684ef40509fb71d151a727 Mon Sep 17 00:00:00 2001 From: Ruben Seyer Date: Mon, 8 Jan 2024 15:47:27 +0100 Subject: [PATCH 17/26] Use expm1 Co-authored-by: David Widmann --- src/univariate/continuous/gumbel.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/univariate/continuous/gumbel.jl b/src/univariate/continuous/gumbel.jl index 8f062fa8d5..7977a17027 100644 --- a/src/univariate/continuous/gumbel.jl +++ b/src/univariate/continuous/gumbel.jl @@ -102,4 +102,4 @@ logcdf(d::Gumbel, x::Real) = -exp(-zval(d, x)) quantile(d::Gumbel, p::Real) = d.μ - d.θ * log(-log(p)) -gradlogpdf(d::Gumbel, x::Real) = (-1 + exp(-zval(d, x))) / d.θ +gradlogpdf(d::Gumbel, x::Real) = expm1(-zval(d, x)) / d.θ From c1705a3015d438f7e841e82ef5148224813831e8 Mon Sep 17 00:00:00 2001 From: David Widmann Date: Mon, 8 Jan 2024 19:26:57 +0100 Subject: [PATCH 18/26] Update Project.toml --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index 0633a6f0b0..23fbf08012 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "Distributions" uuid = "31c24e10-a181-5473-b8eb-7969acd0382f" authors = ["JuliaStats"] -version = "0.25.106" +version = "0.25.107" [deps] ChainRulesCore = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4" From a59323922d838d836bcf11a2d06930999a081d7c Mon Sep 17 00:00:00 2001 From: "dependabot[bot]" <49699333+dependabot[bot]@users.noreply.github.com> Date: Mon, 22 Jan 2024 16:42:28 +0000 Subject: [PATCH 19/26] Bump actions/cache from 3 to 4 Bumps [actions/cache](https://github.com/actions/cache) from 3 to 4. - [Release notes](https://github.com/actions/cache/releases) - [Changelog](https://github.com/actions/cache/blob/main/RELEASES.md) - [Commits](https://github.com/actions/cache/compare/v3...v4) --- updated-dependencies: - dependency-name: actions/cache dependency-type: direct:production update-type: version-update:semver-major ... Signed-off-by: dependabot[bot] --- .github/workflows/CI.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.github/workflows/CI.yml b/.github/workflows/CI.yml index ec2f9a00c9..63bf812177 100644 --- a/.github/workflows/CI.yml +++ b/.github/workflows/CI.yml @@ -39,7 +39,7 @@ jobs: with: version: ${{ matrix.version }} arch: ${{ matrix.arch }} - - uses: actions/cache@v3 + - uses: actions/cache@v4 env: cache-name: cache-artifacts with: From 6bc66b957375596fde599f361264cbee4e49a517 Mon Sep 17 00:00:00 2001 From: "dependabot[bot]" <49699333+dependabot[bot]@users.noreply.github.com> Date: Mon, 5 Feb 2024 17:00:23 +0000 Subject: [PATCH 20/26] Bump codecov/codecov-action from 3 to 4 Bumps [codecov/codecov-action](https://github.com/codecov/codecov-action) from 3 to 4. - [Release notes](https://github.com/codecov/codecov-action/releases) - [Changelog](https://github.com/codecov/codecov-action/blob/main/CHANGELOG.md) - [Commits](https://github.com/codecov/codecov-action/compare/v3...v4) --- updated-dependencies: - dependency-name: codecov/codecov-action dependency-type: direct:production update-type: version-update:semver-major ... Signed-off-by: dependabot[bot] --- .github/workflows/CI.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.github/workflows/CI.yml b/.github/workflows/CI.yml index ec2f9a00c9..5e7c348ada 100644 --- a/.github/workflows/CI.yml +++ b/.github/workflows/CI.yml @@ -57,7 +57,7 @@ jobs: Pkg.instantiate()' - run: julia --project=perf perf/samplers.jl - uses: julia-actions/julia-processcoverage@v1 - - uses: codecov/codecov-action@v3 + - uses: codecov/codecov-action@v4 with: files: lcov.info docs: From 3d71a2b0431892d739e3b861d2dd1baf9fdc46d3 Mon Sep 17 00:00:00 2001 From: adienes <51664769+adienes@users.noreply.github.com> Date: Mon, 5 Feb 2024 17:18:30 -0500 Subject: [PATCH 21/26] microoptimization on rand(::AliasTable) (#1831) * microoptimization on rand(::AliasTable) * Update src/samplers/aliastable.jl Co-authored-by: David Widmann --------- Co-authored-by: David Widmann --- src/samplers/aliastable.jl | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/src/samplers/aliastable.jl b/src/samplers/aliastable.jl index ddebbce57a..56042930d7 100644 --- a/src/samplers/aliastable.jl +++ b/src/samplers/aliastable.jl @@ -15,9 +15,8 @@ end function rand(rng::AbstractRNG, s::AliasTable) i = rand(rng, 1:length(s.alias)) % Int - u = rand(rng) - @inbounds r = u < s.accept[i] ? i : s.alias[i] - r + # using `ifelse` improves performance here: github.com/JuliaStats/Distributions.jl/pull/1831/ + ifelse(rand(rng) < s.accept[i], i, s.alias[i]) end show(io::IO, s::AliasTable) = @printf(io, "AliasTable with %d entries", ncategories(s)) From 5ca92eaf002c060d32d763f697c1b1ab33f8725d Mon Sep 17 00:00:00 2001 From: Andreas Noack Date: Tue, 6 Feb 2024 11:43:27 +0100 Subject: [PATCH 22/26] Update .github/workflows/CI.yml Co-authored-by: David Widmann --- .github/workflows/CI.yml | 2 ++ 1 file changed, 2 insertions(+) diff --git a/.github/workflows/CI.yml b/.github/workflows/CI.yml index 5e7c348ada..eea6ae8db6 100644 --- a/.github/workflows/CI.yml +++ b/.github/workflows/CI.yml @@ -59,6 +59,8 @@ jobs: - uses: julia-actions/julia-processcoverage@v1 - uses: codecov/codecov-action@v4 with: + token: ${{ secrets.CODECOV_TOKEN }} # required + fail_ci_if_error: true files: lcov.info docs: name: Documentation From 1376ee35aa54ddc04066c67ecebf1a56c422849e Mon Sep 17 00:00:00 2001 From: Andreas Noack Date: Thu, 22 Feb 2024 15:05:09 +0100 Subject: [PATCH 23/26] Remove TruncatedNormal (#1837) * Remove TruncatedNormal It's been deprecated for five years. * Merge TruncatedNormal docstring into Truncated docstring --- src/Distributions.jl | 2 +- src/truncate.jl | 3 +++ src/truncated/normal.jl | 12 ------------ test/ref/continuous_test.lst | 10 +++++----- test/ref/continuous_test.ref.json | 10 +++++----- test/truncate.jl | 2 +- test/univariates.jl | 2 +- 7 files changed, 16 insertions(+), 25 deletions(-) diff --git a/src/Distributions.jl b/src/Distributions.jl index ceb6063c70..cf3ec32886 100644 --- a/src/Distributions.jl +++ b/src/Distributions.jl @@ -365,7 +365,7 @@ Supported distributions: NoncentralF, NoncentralHypergeometric, NoncentralT, Normal, NormalCanon, NormalInverseGaussian, Pareto, PGeneralizedGaussian, Poisson, PoissonBinomial, QQPair, Rayleigh, Rician, Skellam, Soliton, StudentizedRange, SymTriangularDist, TDist, TriangularDist, - Triweight, Truncated, TruncatedNormal, Uniform, UnivariateGMM, + Triweight, Truncated, Uniform, UnivariateGMM, VonMises, VonMisesFisher, WalleniusNoncentralHypergeometric, Weibull, Wishart, ZeroMeanIsoNormal, ZeroMeanIsoNormalCanon, ZeroMeanDiagNormal, ZeroMeanDiagNormalCanon, ZeroMeanFullNormal, diff --git a/src/truncate.jl b/src/truncate.jl index a7b036fc82..45709f6b53 100644 --- a/src/truncate.jl +++ b/src/truncate.jl @@ -82,6 +82,9 @@ end Truncated Generic wrapper for a truncated distribution. + +The *truncated normal distribution* is a particularly important one in the family of truncated distributions. +Unlike the general case, truncated normal distributions support `mean`, `mode`, `modes`, `var`, `std`, and `entropy`. """ struct Truncated{D<:UnivariateDistribution, S<:ValueSupport, T<: Real, TL<:Union{T,Nothing}, TU<:Union{T,Nothing}} <: UnivariateDistribution{S} untruncated::D # the original distribution (untruncated) diff --git a/src/truncated/normal.jl b/src/truncated/normal.jl index a3ff33e1e1..6fb3342731 100644 --- a/src/truncated/normal.jl +++ b/src/truncated/normal.jl @@ -1,15 +1,3 @@ -# Truncated normal distribution -""" - TruncatedNormal(mu, sigma, l, u) - -The *truncated normal distribution* is a particularly important one in the family of truncated distributions. -We provide additional support for this type with `TruncatedNormal` which calls `Truncated(Normal(mu, sigma), l, u)`. -Unlike the general case, truncated normal distributions support `mean`, `mode`, `modes`, `var`, `std`, and `entropy`. -""" -TruncatedNormal - -@deprecate TruncatedNormal(mu::Real, sigma::Real, a::Real, b::Real) truncated(Normal(mu, sigma), a, b) - ### statistics function mode(d::Truncated{<:Normal{<:Real},Continuous,T}) where {T<:Real} diff --git a/test/ref/continuous_test.lst b/test/ref/continuous_test.lst index 0ada1b1834..d403be3e65 100644 --- a/test/ref/continuous_test.lst +++ b/test/ref/continuous_test.lst @@ -195,11 +195,11 @@ TriangularDist(-4, 14, 3) TriangularDist(2, 2000, 500) TriangularDist(1, 3, 2) -TruncatedNormal(0, 1, -2, 2) -TruncatedNormal(3, 10, 7, 8) -TruncatedNormal(27, 3, 0, Inf) -TruncatedNormal(-5, 1, -Inf, -10) -TruncatedNormal(1.8, 1.2, -Inf, 0) +truncated(Normal(0, 1), lower = -2, upper = 2) +truncated(Normal(3, 10), lower = 7, upper = 8) +truncated(Normal(27, 3), lower = 0) +truncated(Normal(-5, 1), upper = -10) +truncated(Normal(1.8, 1.2), upper = 0) Uniform() Uniform(0.0, 2.0) diff --git a/test/ref/continuous_test.ref.json b/test/ref/continuous_test.ref.json index bf3ceca213..ae86f2e8a7 100644 --- a/test/ref/continuous_test.ref.json +++ b/test/ref/continuous_test.ref.json @@ -5157,7 +5157,7 @@ }, { "expr": "truncated(Normal(0, 1), -2, 2)", - "dtype": "TruncatedNormal", + "dtype": "Truncated", "minimum": -2, "maximum": 2, "properties": { @@ -5187,7 +5187,7 @@ }, { "expr": "truncated(Normal(3, 10), 7, 8)", - "dtype": "TruncatedNormal", + "dtype": "Truncated", "minimum": 7, "maximum": 8, "properties": { @@ -5217,7 +5217,7 @@ }, { "expr": "truncated(Normal(27, 3); lower=0)", - "dtype": "TruncatedNormal", + "dtype": "Truncated", "minimum": 0, "maximum": "inf", "properties": { @@ -5247,7 +5247,7 @@ }, { "expr": "truncated(Normal(-5, 1); upper=-10)", - "dtype": "TruncatedNormal", + "dtype": "Truncated", "minimum": "-inf", "maximum": -10, "properties": { @@ -5277,7 +5277,7 @@ }, { "expr": "truncated(Normal(1.8, 1.2); upper=0)", - "dtype": "TruncatedNormal", + "dtype": "Truncated", "minimum": "-inf", "maximum": 0, "properties": { diff --git a/test/truncate.jl b/test/truncate.jl index 4481cc48ed..b409cba93e 100644 --- a/test/truncate.jl +++ b/test/truncate.jl @@ -39,7 +39,7 @@ function verify_and_test_drive(jsonfile, selected, n_tsamples::Int,lower::Int,up println(" testing truncated($(ex),$lower,$upper)") d = truncated(eval(Meta.parse(ex)),lower,upper) - if dtype != Uniform && dtype != DiscreteUniform && dtype != TruncatedNormal # Uniform is truncated to Uniform + if dtype != Uniform && dtype != DiscreteUniform # Uniform is truncated to Uniform @assert isa(dtype, Type) && dtype <: UnivariateDistribution @test isa(d, dtypet) # verification and testing diff --git a/test/univariates.jl b/test/univariates.jl index 669ddfe0e7..66fb2df368 100644 --- a/test/univariates.jl +++ b/test/univariates.jl @@ -54,7 +54,7 @@ function verify_and_test(D::Union{Type,Function}, d::UnivariateDistribution, dct # Note: properties include all applicable params and stats # - # D can be a function, e.g. TruncatedNormal + # D can be a function if isa(D, Type) @assert isa(d, D) end From f33af97da34a2e772837b63891c0ed52a3a1a989 Mon Sep 17 00:00:00 2001 From: "dependabot[bot]" <49699333+dependabot[bot]@users.noreply.github.com> Date: Mon, 1 Apr 2024 23:28:00 +0200 Subject: [PATCH 24/26] Bump julia-actions/setup-julia from 1 to 2 (#1846) Bumps [julia-actions/setup-julia](https://github.com/julia-actions/setup-julia) from 1 to 2. - [Release notes](https://github.com/julia-actions/setup-julia/releases) - [Commits](https://github.com/julia-actions/setup-julia/compare/v1...v2) --- updated-dependencies: - dependency-name: julia-actions/setup-julia dependency-type: direct:production update-type: version-update:semver-major ... Signed-off-by: dependabot[bot] Co-authored-by: dependabot[bot] <49699333+dependabot[bot]@users.noreply.github.com> --- .github/workflows/CI.yml | 4 ++-- .github/workflows/IntegrationTest.yml | 2 +- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/.github/workflows/CI.yml b/.github/workflows/CI.yml index fcd3150f9b..ccf2aea36f 100644 --- a/.github/workflows/CI.yml +++ b/.github/workflows/CI.yml @@ -35,7 +35,7 @@ jobs: - x64 steps: - uses: actions/checkout@v4 - - uses: julia-actions/setup-julia@v1 + - uses: julia-actions/setup-julia@v2 with: version: ${{ matrix.version }} arch: ${{ matrix.arch }} @@ -67,7 +67,7 @@ jobs: runs-on: ubuntu-latest steps: - uses: actions/checkout@v4 - - uses: julia-actions/setup-julia@v1 + - uses: julia-actions/setup-julia@v2 with: version: '1' - run: | diff --git a/.github/workflows/IntegrationTest.yml b/.github/workflows/IntegrationTest.yml index 3aec5d9002..38079be1c7 100644 --- a/.github/workflows/IntegrationTest.yml +++ b/.github/workflows/IntegrationTest.yml @@ -32,7 +32,7 @@ jobs: steps: - uses: actions/checkout@v4 - - uses: julia-actions/setup-julia@v1 + - uses: julia-actions/setup-julia@v2 with: version: 1 arch: x64 From b670fee8ae0dff56bcf4af7ddfdaa7e36d67caa1 Mon Sep 17 00:00:00 2001 From: Lilith Orion Hafner Date: Sat, 20 Apr 2024 02:07:17 -0500 Subject: [PATCH 25/26] Use a faster implementation of AliasTables (#1848) * switch to AliasTables.jl * retune heuristic * add test for #832 * add more tests * move alias table import and tighten from using to import * Back out multinomial heuristic adjustment at @adienes's request * Update test/univariate/discrete/categorical.jl (style) Co-authored-by: David Widmann --------- Co-authored-by: David Widmann --- Project.toml | 2 ++ src/Distributions.jl | 2 ++ src/samplers/aliastable.jl | 23 ++++------------------- test/samplers.jl | 2 ++ test/univariate/discrete/categorical.jl | 17 +++++++++++++++-- 5 files changed, 25 insertions(+), 21 deletions(-) diff --git a/Project.toml b/Project.toml index 23fbf08012..1072d0294e 100644 --- a/Project.toml +++ b/Project.toml @@ -4,6 +4,7 @@ authors = ["JuliaStats"] version = "0.25.107" [deps] +AliasTables = "66dad0bd-aa9a-41b7-9441-69ab47430ed8" ChainRulesCore = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4" DensityInterface = "b429d917-457f-4dbc-8f4c-0cc954292b1d" FillArrays = "1a297f60-69ca-5386-bcde-b61e274b549b" @@ -30,6 +31,7 @@ DistributionsDensityInterfaceExt = "DensityInterface" DistributionsTestExt = "Test" [compat] +AliasTables = "1" Aqua = "0.8" Calculus = "0.5" ChainRulesCore = "1" diff --git a/src/Distributions.jl b/src/Distributions.jl index cf3ec32886..8d344c4158 100644 --- a/src/Distributions.jl +++ b/src/Distributions.jl @@ -27,6 +27,8 @@ import PDMats: dim, PDMat, invquad using SpecialFunctions using Base.MathConstants: eulergamma +import AliasTables + export # re-export Statistics mean, median, quantile, std, var, cov, cor, diff --git a/src/samplers/aliastable.jl b/src/samplers/aliastable.jl index 56042930d7..1f633cf716 100644 --- a/src/samplers/aliastable.jl +++ b/src/samplers/aliastable.jl @@ -1,22 +1,7 @@ struct AliasTable <: Sampleable{Univariate,Discrete} - accept::Vector{Float64} - alias::Vector{Int} + at::AliasTables.AliasTable{UInt64, Int} + AliasTable(probs::AbstractVector{<:Real}) = new(AliasTables.AliasTable(probs)) end -ncategories(s::AliasTable) = length(s.alias) - -function AliasTable(probs::AbstractVector) - n = length(probs) - n > 0 || throw(ArgumentError("The input probability vector is empty.")) - accp = Vector{Float64}(undef, n) - alias = Vector{Int}(undef, n) - StatsBase.make_alias_table!(probs, 1.0, accp, alias) - AliasTable(accp, alias) -end - -function rand(rng::AbstractRNG, s::AliasTable) - i = rand(rng, 1:length(s.alias)) % Int - # using `ifelse` improves performance here: github.com/JuliaStats/Distributions.jl/pull/1831/ - ifelse(rand(rng) < s.accept[i], i, s.alias[i]) -end - +ncategories(s::AliasTable) = length(s.at) +rand(rng::AbstractRNG, s::AliasTable) = rand(rng, s.at) show(io::IO, s::AliasTable) = @printf(io, "AliasTable with %d entries", ncategories(s)) diff --git a/test/samplers.jl b/test/samplers.jl index 2744ae9acb..749b45d0d6 100644 --- a/test/samplers.jl +++ b/test/samplers.jl @@ -31,9 +31,11 @@ import Distributions: @testset "p=$p" for p in Any[[1.0], [0.3, 0.7], [0.2, 0.3, 0.4, 0.1]] test_samples(S(p), Categorical(p), n_tsamples) test_samples(S(p), Categorical(p), n_tsamples, rng=rng) + @test ncategories(S(p)) == length(p) end end + @test string(AliasTable(Float16[1,2,3])) == "AliasTable with 3 entries" ## Binomial samplers diff --git a/test/univariate/discrete/categorical.jl b/test/univariate/discrete/categorical.jl index a835f7c13c..6d87d4dc86 100644 --- a/test/univariate/discrete/categorical.jl +++ b/test/univariate/discrete/categorical.jl @@ -93,9 +93,9 @@ end end @testset "reproducibility across julia versions" begin - d= Categorical([0.1, 0.2, 0.7]) + d = Categorical([0.1, 0.2, 0.7]) rng = StableRNGs.StableRNG(600) - @test rand(rng, d, 10) == [2, 1, 3, 3, 2, 3, 3, 3, 3, 3] + @test rand(rng, d, 10) == [3, 1, 1, 2, 3, 2, 3, 3, 2, 3] end @testset "comparisons" begin @@ -124,4 +124,17 @@ end @test Categorical([0.5, 0.5]) ≈ Categorical([0.5f0, 0.5f0]) end +@testset "issue #832" begin + priorities = collect(Float64, 1:1000) + priorities[1:50] .= 1e8 + + at = Distributions.AliasTable(priorities) + iat = rand(at, 16) + + # failure rate of a single sample is sum(51:1000)/50e8 = 9.9845e-5 + # failure rate of 4 out of 16 samples is 1-cdf(Binomial(16, 9.9845e-5), 3) = 1.8074430840897548e-13 + # this test should randomly fail with a probability of 1.8074430840897548e-13 + @test count(==(1e8), priorities[iat]) >= 13 +end + end From 818814fd4ebd46b53e37a6768d33877b546bf8fa Mon Sep 17 00:00:00 2001 From: David Widmann Date: Sat, 20 Apr 2024 09:08:00 +0200 Subject: [PATCH 26/26] Update Project.toml --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index 1072d0294e..39def97c09 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "Distributions" uuid = "31c24e10-a181-5473-b8eb-7969acd0382f" authors = ["JuliaStats"] -version = "0.25.107" +version = "0.25.108" [deps] AliasTables = "66dad0bd-aa9a-41b7-9441-69ab47430ed8"