Extend pem to multivariate independent distributions#6
Draft
Conversation
…variate results Agent-Logs-Url: https://github.com/davide-f/PointEstimateMethod.jl/sessions/7834bdd4-da76-4ccf-b55b-66a2a17c96e2 Co-authored-by: davide-f <67809479+davide-f@users.noreply.github.com>
…ult value Agent-Logs-Url: https://github.com/davide-f/PointEstimateMethod.jl/sessions/7834bdd4-da76-4ccf-b55b-66a2a17c96e2 Co-authored-by: davide-f <67809479+davide-f@users.noreply.github.com>
Copilot
AI
changed the title
[WIP] Extend PointEstimateModel to support multivariate distributions
Extend pem to multivariate independent distributions
Apr 20, 2026
davide-f
reviewed
Apr 20, 2026
Owner
There was a problem hiding this comment.
@copilot No: the result shall not be the product of the combos.
Let K be the number of univariate distributions and N the number of concentrations, then the total number of points is K*N.
In particular, each concentration (or point) is in the form of:
(\mu_1, \mu_2, ..., x_{k,n}, ..., \mu_{N-1}, \mu_N) with n=1, ...N and k=1,...,K
There was a problem hiding this comment.
Pull request overview
This PR extends pem beyond univariate usage by adding a multivariate overload for K independent univariate distributions, producing a Cartesian-product grid of point estimates and joint probabilities. It also updates the test suite to validate multivariate outputs via YAML fixtures.
Changes:
- Added
pem(::Vector{<:UnivariateDistribution}, ::Union{Integer, Vector{<:Integer}})to buildx::K×Mandp::Mvia Cartesian products. - Narrowed the experimental-sample overload from
pem(d::Vector, N::Integer)topem(d::Vector{<:Real}, N::Integer)to avoid dispatch ambiguity. - Added multivariate examples/tests plus golden YAML fixtures under
test/testcases/pem_multivariate/.
Reviewed changes
Copilot reviewed 7 out of 7 changed files in this pull request and generated 3 comments.
Show a summary per file
| File | Description |
|---|---|
src/pem.jl |
Adds multivariate pem overload; narrows vector-sample overload signature. |
test/runtests.jl |
Adds test_multivariate_example and runs multivariate examples. |
test/Examples.jl |
Introduces MultivariateExample and a list of multivariate test cases. |
test/testcases/pem_multivariate/Normal2_N2.yml |
Golden fixture for 2D Normal with N=2. |
test/testcases/pem_multivariate/Normal2_N3.yml |
Golden fixture for 2D Normal with N=3. |
test/testcases/pem_multivariate/Normal_TruncNormal_N2.yml |
Golden fixture for Normal × TruncNormal with N=2. |
test/testcases/pem_multivariate/Normal_TruncNormal_mixed_N.yml |
Golden fixture for Normal × TruncNormal with mixed N=[2,3]. |
💡 Add Copilot custom instructions for smarter, more guided reviews. Learn how to get started.
| path_solution = joinpath(BASE_FOLDER, "test", "testcases", string(testing_function) * "_multivariate", example_name * ".yml") | ||
|
|
||
| # Sort columns lexicographically for deterministic comparison | ||
| I_calc = sortperm(eachcol(calc_solution.x), by=col -> collect(col)) |
There was a problem hiding this comment.
sortperm(eachcol(calc_solution.x), by=col -> collect(col)) will attempt to order Vectors using isless, which is not defined for arrays in Julia and typically throws a MethodError. Convert columns to an orderable key (e.g., Tuple(col) or NTuple) or provide an explicit lt that compares tuples lexicographically so the multivariate tests are deterministic and runnable.
Suggested change
| I_calc = sortperm(eachcol(calc_solution.x), by=col -> collect(col)) | |
| I_calc = sortperm(eachcol(calc_solution.x), by=col -> Tuple(col)) |
Comment on lines
+273
to
+280
| # Form the Cartesian product of all point index combinations | ||
| combos = vec(collect(Iterators.product([1:Ns[k] for k in 1:K]...))) | ||
|
|
||
| M = length(combos) | ||
| x = Matrix{Float64}(undef, K, M) | ||
| p = Vector{Float64}(undef, M) | ||
|
|
||
| for (j, combo) in enumerate(combos) |
There was a problem hiding this comment.
The multivariate implementation materializes the full Cartesian product twice (collect(Iterators.product(...)) into combos, then iterating it to fill x/p). This adds avoidable memory pressure proportional to M (which can be huge). Iterate the Iterators.product(...) directly inside the enumerate loop and drop combos to reduce allocations without changing behavior.
Suggested change
| # Form the Cartesian product of all point index combinations | |
| combos = vec(collect(Iterators.product([1:Ns[k] for k in 1:K]...))) | |
| M = length(combos) | |
| x = Matrix{Float64}(undef, K, M) | |
| p = Vector{Float64}(undef, M) | |
| for (j, combo) in enumerate(combos) | |
| # Form the Cartesian product of all point index combinations lazily | |
| M = prod(Ns) | |
| x = Matrix{Float64}(undef, K, M) | |
| p = Vector{Float64}(undef, M) | |
| for (j, combo) in enumerate(Iterators.product([1:Ns[k] for k in 1:K]...)) |
Comment on lines
90
to
92
| function pem( | ||
| d::Vector, | ||
| d::Vector{<:Real}, | ||
| N::Integer; |
There was a problem hiding this comment.
The signature was narrowed to d::Vector{<:Real}, but the docstring above still documents the parameter as d :: Vector (generic). Please update the docstring to reflect the new accepted type (and consider AbstractVector{<:Real} if you want to keep supporting views/subarrays while still avoiding the dispatch ambiguity).
Owner
|
@copilot No: the result shall not be the product of the combos. In particular, each concentration (or point) is in the form of: |
…ead of Cartesian product Agent-Logs-Url: https://github.com/davide-f/PointEstimateMethod.jl/sessions/a0b94002-41ae-4fce-b2b3-3727252b4461 Co-authored-by: davide-f <67809479+davide-f@users.noreply.github.com>
Contributor
Author
... Fixed in commit 8698bbd. The implementation now uses the Hong (1998) concentration point scheme: for K distributions with N concentrations each, there are K×N total points where each point
Fixed in commit 8698bbd. The implementation now uses the Hong (1998) concentration point scheme: for K distributions with N concentrations each, there are K×N total points where each point Warning Firewall rules blocked me from connecting to one or more addresses (expand for details)I tried to connect to the following addresses, but was blocked by firewall rules:
If you need me to access, download, or install something from one of these locations, you can either:
|
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
pemonly supported single univariate distributions. This extends it to handle K independent univariate distributions using the Hong (1998) concentration point scheme.Changes
pemoverload (src/pem.jl): acceptsVector{<:UnivariateDistribution}andN::Union{Integer, Vector{<:Integer}}(scalar applies same N to all; vector allows per-distribution N). Returns(x::Matrix{Float64}, p::Vector{Float64})wherexis K×M with M = K×N (each column is one concentration point) andp[j]is the weightp_{k,n} / K.(μ_1, ..., x_{k,n}, ..., μ_K)keeps all variables at their means except variable k at concentration n. Weight isp_{k,n} / K, ensuring all weights sum to 1.pem(d::Vector, N::Integer)topem(d::Vector{<:Real}, N::Integer)to avoid dispatch collision.MultivariateExamplestruct, 4 new test cases,test_multivariate_examplehelper with lexicographic column sorting (Tuple(col)), and YAML fixtures.Usage
Implements the multivariate extension from Hong (1998).