Use randomized sobol sampling

Project.toml

@@ -31,7 +31,7 @@ KernelDensity = "0.6.4"
 LinearAlgebra = "1.10"
 OrdinaryDiffEq = "6.62"
 Parameters = "0.12"
-QuasiMonteCarlo = "0.2.3, 0.3"
+QuasiMonteCarlo = "0.3.3"
 Random = "1.10"
 RecursiveArrayTools = "3.2"
 SafeTestsets = "0.1"

docs/src/tutorials/shapley.md

@@ -128,8 +128,8 @@ barplot(
         xticklabelrotation = 1,
         xticks = (1:54, ["θ$i" for i in 1:54]),
         ylabel = "Shapley Indices",
-        limits = (nothing, (0.0, 0.2)),
-    ),
+        limits = (nothing, (0.0, 0.2))
+    )
 )
 ```
 
@@ -160,7 +160,7 @@ barplot(
         xticklabelrotation = 1,
         xticks = (1:54, ["θ$i" for i in 1:54]),
         ylabel = "Shapley Indices",
-        limits = (nothing, (0.0, 0.2)),
-    ),
+        limits = (nothing, (0.0, 0.2))
+    )
 )
 ```

src/shapley_sensitivity.jl

@@ -203,6 +203,11 @@ function gsa(f, method::Shapley, input_distribution::SklarDist; batch = false)
             sample_complement = rand(
                 Copulas.subsetdims(input_distribution, idx_minus), n_outer)
 
+            if size(sample_complement, 2) == 1
+                sample_complement = reshape(
+                    sample_complement, (1, length(sample_complement)))
+            end
+
             for l in 1:n_outer
                 curr_sample = @view sample_complement[:, l]
                 # Sampling of the set conditionally to the complementary element

src/sobol_sensitivity.jl

@@ -342,10 +342,18 @@ function gsa_sobol_all_y_analysis(method, all_y::AbstractArray{T}, d, n, Ei_esti
         nboot > 1 ? reshape(ST_CI, size_...) : nothing)
 end
 
-function gsa(f, method::Sobol, p_range::AbstractVector; samples, kwargs...)
+function gsa(f, method::Sobol, p_range::AbstractVector; samples,
+        rng::AbstractRNG = Random.default_rng(), kwargs...)
+    log2num = ceil(Int, log2(samples))
+    samples2n = 2^log2num
+    if samples2n != samples
+        samples = samples2n
+        @warn "Passed samples is not a power of 2, number of sample points changed to $samples"
+    end
     AB = QuasiMonteCarlo.generate_design_matrices(samples, [i[1] for i in p_range],
         [i[2] for i in p_range],
-        QuasiMonteCarlo.SobolSample(),
+        QuasiMonteCarlo.SobolSample(;
+            R = QuasiMonteCarlo.OwenScramble(; base = 2, pad = log2num, rng)),
         2 * method.nboot)
     A = reduce(hcat, @view(AB[1:(method.nboot)]))
     B = reduce(hcat, @view(AB[(method.nboot + 1):end]))

test/shapley_method.jl

@@ -25,8 +25,8 @@ n_perms = -1;
 n_var = 10_000;
 n_outer = 1000;
 n_inner = 3;
-dim = 3;
-margins = (Uniform(-pi, pi), Uniform(-pi, pi), Uniform(-pi, pi));
+dim = 4;
+margins = (Uniform(-pi, pi), Uniform(-pi, pi), Uniform(-pi, pi), Uniform(-pi, pi));
 dependency_matrix = Matrix(4 * I, dim, dim);
 C = GaussianCopula(dependency_matrix);
 input_distribution = SklarDist(C, margins);
@@ -41,17 +41,17 @@ method = Shapley(n_perms = n_perms,
 
 @test result.shapley_effects[1]≈0.43813841765976547 atol=1e-1
 @test result.shapley_effects[2]≈0.44673952698721386 atol=1e-1
-@test result.shapley_effects[3]≈0.23144736934254417 atol=1e-1
-# @test result.shapley_effects[4]≈0.0 atol=1e-1
+@test result.shapley_effects[3]≈0.11855122481995543 atol=1e-1
+@test result.shapley_effects[4]≈0.0 atol=1e-1
 #<---- non batch
 
 #---> batch
 result = gsa(ishi_batch, method, input_distribution, batch = true);
 
 @test result.shapley_effects[1]≈0.44080027198796035 atol=1e-1
 @test result.shapley_effects[2]≈0.43029987176805085 atol=1e-1
-@test result.shapley_effects[3]≈0.23144736934254417 atol=1e-1
-# @test result.shapley_effects[4]≈0.0 atol=1e-1
+@test result.shapley_effects[3]≈0.11855122481995543 atol=1e-1
+@test result.shapley_effects[4]≈0.0 atol=1e-1
 #<--- batch
 
 d = 3