Feature/kucherenko indices

demo/sensitivity/kucherenkoindices.jl

@@ -0,0 +1,103 @@
+using UncertaintyQuantification
+using LinearAlgebra
+
+# 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]
+
+σ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
+
+Σ = [
+    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(μ[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(GaussianCopula(R), marginals)
+]
+
+model = Model(df -> df.x1 .* df.x3 .+ df.x2 .* df.x4, :y)
+sim = MonteCarlo(200000)
+
+
+# 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("Sample-based Kucherenko Indices calculation using bins:")
+    println(indices)
+
+    tol = 0.01
+    @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 Portfolio model Kucherenko indices: $e")
+    rethrow(e)
+end
+
+try 
+    indices = indices = kucherenkoindices([model], inputs, [:y], sim)
+    println("Standard Kucherenko Indices calculation:")
+    println(indices)
+
+    tol = 0.01
+    @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 Portfolio model Kucherenko indices: $e")
+    rethrow(e)
+end

src/UncertaintyQuantification.jl

@@ -167,7 +167,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
@@ -220,6 +223,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")

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,6 +70,55 @@ function sample(jd::JointDistribution{<:MultivariateDistribution,<:Symbol}, n::I
     return DataFrame(permutedims(rand(jd.d, n)), jd.m)
 end
 
+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_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)
+    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(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
+
+    Σ_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 (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]))
+    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{<:Copula,<:RandomVariable}, x::DataFrame)
     correlated_cdf = to_copula_space(jd.d, Matrix{Float64}(x[:, names(jd)]))
     for (i, rv) in enumerate(jd.m)