From 6093d7c81bd873396cdae874b078c69877572f08 Mon Sep 17 00:00:00 2001 From: Laurenz Knipper Date: Thu, 7 Aug 2025 10:40:53 +0200 Subject: [PATCH 01/15] sample based estimator for kucherenko indices --- src/UncertaintyQuantification.jl | 2 + src/sensitivity/kucherenkoindices.jl | 221 +++++++++++++++++++++++++++ 2 files changed, 223 insertions(+) create mode 100644 src/sensitivity/kucherenkoindices.jl diff --git a/src/UncertaintyQuantification.jl b/src/UncertaintyQuantification.jl index 4d6be345a..2deda7b6e 100644 --- a/src/UncertaintyQuantification.jl +++ b/src/UncertaintyQuantification.jl @@ -169,6 +169,7 @@ export quadrature_weights export rand export sample export sobolindices +export kucherenkoindices export to_copula_space export to_physical_space! export to_standard_normal_space @@ -221,6 +222,7 @@ include("simulations/importancesampling.jl") include("reliability/probabilityoffailure.jl") include("reliability/probabilityoffailure_imprecise.jl") include("sensitivity/sobolindices.jl") +include("sensitivity/kucherenkoindices.jl") include("util/fourier-transform.jl") include("util/wrap.jl") diff --git a/src/sensitivity/kucherenkoindices.jl b/src/sensitivity/kucherenkoindices.jl new file mode 100644 index 000000000..7f835bf2e --- /dev/null +++ b/src/sensitivity/kucherenkoindices.jl @@ -0,0 +1,221 @@ +const _kucherenko_table_types = [Symbol[], Float64[], Float64[]] +const _kucherenko_table_header = [ + "Variables", "FirstOrder", "TotalEffect" +] + +""" +- `X::Matrix`: Input sample matrix where each row is a sample point and each column is a variable +- `Y::Vector`: Output sample vector corresponding to the input samples +- `num_bins::Int=10`: Number of bins to use for conditioning +""" +function kucherenkoindices(X::Matrix, Y::Vector, num_bins::Int=10) + n_samples, n_vars = size(X) + + if length(Y) != n_samples + throw(ArgumentError("Number of output samples must match number of input samples")) + end + + indices = DataFrame(_kucherenko_table_types, _kucherenko_table_header) + total_var = var(Y) + + for i in 1:n_vars + S_i = _compute_first_order_kucherenko(X, Y, i, num_bins, total_var) + ST_i = _compute_total_effect_kucherenko(X, Y, i, num_bins, total_var) + + push!(indices, [Symbol("X$i"), S_i, ST_i]) + end + + return indices +end + +function _compute_first_order_kucherenko(X::Matrix, Y::Vector, var_idx::Int, num_bins::Int, total_var::Float64) + n_samples = length(Y) + + x_var = X[:, var_idx] + + x_min, x_max = extrema(x_var) + if x_min ≈ x_max + return 0.0 + end + + bin_edges = range(x_min, x_max, length=num_bins+1) + + bin_means = Float64[] + bin_weights = Float64[] + + for j in 1:num_bins + if j == num_bins + mask = (x_var .>= bin_edges[j]) .& (x_var .<= bin_edges[j+1]) + else + mask = (x_var .>= bin_edges[j]) .& (x_var .< bin_edges[j+1]) + end + + n_bin = sum(mask) + + if n_bin > 0 + bin_mean = mean(Y[mask]) + push!(bin_means, bin_mean) + push!(bin_weights, n_bin / n_samples) + end + end + + if isempty(bin_means) + return 0.0 + end + + weighted_mean = sum(bin_weights .* bin_means) + + S_i = sum(bin_weights .* (bin_means .- weighted_mean).^2) / total_var + + return S_i +end + +function _compute_total_effect_kucherenko(X::Matrix, Y::Vector, var_idx::Int, num_bins::Int, total_var::Float64) + n_samples, n_vars = size(X) + + if n_vars == 1 + return _compute_first_order_kucherenko(X, Y, var_idx, num_bins, total_var) + end + + other_vars = setdiff(1:n_vars, var_idx) + X_other = X[:, other_vars] + + bin_assignments = _assign_multidimensional_bins(X_other, num_bins) + + weighted_var_sum = 0.0 + total_weight = 0.0 + + for bin_id in unique(bin_assignments) + mask = bin_assignments .== bin_id + n_bin = sum(mask) + + if n_bin > 1 + bin_var = var(Y[mask]) + weight = n_bin / n_samples + weighted_var_sum += weight * bin_var + total_weight += weight + end + end + + ST_i = total_weight > 0 ? weighted_var_sum / total_var : 0.0 + + return ST_i +end + +function _assign_multidimensional_bins(X::Matrix, num_bins::Int) + n_samples, n_dims = size(X) + + X_norm = similar(X) + for j in 1:n_dims + x_min, x_max = extrema(X[:, j]) + if x_min ≈ x_max + X_norm[:, j] .= 0.5 + else + X_norm[:, j] = (X[:, j] .- x_min) ./ (x_max - x_min) + end + end + + bin_assignments = zeros(Int, n_samples) + + for i in 1:n_samples + bin_id = 0 + for j in 1:n_dims + bin_idx = min(floor(Int, X_norm[i, j] * num_bins), num_bins - 1) + bin_id += bin_idx * (num_bins ^ (j - 1)) + end + bin_assignments[i] = bin_id + 1 + end + + return bin_assignments +end + +""" +- `models::Vector{<:UQModel}`: Vector of UQ models to evaluate +- `inputs::Vector{<:UQInput}`: Vector of input distributions +- `outputs::Vector{Symbol}`: Vector of output quantity names +- `sim::AbstractMonteCarlo`: Monte Carlo simulation parameters +- `num_bins::Int=10`: Number of bins for conditioning (default: 10) +""" +function kucherenkoindices( + models::Vector{<:UQModel}, + inputs::Vector{<:UQInput}, + outputs::Vector{Symbol}, + sim::AbstractMonteCarlo; + num_bins::Int=10 +) + samples = sample(inputs, sim) + + evaluate!(models, samples) + + random_names = names(filter(i -> isa(i, RandomUQInput), inputs)) + + indices = Dict([ + (name, DataFrame(_kucherenko_table_types, _kucherenko_table_header)) for name in outputs + ]) + + X = Matrix(samples[:, random_names]) + + for output in outputs + Y = Vector(samples[:, output]) + + output_indices = kucherenkoindices(X, Y, num_bins) + + for (i, var_name) in enumerate(random_names) + output_indices[i, :Variables] = var_name + end + + indices[output] = output_indices + end + + return length(outputs) > 1 ? indices : indices[outputs[1]] +end + +function kucherenkoindices( + models::Vector{<:UQModel}, + inputs::UQInput, + outputs::Vector{Symbol}, + sim::AbstractMonteCarlo; + num_bins::Int=10 +) + return kucherenkoindices(models, [inputs], outputs, sim; num_bins=num_bins) +end + +function kucherenkoindices( + models::Vector{<:UQModel}, + inputs::Vector{<:UQInput}, + outputs::Symbol, + sim::AbstractMonteCarlo; + num_bins::Int=10 +) + return kucherenkoindices(models, inputs, [outputs], sim; num_bins=num_bins) +end + +function kucherenkoindices( + models::UQModel, + inputs::Vector{<:UQInput}, + outputs::Symbol, + sim::AbstractMonteCarlo; + num_bins::Int=10 +) + return kucherenkoindices([models], inputs, [outputs], sim; num_bins=num_bins) +end + +function kucherenkoindices( + models::Vector{<:UQModel}, + inputs::UQInput, + outputs::Symbol, + sim::AbstractMonteCarlo; + num_bins::Int=10 +) + return kucherenkoindices(models, [inputs], [outputs], sim; num_bins=num_bins) +end + +function kucherenkoindices( + models::UQModel, + inputs::UQInput, + outputs::Symbol, + sim::AbstractMonteCarlo; + num_bins::Int=10 +) + return kucherenkoindices([models], [inputs], [outputs], sim; num_bins=num_bins) +end From ceb175dafbb51c8e47c611c0c8d960624be4a665 Mon Sep 17 00:00:00 2001 From: Laurenz Knipper Date: Thu, 7 Aug 2025 10:45:50 +0200 Subject: [PATCH 02/15] Linear example from Kucherenko et al. (2012) --- demo/sensitivity/kucherenkoindices.jl | 68 +++++++++++++++++++++++++++ 1 file changed, 68 insertions(+) create mode 100644 demo/sensitivity/kucherenkoindices.jl diff --git a/demo/sensitivity/kucherenkoindices.jl b/demo/sensitivity/kucherenkoindices.jl new file mode 100644 index 000000000..d463e5361 --- /dev/null +++ b/demo/sensitivity/kucherenkoindices.jl @@ -0,0 +1,68 @@ +include("../../src/UncertaintyQuantification.jl") +using .UncertaintyQuantification +using LinearAlgebra + +# Testing Kucherenko indices - Test Case 1: Linear model with correlated variables from Kucherenko et al. (2012) + +ρ = 0.5 # correlation coefficient +σ = 2.0 # standard deviation for x3 + + +Σ = [1.0 0.0 0.0; + 0.0 1.0 ρ*σ; + 0.0 ρ*σ σ^2] + +D = sqrt.(diag(Σ)) +R = Σ ./ (D * D') + +marginals = UncertaintyQuantification.RandomVariable[ + UncertaintyQuantification.RandomVariable(Normal(0, 1), :x1), # std = 1 + UncertaintyQuantification.RandomVariable(Normal(0, 1), :x2), # std = 1 + UncertaintyQuantification.RandomVariable(Normal(0, σ), :x3) # std = σ +] + +inputs = [ + UncertaintyQuantification.JointDistribution(marginals, UncertaintyQuantification.GaussianCopula(R)) +] + +model = UncertaintyQuantification.Model(df -> df.x1 .+ df.x2 .+ df.x3, :y) + +sim = UncertaintyQuantification.MonteCarlo(500000) + +try + indices = UncertaintyQuantification.kucherenkoindices([model], inputs, [:y], sim, num_bins=60) + + println(indices) + + # Analytical calculations + denom = 2 + σ^2 + 2*ρ*σ + S1_analytical = 1 / denom + S2_analytical = (1 + ρ*σ)^2 / denom + S3_analytical = (σ + ρ)^2 / denom + ST1_analytical = 1 / denom + ST2_analytical = (1 - ρ^2) / denom + ST3_analytical = σ^2 * (1 - ρ^2) / denom + + + computed_values = Dict() + for i in 1:size(indices, 1) + var_name = string(indices[i, :Variables]) + computed_values[var_name] = (indices[i, :FirstOrder], indices[i, :TotalEffect]) + end + + tol = 0.01 + @assert abs(computed_values["x1"][1] - S1_analytical) < tol "S1 error too large" + @assert abs(computed_values["x2"][1] - S2_analytical) < tol "S2 error too large" + @assert abs(computed_values["x3"][1] - S3_analytical) < tol "S3 error too large" + @assert abs(computed_values["x1"][2] - ST1_analytical) < tol "ST1 error too large" + @assert abs(computed_values["x2"][2] - ST2_analytical) < tol "ST2 error too large" + @assert abs(computed_values["x3"][2] - ST3_analytical) < tol "ST3 error too large" + + println("Kucherenko indices validation passed - all values within tolerance.") + +catch e + println("Error computing Kucherenko indices: $e") + rethrow(e) +end + + From 62486f3be7a04e2b5d49c89e78f4fbdccbccdb61 Mon Sep 17 00:00:00 2001 From: Laurenz Knipper Date: Thu, 7 Aug 2025 11:04:17 +0200 Subject: [PATCH 03/15] removed unneccessary UQ namespaces --- demo/sensitivity/kucherenkoindices.jl | 16 ++++++++-------- 1 file changed, 8 insertions(+), 8 deletions(-) diff --git a/demo/sensitivity/kucherenkoindices.jl b/demo/sensitivity/kucherenkoindices.jl index d463e5361..1a3a63930 100644 --- a/demo/sensitivity/kucherenkoindices.jl +++ b/demo/sensitivity/kucherenkoindices.jl @@ -15,22 +15,22 @@ using LinearAlgebra D = sqrt.(diag(Σ)) R = Σ ./ (D * D') -marginals = UncertaintyQuantification.RandomVariable[ - UncertaintyQuantification.RandomVariable(Normal(0, 1), :x1), # std = 1 - UncertaintyQuantification.RandomVariable(Normal(0, 1), :x2), # std = 1 - UncertaintyQuantification.RandomVariable(Normal(0, σ), :x3) # std = σ +marginals = RandomVariable[ + RandomVariable(Normal(0, 1), :x1), # std = 1 + RandomVariable(Normal(0, 1), :x2), # std = 1 + RandomVariable(Normal(0, σ), :x3) # std = σ ] inputs = [ - UncertaintyQuantification.JointDistribution(marginals, UncertaintyQuantification.GaussianCopula(R)) + JointDistribution(marginals, GaussianCopula(R)) ] -model = UncertaintyQuantification.Model(df -> df.x1 .+ df.x2 .+ df.x3, :y) +model = Model(df -> df.x1 .+ df.x2 .+ df.x3, :y) -sim = UncertaintyQuantification.MonteCarlo(500000) +sim = MonteCarlo(500000) try - indices = UncertaintyQuantification.kucherenkoindices([model], inputs, [:y], sim, num_bins=60) + indices = kucherenkoindices([model], inputs, [:y], sim, num_bins=60) println(indices) From 5f1041d2d20c32043e91d969412f1e1f123cb28b Mon Sep 17 00:00:00 2001 From: Laurenz Knipper Date: Thu, 7 Aug 2025 14:57:02 +0200 Subject: [PATCH 04/15] included DOI --- demo/sensitivity/kucherenkoindices.jl | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/demo/sensitivity/kucherenkoindices.jl b/demo/sensitivity/kucherenkoindices.jl index 1a3a63930..82e3b3b59 100644 --- a/demo/sensitivity/kucherenkoindices.jl +++ b/demo/sensitivity/kucherenkoindices.jl @@ -2,7 +2,7 @@ include("../../src/UncertaintyQuantification.jl") using .UncertaintyQuantification using LinearAlgebra -# Testing Kucherenko indices - Test Case 1: Linear model with correlated variables from Kucherenko et al. (2012) +# Testing Kucherenko indices - Test Case 1: Linear model with correlated variables from Kucherenko et al. (2012) (DOI: 10.1016/j.cpc.2011.12.020) ρ = 0.5 # correlation coefficient σ = 2.0 # standard deviation for x3 @@ -16,9 +16,9 @@ D = sqrt.(diag(Σ)) R = Σ ./ (D * D') marginals = RandomVariable[ - RandomVariable(Normal(0, 1), :x1), # std = 1 - RandomVariable(Normal(0, 1), :x2), # std = 1 - RandomVariable(Normal(0, σ), :x3) # std = σ + RandomVariable(Normal(0, 1), :x1), + RandomVariable(Normal(0, 1), :x2), + RandomVariable(Normal(0, σ), :x3) ] inputs = [ From 654a9b968f7cb28426847d3e472da02a62577dc0 Mon Sep 17 00:00:00 2001 From: Laurenz Knipper Date: Thu, 7 Aug 2025 15:09:56 +0200 Subject: [PATCH 05/15] pushed test for single variable up so that first order indices do not get calculated twice --- src/sensitivity/kucherenkoindices.jl | 9 +++------ 1 file changed, 3 insertions(+), 6 deletions(-) diff --git a/src/sensitivity/kucherenkoindices.jl b/src/sensitivity/kucherenkoindices.jl index 7f835bf2e..4766ee34f 100644 --- a/src/sensitivity/kucherenkoindices.jl +++ b/src/sensitivity/kucherenkoindices.jl @@ -20,7 +20,8 @@ function kucherenkoindices(X::Matrix, Y::Vector, num_bins::Int=10) for i in 1:n_vars S_i = _compute_first_order_kucherenko(X, Y, i, num_bins, total_var) - ST_i = _compute_total_effect_kucherenko(X, Y, i, num_bins, total_var) + + ST_i = n_vars == 1 ? S_i : _compute_total_effect_kucherenko(X, Y, i, num_bins, total_var) push!(indices, [Symbol("X$i"), S_i, ST_i]) end @@ -73,10 +74,6 @@ end function _compute_total_effect_kucherenko(X::Matrix, Y::Vector, var_idx::Int, num_bins::Int, total_var::Float64) n_samples, n_vars = size(X) - if n_vars == 1 - return _compute_first_order_kucherenko(X, Y, var_idx, num_bins, total_var) - end - other_vars = setdiff(1:n_vars, var_idx) X_other = X[:, other_vars] @@ -133,7 +130,7 @@ end - `models::Vector{<:UQModel}`: Vector of UQ models to evaluate - `inputs::Vector{<:UQInput}`: Vector of input distributions - `outputs::Vector{Symbol}`: Vector of output quantity names -- `sim::AbstractMonteCarlo`: Monte Carlo simulation parameters +- `sim::AbstractMonteCarlo`: Total MC samples - `num_bins::Int=10`: Number of bins for conditioning (default: 10) """ function kucherenkoindices( From 5e3d5455df19af5caca01e29b6fb88a8e728260f Mon Sep 17 00:00:00 2001 From: Laurenz Knipper Date: Tue, 12 Aug 2025 11:15:43 +0200 Subject: [PATCH 06/15] added functionality to generate conditional samples based on gaussian copula --- src/inputs/jointdistribution.jl | 53 +++++++++++++++++++++++++++++++++ 1 file changed, 53 insertions(+) diff --git a/src/inputs/jointdistribution.jl b/src/inputs/jointdistribution.jl index c7a51504a..7a77e0674 100644 --- a/src/inputs/jointdistribution.jl +++ b/src/inputs/jointdistribution.jl @@ -22,6 +22,59 @@ function sample(jd::JointDistribution, n::Integer=1) return samples end + +function sample_conditional_copula(joint::JointDistribution, var_names::Vector{Symbol}, x_values::Vector{Float64}, N::Int) + length(var_names) != length(x_values) && error("Number of variable names ($(length(var_names))) must match number of values ($(length(x_values)))") + !isa(joint.copula, GaussianCopula) && error("This function for conditional sampling is only implemented for Gaussian copulas.") + + marginals, copula, R = joint.marginals, joint.copula, joint.copula.correlation + all_var_names, d = [marginal.name for marginal in marginals], length(marginals) + + v_indices = Int[] + for var_name in var_names + idx = findfirst(==(var_name), all_var_names) + idx === nothing && error("Variable $var_name not found in joint distribution. Available variables: $all_var_names") + push!(v_indices, idx) + end + + length(v_indices) >= d && error("All $(d) variables are fixed - need at least one free variable to sample") + + w_indices = setdiff(1:d, v_indices) + + z_v = zeros(length(v_indices)) + for (i, (idx, x_val)) in enumerate(zip(v_indices, x_values)) + z_v[i] = quantile(Normal(0,1), cdf(marginals[idx].dist, x_val)) + end + + Σ_vv, Σ_wv, Σ_vw, Σ_ww = R[v_indices, v_indices], R[w_indices, v_indices], R[v_indices, w_indices], R[w_indices, w_indices] + + if length(v_indices) == 1 + μ_cond = (Σ_wv / Σ_vv[1,1]) * z_v[1] + Σ_cond = Σ_ww .- (Σ_wv * Σ_vw) / Σ_vv[1,1] + else + Σ_vv_inv = inv(Σ_vv) + μ_cond = Σ_wv * Σ_vv_inv * z_v + Σ_cond = Σ_ww .- Σ_wv * Σ_vv_inv * Σ_vw + end + + μ_cond = vec(μ_cond) + Z_w = rand(MvNormal(μ_cond, Symmetric(Σ_cond)), N)' + + samples = DataFrame() + for (i, var_name) in enumerate(var_names) + samples[!, var_name] = fill(x_values[i], N) + end + for (i, idx) in enumerate(w_indices) + samples[!, all_var_names[idx]] = quantile.(marginals[idx].dist, cdf.(Normal(), Z_w[:,i])) + end + + return samples +end + +function sample_conditional_copula(joint::JointDistribution, var_name::Symbol, x_v::Float64, N::Int) + return sample_conditional_copula(joint, [var_name], [x_v], N) +end + function to_physical_space!(jd::JointDistribution, x::DataFrame) correlated_cdf = to_copula_space(jd.copula, Matrix{Float64}(x[:, names(jd)])) for (i, rv) in enumerate(jd.marginals) From 7467270f241c79aaae0372e9ec7a3b9ad02ee071 Mon Sep 17 00:00:00 2001 From: Laurenz Knipper Date: Tue, 12 Aug 2025 11:17:43 +0200 Subject: [PATCH 07/15] Added conditional sampling based Kucherenko indices - old method now named kucherenkoindices_bin --- demo/sensitivity/kucherenkoindices.jl | 42 ++++- src/UncertaintyQuantification.jl | 2 + src/sensitivity/kucherenkoindices.jl | 212 +++++++++++++++++++++++--- 3 files changed, 235 insertions(+), 21 deletions(-) diff --git a/demo/sensitivity/kucherenkoindices.jl b/demo/sensitivity/kucherenkoindices.jl index 82e3b3b59..84ca037ea 100644 --- a/demo/sensitivity/kucherenkoindices.jl +++ b/demo/sensitivity/kucherenkoindices.jl @@ -27,10 +27,11 @@ inputs = [ model = Model(df -> df.x1 .+ df.x2 .+ df.x3, :y) -sim = MonteCarlo(500000) +sim = MonteCarlo(50000) +println("Using Bins") try - indices = kucherenkoindices([model], inputs, [:y], sim, num_bins=60) + indices = kucherenkoindices_bin([model], inputs, [:y], sim, num_bins=60) println(indices) @@ -66,3 +67,40 @@ catch e end +println("\nUsing conditional sampling") + +try + indices = kucherenkoindices([model], inputs, [:y], sim) + + println(indices) + + # Analytical calculations + denom = 2 + σ^2 + 2*ρ*σ + S1_analytical = 1 / denom + S2_analytical = (1 + ρ*σ)^2 / denom + S3_analytical = (σ + ρ)^2 / denom + ST1_analytical = 1 / denom + ST2_analytical = (1 - ρ^2) / denom + ST3_analytical = σ^2 * (1 - ρ^2) / denom + + + computed_values = Dict() + for i in 1:size(indices, 1) + var_name = string(indices[i, :Variables]) + computed_values[var_name] = (indices[i, :FirstOrder], indices[i, :TotalEffect]) + end + + tol = 0.01 + @assert abs(computed_values["x1"][1] - S1_analytical) < tol "S1 error too large" + @assert abs(computed_values["x2"][1] - S2_analytical) < tol "S2 error too large" + @assert abs(computed_values["x3"][1] - S3_analytical) < tol "S3 error too large" + @assert abs(computed_values["x1"][2] - ST1_analytical) < tol "ST1 error too large" + @assert abs(computed_values["x2"][2] - ST2_analytical) < tol "ST2 error too large" + @assert abs(computed_values["x3"][2] - ST3_analytical) < tol "ST3 error too large" + + println("Kucherenko indices validation passed - all values within tolerance.") + +catch e + println("Error computing Kucherenko indices: $e") + rethrow(e) +end diff --git a/src/UncertaintyQuantification.jl b/src/UncertaintyQuantification.jl index 2deda7b6e..7ed147052 100644 --- a/src/UncertaintyQuantification.jl +++ b/src/UncertaintyQuantification.jl @@ -168,8 +168,10 @@ export quadrature_nodes export quadrature_weights export rand export sample +export sample_conditional_copula export sobolindices export kucherenkoindices +export kucherenkoindices_bin export to_copula_space export to_physical_space! export to_standard_normal_space diff --git a/src/sensitivity/kucherenkoindices.jl b/src/sensitivity/kucherenkoindices.jl index 4766ee34f..41ec73f7a 100644 --- a/src/sensitivity/kucherenkoindices.jl +++ b/src/sensitivity/kucherenkoindices.jl @@ -3,12 +3,176 @@ const _kucherenko_table_header = [ "Variables", "FirstOrder", "TotalEffect" ] +""" +- `models::Vector{<:UQModel}`: Vector of UQ models to evaluate +- `inputs::Vector{<:UQInput}`: Vector of input distributions +- `output::Symbol`: Output quantity name +- `sim::AbstractMonteCarlo`: == N : Total MC samples = N*(2M+1) +""" +function kucherenkoindices( + models::Vector{<:UQModel}, + inputs::Vector{<:UQInput}, + output::Symbol, + sim::AbstractMonteCarlo +) + indices = DataFrame(_kucherenko_table_types, _kucherenko_table_header) + + samples = sample(inputs, sim) + evaluate!(models, samples) + + random_names = names(filter(i -> isa(i, RandomUQInput), inputs)) + + Y_orig = Vector(samples[:, output]) + total_var = var(Y_orig) + + for i in random_names + + i_cond_samples = _generate_conditional_samples(samples, inputs[1], [i]) + evaluate!(models, i_cond_samples) + S_i = _compute_first_order_kucherenko(samples, i_cond_samples, output, total_var) + + other_vars = setdiff(random_names, [i]) + other_cond_samples = _generate_conditional_samples(samples, inputs[1], other_vars) + evaluate!(models, other_cond_samples) + ST_i = length(random_names) == 1 ? S_i : _compute_total_effect_kucherenko(samples, other_cond_samples, output, total_var) + + push!(indices, [i, S_i, ST_i]) + end + + return indices +end + + +function _compute_first_order_kucherenko( + samples::DataFrame, + cond_samples::DataFrame, + output::Symbol, + total_var::Float64 +) + Y_orig = Vector(samples[:, output]) + Y_cond = Vector(cond_samples[:, output]) + + S_i = (mean(Y_orig .* Y_cond) - mean(Y_orig)^2) / total_var # Kucherenko et. al. 2012 Eq. 5.3 + + return S_i +end + +function _compute_total_effect_kucherenko( + samples::DataFrame, + cond_samples::DataFrame, + output::Symbol, + total_var::Float64 +) + Y_orig = Vector(samples[:, output]) + Y_cond = Vector(cond_samples[:, output]) + + ST_i = mean((Y_orig .- Y_cond).^2) / (2 * total_var) # Kucherenko et. al. 2012 Eq. 5.4 + + return ST_i +end + +function _generate_conditional_samples( + samples::DataFrame, + joint_dist::JointDistribution, + var_names::Vector{Symbol} +) + input_var_names = [marginal.name for marginal in joint_dist.marginals] + conditional_samples = map(1:nrow(samples)) do i + x_values = [samples[i, var_name] for var_name in var_names] + sample_conditional_copula(joint_dist, var_names, x_values, 1) + end + + result = vcat(conditional_samples...) + return select(result, input_var_names) +end + +""" +- `models::Vector{<:UQModel}`: Vector of UQ models to evaluate +- `inputs::Vector{<:UQInput}`: Vector of input distributions +- `outputs::Vector{Symbol}`: Vector of output quantity names +- `sim::AbstractMonteCarlo`: == N : Total MC samples = N*(2M+1) +""" +function kucherenkoindices( + models::Vector{<:UQModel}, + inputs::Vector{<:UQInput}, + outputs::Vector{Symbol}, + sim::AbstractMonteCarlo +) + + random_names = names(filter(i -> isa(i, RandomUQInput), inputs)) + + indices = Dict([ + (name, DataFrame(_kucherenko_table_types, _kucherenko_table_header)) for name in outputs + ]) + + for output in outputs + + output_indices = kucherenkoindices(models, inputs, output, sim) + + for (i, var_name) in enumerate(random_names) + output_indices[i, :Variables] = var_name + end + + indices[output] = output_indices + end + + return length(outputs) > 1 ? indices : indices[outputs[1]] +end + +function kucherenkoindices( + models::Vector{<:UQModel}, + inputs::UQInput, + outputs::Vector{Symbol}, + sim::AbstractMonteCarlo +) + return kucherenkoindices(models, [inputs], outputs, sim) +end + + +function kucherenkoindices( + models::UQModel, + inputs::Vector{<:UQInput}, + outputs::Symbol, + sim::AbstractMonteCarlo +) + return kucherenkoindices([models], inputs, [outputs], sim) +end + +function kucherenkoindices( + models::Vector{<:UQModel}, + inputs::UQInput, + outputs::Symbol, + sim::AbstractMonteCarlo +) + return kucherenkoindices(models, [inputs], [outputs], sim) +end + +function kucherenkoindices( + models::UQModel, + inputs::UQInput, + outputs::Symbol, + sim::AbstractMonteCarlo +) + return kucherenkoindices([models], [inputs], [outputs], sim) +end + +function kucherenkoindices( + models::UQModel, + inputs::Vector{<:UQInput}, + outputs::Vector{Symbol}, + sim::AbstractMonteCarlo +) + return kucherenkoindices([models], inputs, outputs, sim) +end + + + """ - `X::Matrix`: Input sample matrix where each row is a sample point and each column is a variable - `Y::Vector`: Output sample vector corresponding to the input samples - `num_bins::Int=10`: Number of bins to use for conditioning """ -function kucherenkoindices(X::Matrix, Y::Vector, num_bins::Int=10) +function kucherenkoindices_bin(X::Matrix, Y::Vector, num_bins::Int=10) n_samples, n_vars = size(X) if length(Y) != n_samples @@ -19,9 +183,9 @@ function kucherenkoindices(X::Matrix, Y::Vector, num_bins::Int=10) total_var = var(Y) for i in 1:n_vars - S_i = _compute_first_order_kucherenko(X, Y, i, num_bins, total_var) - - ST_i = n_vars == 1 ? S_i : _compute_total_effect_kucherenko(X, Y, i, num_bins, total_var) + S_i = _compute_first_order_kucherenko_bins(X, Y, i, num_bins, total_var) + + ST_i = n_vars == 1 ? S_i : _compute_total_effect_kucherenko_bins(X, Y, i, num_bins, total_var) push!(indices, [Symbol("X$i"), S_i, ST_i]) end @@ -29,7 +193,7 @@ function kucherenkoindices(X::Matrix, Y::Vector, num_bins::Int=10) return indices end -function _compute_first_order_kucherenko(X::Matrix, Y::Vector, var_idx::Int, num_bins::Int, total_var::Float64) +function _compute_first_order_kucherenko_bins(X::Matrix, Y::Vector, var_idx::Int, num_bins::Int, total_var::Float64) n_samples = length(Y) x_var = X[:, var_idx] @@ -71,7 +235,7 @@ function _compute_first_order_kucherenko(X::Matrix, Y::Vector, var_idx::Int, num return S_i end -function _compute_total_effect_kucherenko(X::Matrix, Y::Vector, var_idx::Int, num_bins::Int, total_var::Float64) +function _compute_total_effect_kucherenko_bins(X::Matrix, Y::Vector, var_idx::Int, num_bins::Int, total_var::Float64) n_samples, n_vars = size(X) other_vars = setdiff(1:n_vars, var_idx) @@ -131,9 +295,9 @@ end - `inputs::Vector{<:UQInput}`: Vector of input distributions - `outputs::Vector{Symbol}`: Vector of output quantity names - `sim::AbstractMonteCarlo`: Total MC samples -- `num_bins::Int=10`: Number of bins for conditioning (default: 10) +- `num_bins::Int`: Number of bins for conditioning """ -function kucherenkoindices( +function kucherenkoindices_bin( models::Vector{<:UQModel}, inputs::Vector{<:UQInput}, outputs::Vector{Symbol}, @@ -155,7 +319,7 @@ function kucherenkoindices( for output in outputs Y = Vector(samples[:, output]) - output_indices = kucherenkoindices(X, Y, num_bins) + output_indices = kucherenkoindices_bin(X, Y, num_bins) for (i, var_name) in enumerate(random_names) output_indices[i, :Variables] = var_name @@ -167,52 +331,62 @@ function kucherenkoindices( return length(outputs) > 1 ? indices : indices[outputs[1]] end -function kucherenkoindices( +function kucherenkoindices_bin( models::Vector{<:UQModel}, inputs::UQInput, outputs::Vector{Symbol}, sim::AbstractMonteCarlo; num_bins::Int=10 ) - return kucherenkoindices(models, [inputs], outputs, sim; num_bins=num_bins) + return kucherenkoindices_bin(models, [inputs], outputs, sim; num_bins=num_bins) end -function kucherenkoindices( +function kucherenkoindices_bin( models::Vector{<:UQModel}, inputs::Vector{<:UQInput}, outputs::Symbol, sim::AbstractMonteCarlo; num_bins::Int=10 ) - return kucherenkoindices(models, inputs, [outputs], sim; num_bins=num_bins) + return kucherenkoindices_bin(models, inputs, [outputs], sim; num_bins=num_bins) end -function kucherenkoindices( +function kucherenkoindices_bin( models::UQModel, inputs::Vector{<:UQInput}, outputs::Symbol, sim::AbstractMonteCarlo; num_bins::Int=10 ) - return kucherenkoindices([models], inputs, [outputs], sim; num_bins=num_bins) + return kucherenkoindices_bin([models], inputs, [outputs], sim; num_bins=num_bins) end -function kucherenkoindices( +function kucherenkoindices_bin( models::Vector{<:UQModel}, inputs::UQInput, outputs::Symbol, sim::AbstractMonteCarlo; num_bins::Int=10 ) - return kucherenkoindices(models, [inputs], [outputs], sim; num_bins=num_bins) + return kucherenkoindices_bin(models, [inputs], [outputs], sim; num_bins=num_bins) end -function kucherenkoindices( +function kucherenkoindices_bin( models::UQModel, inputs::UQInput, outputs::Symbol, sim::AbstractMonteCarlo; num_bins::Int=10 ) - return kucherenkoindices([models], [inputs], [outputs], sim; num_bins=num_bins) + return kucherenkoindices_bin([models], [inputs], [outputs], sim; num_bins=num_bins) +end + +function kucherenkoindices_bin( + models::UQModel, + inputs::Vector{<:UQInput}, + outputs::Vector{Symbol}, + sim::AbstractMonteCarlo; + num_bins::Int=10 +) + return kucherenkoindices_bin([models], inputs, outputs, sim; num_bins=num_bins) end From 03c82bd224864f58377fee19d8793e95d300606c Mon Sep 17 00:00:00 2001 From: sitoryu Date: Tue, 12 Aug 2025 13:39:39 +0200 Subject: [PATCH 08/15] kucherenkoindices_bin pass calculated samples as output --- demo/sensitivity/kucherenkoindices.jl | 2 +- src/sensitivity/kucherenkoindices.jl | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/demo/sensitivity/kucherenkoindices.jl b/demo/sensitivity/kucherenkoindices.jl index 84ca037ea..39e1b1fc2 100644 --- a/demo/sensitivity/kucherenkoindices.jl +++ b/demo/sensitivity/kucherenkoindices.jl @@ -31,7 +31,7 @@ sim = MonteCarlo(50000) println("Using Bins") try - indices = kucherenkoindices_bin([model], inputs, [:y], sim, num_bins=60) + indices, bin_samples = kucherenkoindices_bin([model], inputs, [:y], sim, num_bins=60) println(indices) diff --git a/src/sensitivity/kucherenkoindices.jl b/src/sensitivity/kucherenkoindices.jl index 41ec73f7a..8d7195bf1 100644 --- a/src/sensitivity/kucherenkoindices.jl +++ b/src/sensitivity/kucherenkoindices.jl @@ -328,7 +328,7 @@ function kucherenkoindices_bin( indices[output] = output_indices end - return length(outputs) > 1 ? indices : indices[outputs[1]] + return length(outputs) > 1 ? indices : indices[outputs[1]], samples end function kucherenkoindices_bin( From 2e111ab2b2cabfe9d7e8716e3b4414162b33cfcd Mon Sep 17 00:00:00 2001 From: sitoryu Date: Thu, 21 Aug 2025 11:26:33 +0200 Subject: [PATCH 09/15] added quantile bins and sperated 1D and MultiD --- demo/sensitivity/kucherenkoindices.jl | 4 +- src/sensitivity/kucherenkoindices.jl | 87 ++++++++++++++------------- 2 files changed, 46 insertions(+), 45 deletions(-) diff --git a/demo/sensitivity/kucherenkoindices.jl b/demo/sensitivity/kucherenkoindices.jl index 39e1b1fc2..03937cd49 100644 --- a/demo/sensitivity/kucherenkoindices.jl +++ b/demo/sensitivity/kucherenkoindices.jl @@ -27,11 +27,11 @@ inputs = [ model = Model(df -> df.x1 .+ df.x2 .+ df.x3, :y) -sim = MonteCarlo(50000) +sim = MonteCarlo(100000) println("Using Bins") try - indices, bin_samples = kucherenkoindices_bin([model], inputs, [:y], sim, num_bins=60) + indices, bin_samples = kucherenkoindices_bin([model], inputs, [:y], sim; min_bin_sample=2000, min_bin_sample_multi_dims=10) println(indices) diff --git a/src/sensitivity/kucherenkoindices.jl b/src/sensitivity/kucherenkoindices.jl index 8d7195bf1..ac23b4e7e 100644 --- a/src/sensitivity/kucherenkoindices.jl +++ b/src/sensitivity/kucherenkoindices.jl @@ -170,14 +170,17 @@ end """ - `X::Matrix`: Input sample matrix where each row is a sample point and each column is a variable - `Y::Vector`: Output sample vector corresponding to the input samples -- `num_bins::Int=10`: Number of bins to use for conditioning +- `min_bin_sample::Int=25`: Minimum samples per bin for 1 dimension for conditioning; Recommended amount is at least 25 samples per bin with bin amount around 100 +- `min_bin_sample_multi_dims::Int=25`: Minimum samples per bin-dimension in multiple dimensions for conditioning; Recommended amount is about 10-25 samples per bin """ -function kucherenkoindices_bin(X::Matrix, Y::Vector, num_bins::Int=10) +function kucherenkoindices_bin(X::Matrix, Y::Vector; min_bin_sample=nothing, min_bin_sample_multi_dims::Int=25) n_samples, n_vars = size(X) + + num_bins = min_bin_sample === nothing ? min(100, floor(Int, n_samples / 25)) : floor(Int, n_samples / min_bin_sample) + num_bins_multi = floor(Int, n_samples / min_bin_sample_multi_dims) - if length(Y) != n_samples - throw(ArgumentError("Number of output samples must match number of input samples")) - end + if num_bins > 500 @warn "More than 500 bins in single dimension: $num_bins" end + if length(Y) != n_samples throw(ArgumentError("Number of output samples must match number of input samples")) end indices = DataFrame(_kucherenko_table_types, _kucherenko_table_header) total_var = var(Y) @@ -185,7 +188,7 @@ function kucherenkoindices_bin(X::Matrix, Y::Vector, num_bins::Int=10) for i in 1:n_vars S_i = _compute_first_order_kucherenko_bins(X, Y, i, num_bins, total_var) - ST_i = n_vars == 1 ? S_i : _compute_total_effect_kucherenko_bins(X, Y, i, num_bins, total_var) + ST_i = n_vars == 1 ? S_i : _compute_total_effect_kucherenko_bins(X, Y, i, num_bins_multi, total_var) push!(indices, [Symbol("X$i"), S_i, ST_i]) end @@ -195,15 +198,10 @@ end function _compute_first_order_kucherenko_bins(X::Matrix, Y::Vector, var_idx::Int, num_bins::Int, total_var::Float64) n_samples = length(Y) - x_var = X[:, var_idx] - - x_min, x_max = extrema(x_var) - if x_min ≈ x_max - return 0.0 - end - - bin_edges = range(x_min, x_max, length=num_bins+1) + + quantiles = range(0, 1; length=num_bins+1) + bin_edges = quantile(x_var, quantiles) bin_means = Float64[] bin_weights = Float64[] @@ -265,24 +263,20 @@ end function _assign_multidimensional_bins(X::Matrix, num_bins::Int) n_samples, n_dims = size(X) - - X_norm = similar(X) - for j in 1:n_dims - x_min, x_max = extrema(X[:, j]) - if x_min ≈ x_max - X_norm[:, j] .= 0.5 - else - X_norm[:, j] = (X[:, j] .- x_min) ./ (x_max - x_min) - end - end + bins_per_dim = floor(Int, num_bins^(1/n_dims)) + + quantile_edges = [quantile(X[:, j], range(0, 1; length=bins_per_dim+1)) for j in 1:n_dims] bin_assignments = zeros(Int, n_samples) for i in 1:n_samples bin_id = 0 for j in 1:n_dims - bin_idx = min(floor(Int, X_norm[i, j] * num_bins), num_bins - 1) - bin_id += bin_idx * (num_bins ^ (j - 1)) + edges = quantile_edges[j] + x = X[i, j] + bin_idx = searchsortedlast(edges, x) + bin_idx = clamp(bin_idx, 1, bins_per_dim) + bin_id += (bin_idx - 1) * (bins_per_dim ^ (j - 1)) end bin_assignments[i] = bin_id + 1 end @@ -294,18 +288,19 @@ end - `models::Vector{<:UQModel}`: Vector of UQ models to evaluate - `inputs::Vector{<:UQInput}`: Vector of input distributions - `outputs::Vector{Symbol}`: Vector of output quantity names -- `sim::AbstractMonteCarlo`: Total MC samples -- `num_bins::Int`: Number of bins for conditioning +- `sim::AbstractMonteCarlo`: Total MC samples +- `min_bin_sample::Int=25`: Minimum samples per bin for 1 dimension for conditioning; Recommended amount is at least 25 samples per bin but bin amount not higher than 100 +- `min_bin_sample_multi_dims::Int=25`: Minimum samples per bin-dimension in multiple dimensions for conditioning; Recommended amount is about 10-25 samples per bin """ function kucherenkoindices_bin( models::Vector{<:UQModel}, inputs::Vector{<:UQInput}, outputs::Vector{Symbol}, sim::AbstractMonteCarlo; - num_bins::Int=10 + min_bin_sample=nothing, + min_bin_sample_multi_dims::Int=25 ) samples = sample(inputs, sim) - evaluate!(models, samples) random_names = names(filter(i -> isa(i, RandomUQInput), inputs)) @@ -319,7 +314,7 @@ function kucherenkoindices_bin( for output in outputs Y = Vector(samples[:, output]) - output_indices = kucherenkoindices_bin(X, Y, num_bins) + output_indices = kucherenkoindices_bin(X, Y; min_bin_sample = min_bin_sample, min_bin_sample_multi_dims = min_bin_sample_multi_dims) for (i, var_name) in enumerate(random_names) output_indices[i, :Variables] = var_name @@ -336,9 +331,10 @@ function kucherenkoindices_bin( inputs::UQInput, outputs::Vector{Symbol}, sim::AbstractMonteCarlo; - num_bins::Int=10 + min_bin_sample=nothing, + min_bin_sample_multi_dims::Int=25 ) - return kucherenkoindices_bin(models, [inputs], outputs, sim; num_bins=num_bins) + return kucherenkoindices_bin(models, [inputs], outputs, sim; min_bin_sample=min_bin_sample, min_bin_sample_multi_dims=min_bin_sample_multi_dims) end function kucherenkoindices_bin( @@ -346,9 +342,10 @@ function kucherenkoindices_bin( inputs::Vector{<:UQInput}, outputs::Symbol, sim::AbstractMonteCarlo; - num_bins::Int=10 + min_bin_sample=nothing, + min_bin_sample_multi_dims::Int=25 ) - return kucherenkoindices_bin(models, inputs, [outputs], sim; num_bins=num_bins) + return kucherenkoindices_bin(models, inputs, [outputs], sim; min_bin_sample=min_bin_sample, min_bin_sample_multi_dims=min_bin_sample_multi_dims) end function kucherenkoindices_bin( @@ -356,9 +353,10 @@ function kucherenkoindices_bin( inputs::Vector{<:UQInput}, outputs::Symbol, sim::AbstractMonteCarlo; - num_bins::Int=10 + min_bin_sample=nothing, + min_bin_sample_multi_dims::Int=25 ) - return kucherenkoindices_bin([models], inputs, [outputs], sim; num_bins=num_bins) + return kucherenkoindices_bin([models], inputs, [outputs], sim; min_bin_sample=min_bin_sample, min_bin_sample_multi_dims=min_bin_sample_multi_dims) end function kucherenkoindices_bin( @@ -366,9 +364,10 @@ function kucherenkoindices_bin( inputs::UQInput, outputs::Symbol, sim::AbstractMonteCarlo; - num_bins::Int=10 + min_bin_sample=nothing, + min_bin_sample_multi_dims::Int=25 ) - return kucherenkoindices_bin(models, [inputs], [outputs], sim; num_bins=num_bins) + return kucherenkoindices_bin(models, [inputs], [outputs], sim; min_bin_sample=min_bin_sample, min_bin_sample_multi_dims=min_bin_sample_multi_dims) end function kucherenkoindices_bin( @@ -376,9 +375,10 @@ function kucherenkoindices_bin( inputs::UQInput, outputs::Symbol, sim::AbstractMonteCarlo; - num_bins::Int=10 + min_bin_sample=nothing, + min_bin_sample_multi_dims::Int=25 ) - return kucherenkoindices_bin([models], [inputs], [outputs], sim; num_bins=num_bins) + return kucherenkoindices_bin([models], [inputs], [outputs], sim; min_bin_sample=min_bin_sample, min_bin_sample_multi_dims=min_bin_sample_multi_dims) end function kucherenkoindices_bin( @@ -386,7 +386,8 @@ function kucherenkoindices_bin( inputs::Vector{<:UQInput}, outputs::Vector{Symbol}, sim::AbstractMonteCarlo; - num_bins::Int=10 + min_bin_sample=nothing, + min_bin_sample_multi_dims::Int=25 ) - return kucherenkoindices_bin([models], inputs, outputs, sim; num_bins=num_bins) + return kucherenkoindices_bin([models], inputs, outputs, sim; min_bin_sample=min_bin_sample, min_bin_sample_multi_dims=min_bin_sample_multi_dims) end From 1bd0f239c6d83c60ce67503f6f045edfa488f97a Mon Sep 17 00:00:00 2001 From: sitoryu Date: Mon, 25 Aug 2025 15:32:30 +0200 Subject: [PATCH 10/15] changed Demo example to Portfolio model --- demo/sensitivity/kucherenkoindices.jl | 148 +++++++++++++------------- 1 file changed, 72 insertions(+), 76 deletions(-) diff --git a/demo/sensitivity/kucherenkoindices.jl b/demo/sensitivity/kucherenkoindices.jl index 03937cd49..5f7243d8c 100644 --- a/demo/sensitivity/kucherenkoindices.jl +++ b/demo/sensitivity/kucherenkoindices.jl @@ -2,105 +2,101 @@ include("../../src/UncertaintyQuantification.jl") using .UncertaintyQuantification using LinearAlgebra -# Testing Kucherenko indices - Test Case 1: Linear model with correlated variables from Kucherenko et al. (2012) (DOI: 10.1016/j.cpc.2011.12.020) +# --- Test Case 2: Portfolio model with analytical indices --- -ρ = 0.5 # correlation coefficient -σ = 2.0 # standard deviation for x3 +μ = [0.0, 0.0, 250.0, 400.0] +σ1, σ2, σ3, σ4 = sqrt(16.0), sqrt(4.0), sqrt(4e4), sqrt(9e4) +σ12, σ21 = 2.4, 2.4 +σ34, σ43 = -1.8e4, -1.8e4 -Σ = [1.0 0.0 0.0; - 0.0 1.0 ρ*σ; - 0.0 ρ*σ σ^2] - -D = sqrt.(diag(Σ)) -R = Σ ./ (D * D') +Σ = [ + 16.0 2.4 0.0 0.0; + 2.4 4.0 0.0 0.0; + 0.0 0.0 4e4 -1.8e4; + 0.0 0.0 -1.8e4 9e4 +] +# Marginals marginals = RandomVariable[ - RandomVariable(Normal(0, 1), :x1), - RandomVariable(Normal(0, 1), :x2), - RandomVariable(Normal(0, σ), :x3) + RandomVariable(Normal(μ[1], σ1), :x1), + RandomVariable(Normal(μ[2], σ2), :x2), + RandomVariable(Normal(μ[3], σ3), :x3), + RandomVariable(Normal(μ[4], σ4), :x4) ] +# Correlation matrix for Gaussian copula +Dvec = sqrt.(diag(Σ)) +R = Σ ./ (Dvec * Dvec') + inputs = [ JointDistribution(marginals, GaussianCopula(R)) ] -model = Model(df -> df.x1 .+ df.x2 .+ df.x3, :y) +model = Model(df -> df.x1 .* df.x3 .+ df.x2 .* df.x4, :y) +sim = MonteCarlo(200000) -sim = MonteCarlo(100000) -println("Using Bins") -try - indices, bin_samples = kucherenkoindices_bin([model], inputs, [:y], sim; min_bin_sample=2000, min_bin_sample_multi_dims=10) +# Analytical values +μ3, μ4 = μ[3], μ[4] +σ2_1, σ2_2, σ2_3, σ2_4 = 16.0, 4.0, 4e4, 9e4 +σ12, σ34 = 2.4, -1.8e4 +ρ12 = σ12 / (σ1 * σ2) +ρ34 = σ34 / (σ3 * σ4) + +D = σ2_1 * (σ2_3 + μ3^2) + σ2_2 * (σ2_4 + μ4^2) + 2 * σ12 * (σ34 + μ3 * μ4) + +S1_analytical = σ2_1 * (μ3 + μ4 * ρ12 * σ2 / σ1)^2 / D +ST1_analytical = σ2_1 * (1 - ρ12^2) * (σ2_3 + μ3^2) / D + +S2_analytical = σ2_2 * (μ4 + μ3 * ρ12 * σ1 / σ2)^2 / D +ST2_analytical = σ2_2 * (1 - ρ12^2) * (σ2_4 + μ4^2) / D + +S3_analytical = 0.0 +ST3_analytical = σ2_1 * σ2_3 * (1 - ρ34^2) / D + +S4_analytical = 0.0 +ST4_analytical = σ2_2 * σ2_4 * (1 - ρ34^2) / D + + +try + indices, bin_samples = kucherenkoindices_bin([model], inputs, [:y], sim; min_bin_sample_multi_dims=5) println(indices) - # Analytical calculations - denom = 2 + σ^2 + 2*ρ*σ - S1_analytical = 1 / denom - S2_analytical = (1 + ρ*σ)^2 / denom - S3_analytical = (σ + ρ)^2 / denom - ST1_analytical = 1 / denom - ST2_analytical = (1 - ρ^2) / denom - ST3_analytical = σ^2 * (1 - ρ^2) / denom - - - computed_values = Dict() - for i in 1:size(indices, 1) - var_name = string(indices[i, :Variables]) - computed_values[var_name] = (indices[i, :FirstOrder], indices[i, :TotalEffect]) - end - tol = 0.01 - @assert abs(computed_values["x1"][1] - S1_analytical) < tol "S1 error too large" - @assert abs(computed_values["x2"][1] - S2_analytical) < tol "S2 error too large" - @assert abs(computed_values["x3"][1] - S3_analytical) < tol "S3 error too large" - @assert abs(computed_values["x1"][2] - ST1_analytical) < tol "ST1 error too large" - @assert abs(computed_values["x2"][2] - ST2_analytical) < tol "ST2 error too large" - @assert abs(computed_values["x3"][2] - ST3_analytical) < tol "ST3 error too large" - - println("Kucherenko indices validation passed - all values within tolerance.") - + @assert abs(indices.FirstOrder[1] - S1_analytical) < tol "S1 error too large" + @assert abs(indices.FirstOrder[2] - S2_analytical) < tol "S2 error too large" + @assert abs(indices.FirstOrder[3] - S3_analytical) < tol "S3 error too large" + @assert abs(indices.FirstOrder[4] - S4_analytical) < tol "S4 error too large" + @assert abs(indices.TotalEffect[1] - ST1_analytical) < tol "ST1 error too large" + @assert abs(indices.TotalEffect[2] - ST2_analytical) < tol "ST2 error too large" + @assert abs(indices.TotalEffect[3] - ST3_analytical) < tol "ST3 error too large" + @assert abs(indices.TotalEffect[4] - ST4_analytical) < tol "ST4 error too large" + + println("Portfolio model Kucherenko indices validation passed - all values within tolerance.") + catch e - println("Error computing Kucherenko indices: $e") + println("Error computing Portfolio model Kucherenko indices: $e") rethrow(e) end - -println("\nUsing conditional sampling") - -try - indices = kucherenkoindices([model], inputs, [:y], sim) - +try + indices = indices = kucherenkoindices([model], inputs, [:y], sim) println(indices) - # Analytical calculations - denom = 2 + σ^2 + 2*ρ*σ - S1_analytical = 1 / denom - S2_analytical = (1 + ρ*σ)^2 / denom - S3_analytical = (σ + ρ)^2 / denom - ST1_analytical = 1 / denom - ST2_analytical = (1 - ρ^2) / denom - ST3_analytical = σ^2 * (1 - ρ^2) / denom - - - computed_values = Dict() - for i in 1:size(indices, 1) - var_name = string(indices[i, :Variables]) - computed_values[var_name] = (indices[i, :FirstOrder], indices[i, :TotalEffect]) - end - tol = 0.01 - @assert abs(computed_values["x1"][1] - S1_analytical) < tol "S1 error too large" - @assert abs(computed_values["x2"][1] - S2_analytical) < tol "S2 error too large" - @assert abs(computed_values["x3"][1] - S3_analytical) < tol "S3 error too large" - @assert abs(computed_values["x1"][2] - ST1_analytical) < tol "ST1 error too large" - @assert abs(computed_values["x2"][2] - ST2_analytical) < tol "ST2 error too large" - @assert abs(computed_values["x3"][2] - ST3_analytical) < tol "ST3 error too large" - - println("Kucherenko indices validation passed - all values within tolerance.") - + @assert abs(indices.FirstOrder[1] - S1_analytical) < tol "S1 error too large" + @assert abs(indices.FirstOrder[2] - S2_analytical) < tol "S2 error too large" + @assert abs(indices.FirstOrder[3] - S3_analytical) < tol "S3 error too large" + @assert abs(indices.FirstOrder[4] - S4_analytical) < tol "S4 error too large" + @assert abs(indices.TotalEffect[1] - ST1_analytical) < tol "ST1 error too large" + @assert abs(indices.TotalEffect[2] - ST2_analytical) < tol "ST2 error too large" + @assert abs(indices.TotalEffect[3] - ST3_analytical) < tol "ST3 error too large" + @assert abs(indices.TotalEffect[4] - ST4_analytical) < tol "ST4 error too large" + + println("Portfolio model Kucherenko indices validation passed - all values within tolerance.") catch e - println("Error computing Kucherenko indices: $e") + println("Error computing Portfolio model Kucherenko indices: $e") rethrow(e) -end +end \ No newline at end of file From e651f37a17c06226b4fa78e57004dcd4917b2838 Mon Sep 17 00:00:00 2001 From: sitoryu Date: Mon, 25 Aug 2025 16:02:07 +0200 Subject: [PATCH 11/15] added tests for kucherenko inidces and conditional sampling --- test/inputs/copulas/gaussian.jl | 20 ++++ test/runtests.jl | 2 + test/sensitivity/kucherenkoindices.jl | 144 ++++++++++++++++++++++++++ 3 files changed, 166 insertions(+) create mode 100644 test/sensitivity/kucherenkoindices.jl diff --git a/test/inputs/copulas/gaussian.jl b/test/inputs/copulas/gaussian.jl index 0dd0b358a..fb735690a 100644 --- a/test/inputs/copulas/gaussian.jl +++ b/test/inputs/copulas/gaussian.jl @@ -13,6 +13,26 @@ correlation = [1 0.7071; 0.7071 1] @test corkendall(u[:, 1], u[:, 2]) ≈ 0.5 atol = 0.01 end + @testset "sample_conditional_copula" begin + marginals = [RandomVariable(Normal(0,1), :x1), RandomVariable(Normal(0,1), :x2), RandomVariable(Normal(0,1), :x3)] + joint_dist = JointDistribution(marginals, GaussianCopula([1 0.7071 0.5; 0.7071 1 0.3; 0.5 0.3 1])) + unconditional = sample(joint_dist, 500000) + + x_val = 0.5 + cond1 = sample_conditional_copula(joint_dist, [:x1], [x_val], 10000) + filtered1 = unconditional[abs.(unconditional[:, :x1] .- x_val) .< 0.1, :] + @test mean(filtered1[:, :x2]) ≈ mean(cond1[:, :x2]) rtol=0.1 + @test std(filtered1[:, :x2]) ≈ std(cond1[:, :x2]) rtol=0.1 + @test mean(filtered1[:, :x3]) ≈ mean(cond1[:, :x3]) rtol=0.1 + @test std(filtered1[:, :x3]) ≈ std(cond1[:, :x3]) rtol=0.1 + + x_val2, x_val3 = -0.5, 0.2 + cond2 = sample_conditional_copula(joint_dist, [:x2, :x3], [x_val2, x_val3], 10000) + filtered2 = unconditional[(abs.(unconditional[:, :x2] .- x_val2) .< 0.1) .& (abs.(unconditional[:, :x3] .- x_val3) .< 0.1), :] + @test mean(filtered2[:, :x1]) ≈ mean(cond2[:, :x1]) rtol=0.1 + @test std(filtered2[:, :x1]) ≈ std(cond2[:, :x1]) rtol=0.1 + end + @testset "to_standard_normal_space" begin u = sample(copula, 10^5) s = to_standard_normal_space(copula, u) diff --git a/test/runtests.jl b/test/runtests.jl index 00a02afe8..1a7bc002f 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -6,6 +6,7 @@ using InteractiveUtils using QuasiMonteCarlo using Random using StatsBase: fit, Histogram, corkendall +using LinearAlgebra: diag using Test using UncertaintyQuantification @@ -41,6 +42,7 @@ include("reliability/probabilityoffailure_imprecise.jl") include("sensitivity/gradient.jl") include("sensitivity/sobolindices.jl") +include("sensitivity/kucherenkoindices.jl") include("simulations/doe.jl") include("simulations/montecarlo.jl") diff --git a/test/sensitivity/kucherenkoindices.jl b/test/sensitivity/kucherenkoindices.jl new file mode 100644 index 000000000..0e424835b --- /dev/null +++ b/test/sensitivity/kucherenkoindices.jl @@ -0,0 +1,144 @@ +@testset "kucherenkoindices" begin + # Testing Kucherenko indices - Test Case 1: Linear model with correlated variables from Kucherenko et al. (2012) (DOI: 10.1016/j.cpc.2011.12.020) + ρ, σ = 0.5, 2.0 + Σ = [1.0 0.0 0.0; 0.0 1.0 ρ*σ; 0.0 ρ*σ σ^2] + R = Σ ./ (sqrt.(diag(Σ)) * sqrt.(diag(Σ))') + + marginals = RandomVariable[ + RandomVariable(Normal(0, 1), :x1), + RandomVariable(Normal(0, 1), :x2), + RandomVariable(Normal(0, σ), :x3) + ] + + inputs = [ JointDistribution(marginals, GaussianCopula(R)) ] + model = Model(df -> df.x1 .+ df.x2 .+ df.x3, :y) + sim = MonteCarlo(50000) + + # Analytical calculations + denom = 2 + σ^2 + 2*ρ*σ + firstorder_analytical = [1 / denom, (1 + ρ*σ)^2 / denom, (σ + ρ)^2 / denom] + totaleffect_analytical = [1 / denom, (1 - ρ^2) / denom, (σ^2 * (1 - ρ^2)) / denom] + + model_samples = sample(inputs, sim) + evaluate!(model, model_samples) + + + @testset "Standard Kucherenko Indices" begin + indices = kucherenkoindices([model], inputs, [:y], sim) + + @test indices.FirstOrder ≈ firstorder_analytical rtol = 0.1 + @test indices.TotalEffect ≈ totaleffect_analytical rtol = 0.1 + end + + + @testset "Standard Kucherenko Indices - First Order" begin + random_names = names(filter(i -> isa(i, RandomUQInput), inputs)) + Y_orig = Vector(model_samples[:, :y]) + total_var = var(Y_orig) + + S_i = zeros(length(random_names)) + + for (j, i) in enumerate(random_names) + i_cond_samples = UncertaintyQuantification._generate_conditional_samples(model_samples, inputs[1], [i]) + evaluate!(model, i_cond_samples) + S_i[j] = UncertaintyQuantification._compute_first_order_kucherenko(model_samples, i_cond_samples, :y, total_var) + end + @test S_i ≈ firstorder_analytical rtol = 0.1 + end + + @testset "Stnadard Kucherenko Indices - Total Effect" begin + random_names = names(filter(i -> isa(i, RandomUQInput), inputs)) + Y_orig = Vector(model_samples[:, :y]) + total_var = var(Y_orig) + + ST_i = zeros(length(random_names)) + + for (j, i) in enumerate(random_names) + not_i_cond_samples = UncertaintyQuantification._generate_conditional_samples(model_samples, inputs[1], setdiff(random_names, [i])) + evaluate!(model, not_i_cond_samples) + ST_i[j] = UncertaintyQuantification._compute_total_effect_kucherenko(model_samples, not_i_cond_samples, :y, total_var) + end + @test ST_i ≈ totaleffect_analytical rtol = 0.1 + end + + @testset "_generate_conditional_samples" begin + marginals = [RandomVariable(Normal(0,1), :x1), RandomVariable(Normal(0,1), :x2)] + R = [1.0 0.5; 0.5 1.0] + joint_dist = JointDistribution(marginals, GaussianCopula(R)) + samples_df = DataFrame(sample([joint_dist], MonteCarlo(100))) + cond_samples = UncertaintyQuantification._generate_conditional_samples(samples_df, joint_dist, [:x1]) + @test names(cond_samples) == ["x1", "x2"] + @test nrow(cond_samples) == 100 + @test cond_samples[!, :x1] ≈ samples_df[!, :x1] + end + + @testset "Kucherenko Indices with bins" begin + indices, bin_samples = kucherenkoindices_bin([model], inputs, [:y], sim; min_bin_sample_multi_dims=10) + + @test indices.FirstOrder ≈ firstorder_analytical rtol = 0.1 + @test indices.TotalEffect ≈ totaleffect_analytical rtol = 0.1 + end + + @testset "Kucherenko Indices with bins - Existing Samples" begin + random_names = names(filter(i -> isa(i, RandomUQInput), inputs)) + X = Matrix(model_samples[:, random_names]) + Y = Vector(model_samples[:, :y]) + + indices = kucherenkoindices_bin(X, Y; min_bin_sample_multi_dims=10) + + @test indices.FirstOrder ≈ firstorder_analytical rtol = 0.1 + @test indices.TotalEffect ≈ totaleffect_analytical rtol = 0.1 + end + + @testset "Kucherenko Indices with bins - First Order" begin + random_names = names(filter(i -> isa(i, RandomUQInput), inputs)) + X = Matrix(model_samples[:, random_names]) + Y = Vector(model_samples[:, :y]) + n_samples, n_vars = size(X) + total_var = var(Y) + + S_i = zeros(n_vars) + for i in 1:n_vars + S_i[i] = UncertaintyQuantification._compute_first_order_kucherenko_bins(X, Y, i, min(100, floor(Int, n_samples / 25)), total_var) + end + @test S_i ≈ firstorder_analytical rtol = 0.1 + end + + @testset "Kucherenko Indices with bins - Total Order" begin + random_names = names(filter(i -> isa(i, RandomUQInput), inputs)) + X = Matrix(model_samples[:, random_names]) + Y = Vector(model_samples[:, :y]) + n_samples, n_vars = size(X) + total_var = var(Y) + + ST_i = zeros(n_vars) + for i in 1:n_vars + ST_i[i] = UncertaintyQuantification._compute_total_effect_kucherenko_bins(X, Y, i, floor(Int, n_samples / 25), total_var) + end + @test ST_i ≈ totaleffect_analytical rtol = 0.1 + end + + @testset "_assign_multidimensional_bins" begin + X = [1 10 100; + 2 20 110; + 3 30 120; + 4 40 130; + 5 50 140; + 6 60 150; + 7 70 160; + 8 80 170; + 9 90 180; + 10 100 190; + 11 110 200; + 12 120 210] + num_bins = 8 + bin_assignments = UncertaintyQuantification._assign_multidimensional_bins(X, num_bins) + @test length(bin_assignments) == size(X, 1) + @test all(bin_assignments .>= 1) + @test all(bin_assignments .<= num_bins) + @test bin_assignments[1] == bin_assignments[2] + @test bin_assignments[end-1] == bin_assignments[end] + @test bin_assignments[1] != bin_assignments[end] + end + +end \ No newline at end of file From ae5e0d8ab6ce1e0bab2552b57bb97ce80160ce65 Mon Sep 17 00:00:00 2001 From: Laurenz Knipper <60561349+sitoryu@users.noreply.github.com> Date: Wed, 27 Aug 2025 10:58:17 +0200 Subject: [PATCH 12/15] Refactor demo/kucherenkoindices.jl for local package usage Updated local package activation and added print statements/comments for clarity. --- demo/sensitivity/kucherenkoindices.jl | 11 +++++++---- 1 file changed, 7 insertions(+), 4 deletions(-) diff --git a/demo/sensitivity/kucherenkoindices.jl b/demo/sensitivity/kucherenkoindices.jl index 5f7243d8c..b0c19c408 100644 --- a/demo/sensitivity/kucherenkoindices.jl +++ b/demo/sensitivity/kucherenkoindices.jl @@ -1,8 +1,9 @@ -include("../../src/UncertaintyQuantification.jl") -using .UncertaintyQuantification +using Pkg +Pkg.activate(".") +using UncertaintyQuantification using LinearAlgebra -# --- Test Case 2: Portfolio model with analytical indices --- +# Kucherenko et al. (2012) (DOI: 10.1016/j.cpc.2011.12.020) --- Test Case 2: Portfolio model with analytical indices --- μ = [0.0, 0.0, 250.0, 400.0] @@ -62,6 +63,7 @@ ST4_analytical = σ2_2 * σ2_4 * (1 - ρ34^2) / D try indices, bin_samples = kucherenkoindices_bin([model], inputs, [:y], sim; min_bin_sample_multi_dims=5) + println("Sample-based Kucherenko Indices calculation using bins:") println(indices) tol = 0.01 @@ -83,6 +85,7 @@ end try indices = indices = kucherenkoindices([model], inputs, [:y], sim) + println("Standard Kucherenko Indices calculation:") println(indices) tol = 0.01 @@ -99,4 +102,4 @@ try catch e println("Error computing Portfolio model Kucherenko indices: $e") rethrow(e) -end \ No newline at end of file +end From 9ec08e94d90d1f0543bae97707e1d566752cc7bf Mon Sep 17 00:00:00 2001 From: sitoryu Date: Wed, 27 Aug 2025 20:27:03 +0200 Subject: [PATCH 13/15] changed jointdistribution calls to the new constructor/changed sample_conditional_copula var values intput to tuples --- demo/sensitivity/kucherenkoindices.jl | 4 +--- src/inputs/jointdistribution.jl | 19 ++++++++----------- src/sensitivity/kucherenkoindices.jl | 5 +++-- test/inputs/copulas/gaussian.jl | 6 +++--- test/sensitivity/kucherenkoindices.jl | 4 ++-- 5 files changed, 17 insertions(+), 21 deletions(-) diff --git a/demo/sensitivity/kucherenkoindices.jl b/demo/sensitivity/kucherenkoindices.jl index b0c19c408..6b1dd28f8 100644 --- a/demo/sensitivity/kucherenkoindices.jl +++ b/demo/sensitivity/kucherenkoindices.jl @@ -1,5 +1,3 @@ -using Pkg -Pkg.activate(".") using UncertaintyQuantification using LinearAlgebra @@ -31,7 +29,7 @@ Dvec = sqrt.(diag(Σ)) R = Σ ./ (Dvec * Dvec') inputs = [ - JointDistribution(marginals, GaussianCopula(R)) + JointDistribution(GaussianCopula(R), marginals) ] model = Model(df -> df.x1 .* df.x3 .+ df.x2 .* df.x4, :y) diff --git a/src/inputs/jointdistribution.jl b/src/inputs/jointdistribution.jl index 21a2b26e4..9719bf7ac 100644 --- a/src/inputs/jointdistribution.jl +++ b/src/inputs/jointdistribution.jl @@ -43,7 +43,7 @@ struct JointDistribution{ function JointDistribution(d::Copula, m::Vector{<:RandomVariable}) dimensions(d) == length(m) || throw(ArgumentError("Dimension mismatch between copula and marginals.")) - return new{Copula,RandomVariable}(d, m) + return new{typeof(d),RandomVariable}(d, m) end # MultivariateDistribution + Symbol @@ -70,15 +70,12 @@ function sample(jd::JointDistribution{<:MultivariateDistribution,<:Symbol}, n::I return DataFrame(permutedims(rand(jd.d, n)), jd.m) end -function sample_conditional_copula(joint::JointDistribution, var_names::Vector{Symbol}, x_values::Vector{Float64}, N::Int) - length(var_names) != length(x_values) && error("Number of variable names ($(length(var_names))) must match number of values ($(length(x_values)))") - !isa(joint.d, GaussianCopula) && error("This function for conditional sampling is only implemented for Gaussian copulas.") - +function sample_conditional_copula(joint::JointDistribution{<:GaussianCopula,<:RandomVariable}, var_values::Vector{Tuple{Symbol,Float64}}, N::Int) marginals, copula, R = joint.m, joint.d, joint.d.correlation all_var_names, d = [marginal.name for marginal in marginals], length(marginals) v_indices = Int[] - for var_name in var_names + for (var_name, _) in var_values idx = findfirst(==(var_name), all_var_names) idx === nothing && error("Variable $var_name not found in joint distribution. Available variables: $all_var_names") push!(v_indices, idx) @@ -88,8 +85,8 @@ function sample_conditional_copula(joint::JointDistribution, var_names::Vector{S w_indices = setdiff(1:d, v_indices) - z_v = zeros(length(v_indices)) - for (i, (idx, x_val)) in enumerate(zip(v_indices, x_values)) + z_v = zeros(length(var_values)) + for (i, (idx, (_, x_val))) in enumerate(zip(v_indices, var_values)) z_v[i] = quantile(Normal(0,1), cdf(marginals[idx].dist, x_val)) end @@ -108,8 +105,8 @@ function sample_conditional_copula(joint::JointDistribution, var_names::Vector{S Z_w = rand(MvNormal(μ_cond, Symmetric(Σ_cond)), N)' samples = DataFrame() - for (i, var_name) in enumerate(var_names) - samples[!, var_name] = fill(x_values[i], N) + for (var_name, x_val) in var_values + samples[!, var_name] = fill(x_val, N) end for (i, idx) in enumerate(w_indices) samples[!, all_var_names[idx]] = quantile.(marginals[idx].dist, cdf.(Normal(), Z_w[:,i])) @@ -119,7 +116,7 @@ function sample_conditional_copula(joint::JointDistribution, var_names::Vector{S end function sample_conditional_copula(joint::JointDistribution, var_name::Symbol, x_v::Float64, N::Int) - return sample_conditional_copula(joint, [var_name], [x_v], N) + return sample_conditional_copula(joint, [(var_name, x_v)], N) end function to_physical_space!(jd::JointDistribution{<:Copula,<:RandomVariable}, x::DataFrame) diff --git a/src/sensitivity/kucherenkoindices.jl b/src/sensitivity/kucherenkoindices.jl index ac23b4e7e..8e90f9d43 100644 --- a/src/sensitivity/kucherenkoindices.jl +++ b/src/sensitivity/kucherenkoindices.jl @@ -76,10 +76,11 @@ function _generate_conditional_samples( joint_dist::JointDistribution, var_names::Vector{Symbol} ) - input_var_names = [marginal.name for marginal in joint_dist.marginals] + input_var_names = [marginal.name for marginal in joint_dist.m] conditional_samples = map(1:nrow(samples)) do i x_values = [samples[i, var_name] for var_name in var_names] - sample_conditional_copula(joint_dist, var_names, x_values, 1) + var_value_tuples = [(var_names[j], x_values[j]) for j in 1:length(var_names)] + sample_conditional_copula(joint_dist, var_value_tuples, 1) end result = vcat(conditional_samples...) diff --git a/test/inputs/copulas/gaussian.jl b/test/inputs/copulas/gaussian.jl index fb735690a..6daf53fda 100644 --- a/test/inputs/copulas/gaussian.jl +++ b/test/inputs/copulas/gaussian.jl @@ -15,11 +15,11 @@ correlation = [1 0.7071; 0.7071 1] @testset "sample_conditional_copula" begin marginals = [RandomVariable(Normal(0,1), :x1), RandomVariable(Normal(0,1), :x2), RandomVariable(Normal(0,1), :x3)] - joint_dist = JointDistribution(marginals, GaussianCopula([1 0.7071 0.5; 0.7071 1 0.3; 0.5 0.3 1])) + joint_dist = JointDistribution(GaussianCopula([1 0.7071 0.5; 0.7071 1 0.3; 0.5 0.3 1]), marginals) unconditional = sample(joint_dist, 500000) x_val = 0.5 - cond1 = sample_conditional_copula(joint_dist, [:x1], [x_val], 10000) + cond1 = sample_conditional_copula(joint_dist, [(:x1, x_val)], 10000) filtered1 = unconditional[abs.(unconditional[:, :x1] .- x_val) .< 0.1, :] @test mean(filtered1[:, :x2]) ≈ mean(cond1[:, :x2]) rtol=0.1 @test std(filtered1[:, :x2]) ≈ std(cond1[:, :x2]) rtol=0.1 @@ -27,7 +27,7 @@ correlation = [1 0.7071; 0.7071 1] @test std(filtered1[:, :x3]) ≈ std(cond1[:, :x3]) rtol=0.1 x_val2, x_val3 = -0.5, 0.2 - cond2 = sample_conditional_copula(joint_dist, [:x2, :x3], [x_val2, x_val3], 10000) + cond2 = sample_conditional_copula(joint_dist, [(:x2, x_val2), (:x3, x_val3)], 10000) filtered2 = unconditional[(abs.(unconditional[:, :x2] .- x_val2) .< 0.1) .& (abs.(unconditional[:, :x3] .- x_val3) .< 0.1), :] @test mean(filtered2[:, :x1]) ≈ mean(cond2[:, :x1]) rtol=0.1 @test std(filtered2[:, :x1]) ≈ std(cond2[:, :x1]) rtol=0.1 diff --git a/test/sensitivity/kucherenkoindices.jl b/test/sensitivity/kucherenkoindices.jl index 0e424835b..a7a0c691c 100644 --- a/test/sensitivity/kucherenkoindices.jl +++ b/test/sensitivity/kucherenkoindices.jl @@ -10,7 +10,7 @@ RandomVariable(Normal(0, σ), :x3) ] - inputs = [ JointDistribution(marginals, GaussianCopula(R)) ] + inputs = [ JointDistribution(GaussianCopula(R), marginals) ] model = Model(df -> df.x1 .+ df.x2 .+ df.x3, :y) sim = MonteCarlo(50000) @@ -64,7 +64,7 @@ @testset "_generate_conditional_samples" begin marginals = [RandomVariable(Normal(0,1), :x1), RandomVariable(Normal(0,1), :x2)] R = [1.0 0.5; 0.5 1.0] - joint_dist = JointDistribution(marginals, GaussianCopula(R)) + joint_dist = JointDistribution(GaussianCopula(R), marginals) samples_df = DataFrame(sample([joint_dist], MonteCarlo(100))) cond_samples = UncertaintyQuantification._generate_conditional_samples(samples_df, joint_dist, [:x1]) @test names(cond_samples) == ["x1", "x2"] From 692947748dff263a155990bd1fd9bd5c7af33fc8 Mon Sep 17 00:00:00 2001 From: sitoryu Date: Thu, 28 Aug 2025 14:29:20 +0200 Subject: [PATCH 14/15] Added LinearAlgebra to test/Project.toml --- test/Project.toml | 1 + 1 file changed, 1 insertion(+) diff --git a/test/Project.toml b/test/Project.toml index f36f601fe..0e4fef12b 100644 --- a/test/Project.toml +++ b/test/Project.toml @@ -4,6 +4,7 @@ Distributed = "8ba89e20-285c-5b6f-9357-94700520ee1b" HCubature = "19dc6840-f33b-545b-b366-655c7e3ffd49" HypothesisTests = "09f84164-cd44-5f33-b23f-e6b0d136a0d5" InteractiveUtils = "b77e0a4c-d291-57a0-90e8-8db25a27a240" +LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" QuasiMonteCarlo = "8a4e6c94-4038-4cdc-81c3-7e6ffdb2a71b" Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" Revise = "295af30f-e4ad-537b-8983-00126c2a3abe" From 6734b30b6f17a44d69490921447d0b877ea7e81f Mon Sep 17 00:00:00 2001 From: sitoryu Date: Thu, 28 Aug 2025 14:53:42 +0200 Subject: [PATCH 15/15] slightly increased tolerance on 2 dimensional binning test due to low samplesize --- test/inputs/copulas/gaussian.jl | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/test/inputs/copulas/gaussian.jl b/test/inputs/copulas/gaussian.jl index 6daf53fda..803e80695 100644 --- a/test/inputs/copulas/gaussian.jl +++ b/test/inputs/copulas/gaussian.jl @@ -28,9 +28,9 @@ correlation = [1 0.7071; 0.7071 1] x_val2, x_val3 = -0.5, 0.2 cond2 = sample_conditional_copula(joint_dist, [(:x2, x_val2), (:x3, x_val3)], 10000) - filtered2 = unconditional[(abs.(unconditional[:, :x2] .- x_val2) .< 0.1) .& (abs.(unconditional[:, :x3] .- x_val3) .< 0.1), :] - @test mean(filtered2[:, :x1]) ≈ mean(cond2[:, :x1]) rtol=0.1 - @test std(filtered2[:, :x1]) ≈ std(cond2[:, :x1]) rtol=0.1 + filtered2 = unconditional[(abs.(unconditional[:, :x2] .- x_val2) .< 0.15) .& (abs.(unconditional[:, :x3] .- x_val3) .< 0.15), :] + @test mean(filtered2[:, :x1]) ≈ mean(cond2[:, :x1]) rtol=0.15 + @test std(filtered2[:, :x1]) ≈ std(cond2[:, :x1]) rtol=0.15 end @testset "to_standard_normal_space" begin