Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Add MvLogitNormal #1774

Merged
merged 21 commits into from
Sep 28, 2023
Merged

Add MvLogitNormal #1774

merged 21 commits into from
Sep 28, 2023

Conversation

sethaxen
Copy link
Contributor

@sethaxen sethaxen commented Sep 8, 2023

Implements #1771

@sethaxen
Copy link
Contributor Author

sethaxen commented Sep 8, 2023

Currently it takes a different approach than MvLogNormal. MvLogNormal wraps an MvNormal, so it cannot be parameterized like an MvNormalCanon or other AbstractMvNormal types. MvLogitNormal instead wraps an AbstractMvNormal, which allows it to also wrap MvNormalCanon or any user-defined AbstractMvNormal type.

@codecov-commenter
Copy link

codecov-commenter commented Sep 8, 2023

Codecov Report

Attention: 2 lines in your changes are missing coverage. Please review.

Files Coverage Δ
src/Distributions.jl 100.00% <ø> (ø)
src/multivariates.jl 44.82% <ø> (ø)
src/multivariate/mvlogitnormal.jl 97.14% <97.14%> (ø)

... and 1 file with indirect coverage changes

📢 Thoughts on this report? Let us know!.

@sethaxen sethaxen marked this pull request as ready for review September 10, 2023 08:48
MvLogitNormal(d::AbstractMvNormal) = MvLogitNormal{typeof(d)}(d)
MvLogitNormal(args...) = MvLogitNormal(MvNormal(args...))

function Base.show(io::IO, d::MvLogitNormal)
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Is there a reason the generic methods don't work? I think we should probably be using something more generic.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The generic methods treat the field names as parameters, which here looks weird because the field contains another distribution:

julia> using Distributions, LinearAlgebra

julia> d = MvLogitNormal(randn(3), exp(Symmetric(randn(3, 3))))
MvLogitNormal{FullNormal}(
dim: 3
μ: [-0.6695883087710813, -0.4262204488438936, 1.5508073295811842]
Σ: [1.8911987804702053 0.038814771826912765 1.2907448948413436; 0.038814771826912765 0.7054789545738211 0.06748072828787553; 1.2907448948413436 0.06748072828787553 0.9560938804684098]
)


julia> invoke(Base.show, Tuple{IO,Distribution}, stdout, d)
MvLogitNormal{FullNormal}(
normal: FullNormal(
dim: 3
μ: [-0.6695883087710813, -0.4262204488438936, 1.5508073295811842]
Σ: [1.8911987804702053 0.038814771826912765 1.2907448948413436; 0.038814771826912765 0.7054789545738211 0.06748072828787553; 1.2907448948413436 0.06748072828787553 0.9560938804684098]
)

)

julia> d = MvLogitNormal{MvNormalCanon}(randn(3), exp(Symmetric(randn(3, 3))))
MvLogitNormal{FullNormalCanon}(
μ: [-12.008515816447346, 9.221467028160104, 0.9303379007669893]
h: [-1.3761114231540177, 2.320365757512754, -0.05647292940406612]
J: [0.7714588053672484 0.9492070250024665 -0.929895944697952; 0.9492070250024665 1.6622992471750913 -1.7304512340037166; -0.929895944697952 -1.7304512340037166 5.088641347774431]
)


julia> invoke(Base.show, Tuple{IO,Distribution}, stdout, d)
MvLogitNormal{FullNormalCanon}(
normal: FullNormalCanon(
μ: [-12.008515816447346, 9.221467028160104, 0.9303379007669893]
h: [-1.3761114231540177, 2.320365757512754, -0.05647292940406612]
J: [0.7714588053672484 0.9492070250024665 -0.929895944697952; 0.9492070250024665 1.6622992471750913 -1.7304512340037166; -0.929895944697952 -1.7304512340037166 5.088641347774431]
)

)

We could use show_multiline, but this requires knowledge of the names of the parameters of the AbstractMvNormal, which we don't know, since Distributions doesn't have an API function for parameter names, just values.

src/multivariate/mvlogitnormal.jl Outdated Show resolved Hide resolved
src/multivariate/mvlogitnormal.jl Outdated Show resolved Hide resolved
@sethaxen sethaxen requested a review from devmotion September 28, 2023 07:50
Copy link
Member

@devmotion devmotion left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Looks very good, thank you! I only had a few minor comments.

src/multivariate/mvlogitnormal.jl Outdated Show resolved Hide resolved
src/multivariate/mvlogitnormal.jl Show resolved Hide resolved
src/multivariate/mvlogitnormal.jl Show resolved Hide resolved
test/multivariate/mvlogitnormal.jl Show resolved Hide resolved
test/multivariate/mvlogitnormal.jl Outdated Show resolved Hide resolved
test/multivariate/mvlogitnormal.jl Outdated Show resolved Hide resolved
Ymean = vec(mean(Y; dims=2))
Ycov = cov(Y; dims=2)
for i in 1:(length(d) - 1)
@test isapprox(
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I'm not completely sure anymore but I assume we don't have any statistical test utilities for e.g. hypothesis tests of samples? IIRC KS tests are used for the univariate distributions but I assume we don't have anything for multi- or matrixvariate distributions.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This is the closest I can find, which just checks that a 1x1 matrix distribution is consistent with the univariate distribution it reduces to:

function test_draws_against_univariate_cdf(D::MatrixDistribution, d::UnivariateDistribution)
α = 0.025
M = 100000
matvardraws = [rand(D)[1] for m in 1:M]
@test pvalue_kolmogorovsmirnoff(matvardraws, d) >= α
nothing
end

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Thank you for checking! Would be useful to have something similar for multivariate distributions (generally, IMO we should provide test utilities such that you can check more easily that you implement the interface correctly) but better to do this in a separate PR.

Co-authored-by: David Widmann <devmotion@users.noreply.github.com>
@devmotion devmotion merged commit cd5d4cc into JuliaStats:master Sep 28, 2023
@sethaxen sethaxen deleted the mvlogitnormal branch September 28, 2023 11:48
@sethaxen sethaxen mentioned this pull request Sep 29, 2023
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging this pull request may close these issues.

4 participants