From 0103accd30dd4b2b1bdbe84870eb58b6e9db1f4c Mon Sep 17 00:00:00 2001 From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com> Date: Mon, 20 Apr 2026 22:16:17 +0000 Subject: [PATCH 1/4] Initial plan From 687c44f30df0f6d0981374b09fcc233e540529f7 Mon Sep 17 00:00:00 2001 From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com> Date: Mon, 20 Apr 2026 22:26:41 +0000 Subject: [PATCH 2/4] Extend pem to multivariate distributions via Cartesian product of univariate 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> --- src/pem.jl | 84 ++++++++++++++++++- test/Examples.jl | 21 +++++ test/runtests.jl | 47 +++++++++++ .../testcases/pem_multivariate/Normal2_N2.yml | 16 ++++ .../testcases/pem_multivariate/Normal2_N3.yml | 31 +++++++ .../Normal_TruncNormal_N2.yml | 16 ++++ .../Normal_TruncNormal_mixed_N.yml | 22 +++++ 7 files changed, 236 insertions(+), 1 deletion(-) create mode 100644 test/testcases/pem_multivariate/Normal2_N2.yml create mode 100644 test/testcases/pem_multivariate/Normal2_N3.yml create mode 100644 test/testcases/pem_multivariate/Normal_TruncNormal_N2.yml create mode 100644 test/testcases/pem_multivariate/Normal_TruncNormal_mixed_N.yml diff --git a/src/pem.jl b/src/pem.jl index b1b8762..32ef83c 100644 --- a/src/pem.jl +++ b/src/pem.jl @@ -88,7 +88,7 @@ Returns """ function pem( - d::Vector, + d::Vector{<:Real}, N::Integer; mean_fun::Function=Distributions.mean, central_moment_fun::Function=Distributions.moment, @@ -202,5 +202,87 @@ function pem( x = ϵ .+ mean_value # results = NamedTuple{(:x, :p)}.(zip(x, p)) + return (x=x, p=p) +end + + +""" + pem(distributions, N; mean_fun=mean, central_moment_fun=moment, optimizer=HiGHS.Optimizer) + +Point Estimate Method to identify estimate points for K independent univariate distributions. + +For K independent distributions, the multivariate representation is constructed by taking +the Cartesian product of the individual point estimates. Each combination has probability +equal to the product of the individual point probabilities. + +This function implements the multivariate extension described in: +- H.P. Hong, An efficient point estimate method for probabilistic analysis, + Reliability Engineering and System Safety, 1998, https://doi.org/10.1016/S0951-8320(97)00071-9 + +Parameters +---------- +- distributions :: Vector{<:UnivariateDistribution} + Vector of K independent univariate distributions +- N :: Union{Integer, Vector{<:Integer}} + Number of desired estimate points per distribution. If an Integer is provided, + the same N is used for all distributions. If a Vector is provided, N[k] is + used for the k-th distribution. +- mean_fun :: Function (optional) + Function used to calculate the mean value of each distribution +- central_moment_fun :: Function (optional) + Function used to calculate the central moment of each distribution +- montecarlo_sampling :: Integer (optional, default 1e6) + Number of Monte Carlo samples used if a non-specific moment function is available +- optimizer (optional) + JuMP optimizer for executing the optimization + +Returns +------- +- (x, p) :: NamedTuple + - x :: Matrix{Float64} + Matrix of shape (K, M) where M = prod(N_1, ..., N_K). Each column is one + combination of point estimates across the K distributions. + - p :: Vector{Float64} + Joint probability of each combination point (product of individual probabilities). + +""" +function pem( + distributions::Vector{<:UnivariateDistribution}, + N::Union{Integer, Vector{<:Integer}}; + mean_fun::Function=Distributions.mean, + central_moment_fun::Function=Distributions.moment, + montecarlo_sampling::Integer=1000000, + optimizer=DEFAULT_SOLVER, + ) + + K = length(distributions) + Ns = N isa Integer ? fill(N, K) : N + + @assert length(Ns) == K "Length of N vector must match the number of distributions (got $(length(Ns)), expected $K)" + + # Apply univariate pem to each distribution independently + results = [ + pem(distributions[k], Ns[k]; + mean_fun=mean_fun, + central_moment_fun=central_moment_fun, + montecarlo_sampling=montecarlo_sampling, + optimizer=optimizer) + for k in 1:K + ] + + # 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) + for k in 1:K + x[k, j] = real(results[k].x[combo[k]]) + end + p[j] = prod(real(results[k].p[combo[k]]) for k in 1:K) + end + return (x=x, p=p) end \ No newline at end of file diff --git a/test/Examples.jl b/test/Examples.jl index 07a2839..e14dec9 100644 --- a/test/Examples.jl +++ b/test/Examples.jl @@ -15,4 +15,25 @@ module Examples Example("Normal_truncated_9", truncated(Normal(1.0, 0.4), 0.0, +Inf), 9), ] + struct MultivariateExample + name # name of test + d # Vector of distributions + N # Number of estimate points (Integer or Vector) + end + + list_multivariate_examples = [ + MultivariateExample("Normal2_N2", [Normal(), Normal()], 2), + MultivariateExample("Normal2_N3", [Normal(), Normal()], 3), + MultivariateExample( + "Normal_TruncNormal_N2", + [Normal(), truncated(Normal(1.0, 0.4), 0.0, +Inf)], + 2, + ), + MultivariateExample( + "Normal_TruncNormal_mixed_N", + [Normal(), truncated(Normal(1.0, 0.4), 0.0, +Inf)], + [2, 3], + ), + ] + end \ No newline at end of file diff --git a/test/runtests.jl b/test/runtests.jl index b56922c..576a0e1 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -50,12 +50,59 @@ function test_example(example_name, testing_function, args...) end +"Function to test multivariate examples" +function test_multivariate_example(example_name, testing_function, args...) + + # for reproducibility + Distributions.Random.seed!(SEED) + + # calculate simulations + calc_solution = testing_function(args...) + + 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)) + + x_sorted = calc_solution.x[:, I_calc] + p_sorted = calc_solution.p[I_calc] + + if isfile(path_solution) + # if the file exists run tests + proven_solution = YAML.load_file(path_solution) + + proven_x = Matrix(hcat([Float64.(row) for row in proven_solution["x"]]...)') + proven_p = Float64.(proven_solution["p"]) + + @test size(x_sorted) == size(proven_x) + @test vector_approx_test(vec(x_sorted), vec(proven_x)) + @test vector_approx_test(p_sorted, proven_p) + else + # otherwise create the tests + mkpath(dirname(path_solution)) + + dict_calc_solution = Dict( + "x" => [x_sorted[k, :] for k in 1:size(x_sorted, 1)], + "p" => p_sorted, + ) + + YAML.write_file(path_solution, dict_calc_solution) + @warn("Preloaded solution not found, then it has been created\nPath: $path_solution") + end + +end + list_examples = Examples.list_examples +list_multivariate_examples = Examples.list_multivariate_examples @testset "tests" begin for test in list_examples test_example(test.name, pem, test.d, test.N) end + + for test in list_multivariate_examples + test_multivariate_example(test.name, pem, test.d, test.N) + end end diff --git a/test/testcases/pem_multivariate/Normal2_N2.yml b/test/testcases/pem_multivariate/Normal2_N2.yml new file mode 100644 index 0000000..c83db20 --- /dev/null +++ b/test/testcases/pem_multivariate/Normal2_N2.yml @@ -0,0 +1,16 @@ +x: + - + - -1.0 + - -1.0 + - 1.0 + - 1.0 + - + - -1.0 + - 1.0 + - -1.0 + - 1.0 +p: + - 0.25 + - 0.25 + - 0.25 + - 0.25 diff --git a/test/testcases/pem_multivariate/Normal2_N3.yml b/test/testcases/pem_multivariate/Normal2_N3.yml new file mode 100644 index 0000000..b182c65 --- /dev/null +++ b/test/testcases/pem_multivariate/Normal2_N3.yml @@ -0,0 +1,31 @@ +x: + - + - -1.7320508075688776 + - -1.7320508075688776 + - -1.7320508075688776 + - 0.0 + - 0.0 + - 0.0 + - 1.7320508075688772 + - 1.7320508075688772 + - 1.7320508075688772 + - + - -1.7320508075688776 + - 0.0 + - 1.7320508075688772 + - -1.7320508075688776 + - 0.0 + - 1.7320508075688772 + - -1.7320508075688776 + - 0.0 + - 1.7320508075688772 +p: + - 0.027777777777777755 + - 0.11111111111111108 + - 0.027777777777777762 + - 0.11111111111111108 + - 0.44444444444444453 + - 0.1111111111111111 + - 0.027777777777777762 + - 0.1111111111111111 + - 0.027777777777777766 diff --git a/test/testcases/pem_multivariate/Normal_TruncNormal_N2.yml b/test/testcases/pem_multivariate/Normal_TruncNormal_N2.yml new file mode 100644 index 0000000..3f61610 --- /dev/null +++ b/test/testcases/pem_multivariate/Normal_TruncNormal_N2.yml @@ -0,0 +1,16 @@ +x: + - + - -1.0 + - -1.0 + - 1.0 + - 1.0 + - + - 0.6355198478077098 + - 1.4184929618467512 + - 0.6355198478077098 + - 1.4184929618467512 +p: + - 0.26274342969642217 + - 0.23725657030357783 + - 0.26274342969642217 + - 0.23725657030357783 diff --git a/test/testcases/pem_multivariate/Normal_TruncNormal_mixed_N.yml b/test/testcases/pem_multivariate/Normal_TruncNormal_mixed_N.yml new file mode 100644 index 0000000..f14aef2 --- /dev/null +++ b/test/testcases/pem_multivariate/Normal_TruncNormal_mixed_N.yml @@ -0,0 +1,22 @@ +x: + - + - -1.0 + - -1.0 + - -1.0 + - 1.0 + - 1.0 + - 1.0 + - + - 0.4089295792763922 + - 1.0532985603054232 + - 1.7295827634987933 + - 0.4089295792763922 + - 1.0532985603054232 + - 1.7295827634987933 +p: + - 0.10944862748508999 + - 0.3204603371606533 + - 0.07009103535425674 + - 0.10944862748508999 + - 0.3204603371606533 + - 0.07009103535425674 From 6b89d217809e232cfe3502b83985f765990506f0 Mon Sep 17 00:00:00 2001 From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com> Date: Mon, 20 Apr 2026 22:27:52 +0000 Subject: [PATCH 3/4] Fix review comments: remove unnecessary real() calls and fix doc default 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> --- src/pem.jl | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/pem.jl b/src/pem.jl index 32ef83c..ec8c6b8 100644 --- a/src/pem.jl +++ b/src/pem.jl @@ -231,7 +231,7 @@ Parameters Function used to calculate the mean value of each distribution - central_moment_fun :: Function (optional) Function used to calculate the central moment of each distribution -- montecarlo_sampling :: Integer (optional, default 1e6) +- montecarlo_sampling :: Integer (optional, default 1000000) Number of Monte Carlo samples used if a non-specific moment function is available - optimizer (optional) JuMP optimizer for executing the optimization @@ -279,9 +279,9 @@ function pem( for (j, combo) in enumerate(combos) for k in 1:K - x[k, j] = real(results[k].x[combo[k]]) + x[k, j] = results[k].x[combo[k]] end - p[j] = prod(real(results[k].p[combo[k]]) for k in 1:K) + p[j] = prod(results[k].p[combo[k]] for k in 1:K) end return (x=x, p=p) From 8698bbd4faef8fd2e3ce119a936acd251e4d1ab7 Mon Sep 17 00:00:00 2001 From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com> Date: Mon, 20 Apr 2026 22:54:00 +0000 Subject: [PATCH 4/4] Fix multivariate pem: use Hong (1998) concentration point scheme instead 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> --- src/pem.jl | 40 ++++++++++++------- test/runtests.jl | 2 +- .../testcases/pem_multivariate/Normal2_N2.yml | 8 ++-- .../testcases/pem_multivariate/Normal2_N3.yml | 25 ++++-------- .../Normal_TruncNormal_N2.yml | 12 +++--- .../Normal_TruncNormal_mixed_N.yml | 17 ++++---- 6 files changed, 52 insertions(+), 52 deletions(-) diff --git a/src/pem.jl b/src/pem.jl index ec8c6b8..6e34461 100644 --- a/src/pem.jl +++ b/src/pem.jl @@ -67,7 +67,7 @@ represented by the vector of elements d. Parameters ---------- -- d :: Vector +- d :: Vector{<:Real} Distribution under interest - N :: Integer Number of desired estimate points @@ -211,9 +211,14 @@ end Point Estimate Method to identify estimate points for K independent univariate distributions. -For K independent distributions, the multivariate representation is constructed by taking -the Cartesian product of the individual point estimates. Each combination has probability -equal to the product of the individual point probabilities. +For K independent distributions each with N concentrations, the multivariate representation +uses K*N total concentration points. Each concentration point for variable k at concentration n +has all other variables set to their means: + + x_{k,n} = (μ_1, ..., μ_{k-1}, x_{k,n}, μ_{k+1}, ..., μ_K) + +with weight p_{k,n} / K, where p_{k,n} is the univariate weight for distribution k at +concentration n. This ensures all weights sum to 1. This function implements the multivariate extension described in: - H.P. Hong, An efficient point estimate method for probabilistic analysis, @@ -240,10 +245,10 @@ Returns ------- - (x, p) :: NamedTuple - x :: Matrix{Float64} - Matrix of shape (K, M) where M = prod(N_1, ..., N_K). Each column is one - combination of point estimates across the K distributions. + Matrix of shape (K, M) where M = sum(N_1, ..., N_K). Each column is one + concentration point: only one variable deviates from its mean. - p :: Vector{Float64} - Joint probability of each combination point (product of individual probabilities). + Weight of each concentration point (p_{k,n} / K). """ function pem( @@ -270,18 +275,25 @@ function pem( for k in 1:K ] - # Form the Cartesian product of all point index combinations - combos = vec(collect(Iterators.product([1:Ns[k] for k in 1:K]...))) + # Mean of each distribution + μ = [mean_fun(distributions[k]) for k in 1:K] - M = length(combos) + # Build K*N concentration points following Hong (1998): + # For variable k at concentration n, all other variables stay at their mean. + M = sum(Ns) x = Matrix{Float64}(undef, K, M) p = Vector{Float64}(undef, M) - for (j, combo) in enumerate(combos) - for k in 1:K - x[k, j] = results[k].x[combo[k]] + j = 1 + for k in 1:K + for n in 1:Ns[k] + # All variables at their mean, except variable k + x[:, j] .= μ + x[k, j] = results[k].x[n] + # Weight: p_{k,n} / K + p[j] = results[k].p[n] / K + j += 1 end - p[j] = prod(results[k].p[combo[k]] for k in 1:K) end return (x=x, p=p) diff --git a/test/runtests.jl b/test/runtests.jl index 576a0e1..91e704a 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -62,7 +62,7 @@ function test_multivariate_example(example_name, testing_function, args...) 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)) + I_calc = sortperm(eachcol(calc_solution.x), by=col -> Tuple(col)) x_sorted = calc_solution.x[:, I_calc] p_sorted = calc_solution.p[I_calc] diff --git a/test/testcases/pem_multivariate/Normal2_N2.yml b/test/testcases/pem_multivariate/Normal2_N2.yml index c83db20..f1a6360 100644 --- a/test/testcases/pem_multivariate/Normal2_N2.yml +++ b/test/testcases/pem_multivariate/Normal2_N2.yml @@ -1,14 +1,14 @@ x: - - -1.0 - - -1.0 - - 1.0 + - 0.0 + - 0.0 - 1.0 - + - 0.0 - -1.0 - 1.0 - - -1.0 - - 1.0 + - 0.0 p: - 0.25 - 0.25 diff --git a/test/testcases/pem_multivariate/Normal2_N3.yml b/test/testcases/pem_multivariate/Normal2_N3.yml index b182c65..2c31f90 100644 --- a/test/testcases/pem_multivariate/Normal2_N3.yml +++ b/test/testcases/pem_multivariate/Normal2_N3.yml @@ -1,31 +1,22 @@ x: - - - -1.7320508075688776 - - -1.7320508075688776 - -1.7320508075688776 - 0.0 - 0.0 - 0.0 - - 1.7320508075688772 - - 1.7320508075688772 + - 0.0 - 1.7320508075688772 - - - -1.7320508075688776 - 0.0 - - 1.7320508075688772 - -1.7320508075688776 - 0.0 - - 1.7320508075688772 - - -1.7320508075688776 - 0.0 - 1.7320508075688772 + - 0.0 p: - - 0.027777777777777755 - - 0.11111111111111108 - - 0.027777777777777762 - - 0.11111111111111108 - - 0.44444444444444453 - - 0.1111111111111111 - - 0.027777777777777762 - - 0.1111111111111111 - - 0.027777777777777766 + - 0.0833333333333333 + - 0.0833333333333333 + - 0.33333333333333337 + - 0.33333333333333337 + - 0.08333333333333331 + - 0.08333333333333331 diff --git a/test/testcases/pem_multivariate/Normal_TruncNormal_N2.yml b/test/testcases/pem_multivariate/Normal_TruncNormal_N2.yml index 3f61610..60cf308 100644 --- a/test/testcases/pem_multivariate/Normal_TruncNormal_N2.yml +++ b/test/testcases/pem_multivariate/Normal_TruncNormal_N2.yml @@ -1,16 +1,16 @@ x: - - -1.0 - - -1.0 - - 1.0 + - 0.0 + - 0.0 - 1.0 - + - 1.0070551301947668 - 0.6355198478077098 - 1.4184929618467512 - - 0.6355198478077098 - - 1.4184929618467512 + - 1.0070551301947668 p: + - 0.25 - 0.26274342969642217 - 0.23725657030357783 - - 0.26274342969642217 - - 0.23725657030357783 + - 0.25 diff --git a/test/testcases/pem_multivariate/Normal_TruncNormal_mixed_N.yml b/test/testcases/pem_multivariate/Normal_TruncNormal_mixed_N.yml index f14aef2..593a5c6 100644 --- a/test/testcases/pem_multivariate/Normal_TruncNormal_mixed_N.yml +++ b/test/testcases/pem_multivariate/Normal_TruncNormal_mixed_N.yml @@ -1,22 +1,19 @@ x: - - -1.0 - - -1.0 - - -1.0 - - 1.0 - - 1.0 + - 0.0 + - 0.0 + - 0.0 - 1.0 - + - 1.0070551301947668 - 0.4089295792763922 - 1.0532985603054232 - 1.7295827634987933 - - 0.4089295792763922 - - 1.0532985603054232 - - 1.7295827634987933 + - 1.0070551301947668 p: + - 0.25 - 0.10944862748508999 - 0.3204603371606533 - 0.07009103535425674 - - 0.10944862748508999 - - 0.3204603371606533 - - 0.07009103535425674 + - 0.25