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/juliacon21.md

@@ -5,7 +5,7 @@ We showcase how to use multiple GSA methods, analyze their results and leverage
 perform Global Sensitivity analysis at scale.
 
 ```@example lv
-using GlobalSensitivity, QuasiMonteCarlo, OrdinaryDiffEq, Statistics, CairoMakie
+using GlobalSensitivity, QuasiMonteCarlo, OrdinaryDiffEq, Statistics, CairoMakie, Random
 
 function f(du, u, p, t)
     du[1] = p[1] * u[1] - p[2] * u[1] * u[2] #prey
@@ -54,7 +54,7 @@ fig
 ```
 
 ```@example lv
-sobol_sens = gsa(f1, Sobol(), bounds, samples = 500)
+sobol_sens = gsa(f1, Sobol(), bounds, samples = 512)
 efast_sens = gsa(f1, eFAST(), bounds, samples = 500)
 ```
 
@@ -94,10 +94,10 @@ fig
 
 ```@example lv
 using QuasiMonteCarlo
-samples = 500
+samples = 512
 lb = [1.0, 1.0, 1.0, 1.0]
 ub = [5.0, 5.0, 5.0, 5.0]
-sampler = SobolSample()
+sampler = SobolSample(; R = QuasiMonteCarlo.OwenScramble(; base = 2, pad = 9, rng = _rng))
 A, B = QuasiMonteCarlo.generate_design_matrices(samples, lb, ub, sampler)
 sobol_sens_desmat = gsa(f1, Sobol(), A, B)
 
@@ -129,7 +129,7 @@ f1 = function (p)
     prob1 = remake(prob; p = p)
     sol = solve(prob1, Tsit5(); saveat = t)
 end
-sobol_sens = gsa(f1, Sobol(nboot = 20), bounds, samples = 500)
+sobol_sens = gsa(f1, Sobol(nboot = 20), bounds, samples = 512)
 fig = Figure(resolution = (600, 400))
 ax, hm = CairoMakie.scatter(
     fig[1, 1], sobol_sens.S1[1][1, 2:end], label = "Prey", markersize = 4)

docs/src/tutorials/parallelized_gsa.md

@@ -65,7 +65,7 @@ scatter(
 For the Sobol method, we can similarly do:
 
 ```@example ode
-m = gsa(f1, Sobol(), [[1, 5], [1, 5], [1, 5], [1, 5]], samples = 1000)
+m = gsa(f1, Sobol(), [[1, 5], [1, 5], [1, 5], [1, 5]], samples = 1024)
 ```
 
 ## Direct Use of Design Matrices
@@ -76,10 +76,10 @@ we use [QuasiMonteCarlo.jl](https://docs.sciml.ai/QuasiMonteCarlo/stable/) to ge
 as follows:
 
 ```@example ode
-samples = 500
+samples = 512
 lb = [1.0, 1.0, 1.0, 1.0]
 ub = [5.0, 5.0, 5.0, 5.0]
-sampler = SobolSample()
+sampler = SobolSample(; R = QuasiMonteCarlo.OwenScramble(; base = 2, pad = 9, rng = _rng))
 A, B = QuasiMonteCarlo.generate_design_matrices(samples, lb, ub, sampler)
 ```
 

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

@@ -46,6 +46,10 @@ by dividing other terms in the variance decomposition by `` Var(Y) ``.
 - `:Jansen1999` - [M.J.W. Jansen, 1999, Analysis of variance designs for model output, Computer Physics Communi- cation, 117, 35–43.](https://www.sciencedirect.com/science/article/abs/pii/S0010465598001544)
 - `:Janon2014` - [Janon, A., Klein, T., Lagnoux, A., Nodet, M., & Prieur, C. (2014). Asymptotic normality and efficiency of two Sobol index estimators. ESAIM: Probability and Statistics, 18, 342-364.](https://arxiv.org/abs/1303.6451)
 
+!!! note
+
+    Sobol sampling should be done with $2^k$ points and randomization, take a look at the docs for [QuasiMonteCarlo](https://docs.sciml.ai/QuasiMonteCarlo/stable/randomization/). If the number of samples is not a power of 2, the number of sample points will be changed to the next power of 2.
+
 ### Example
 
 ```julia
@@ -57,7 +61,7 @@ function ishi(X)
     sin(X[1]) + A*sin(X[2])^2+ B*X[3]^4 *sin(X[1])
 end
 
-samples = 600000
+samples = 524288
 lb = -ones(4)*π
 ub = ones(4)*π
 sampler = SobolSample()
@@ -342,10 +346,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...)
+    samples2n = nextpow(2, samples)
+    if samples2n != samples
+        samples = samples2n
+        @warn "Passed samples is not a power of 2, number of sample points changed to $samples"
+    end
+    log2num = round(Int, log2(samples))
     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

test/sobol_method.jl

@@ -22,7 +22,7 @@ function linear(X)
     A * X[1] + B * X[2]
 end
 
-n = 600000
+n = 524288
 lb = -ones(4) * π
 ub = ones(4) * π
 sampler = SobolSample()
@@ -48,10 +48,7 @@ res1 = gsa(ishi, Sobol(order = [0, 1, 2], nboot = 20), A, B)
                0.0 0.0 4.279200031668718e-5 1.2542212940962112e-5;
                0.0 0.0 0.0 -7.998213172266514e-7; 0.0 0.0 0.0 0.0] atol=1e-4
 @test res1.S1_Conf_Int≈[
-    0.00013100970128286063,
-    0.00014730548523359544,
-    7.398816006175431e-5,
-    0.0
+    0.00018057916212640867, 0.0002327582551601999, 9.116874912775071e-5, 0.0
 ] atol=1e-4
 @test res1.ST_Conf_Int≈[
     5.657364947147881e-5,
@@ -94,7 +91,7 @@ ishigami.fun <- function(X) {
   B <- 0.1
   A * X[, 1] + B * X[, 2]
 }
-n <- 6000000
+n <- 524288
 X1 <- data.frame(matrix(runif(4 * n,-pi,pi), nrow = n))
 X2 <- data.frame(matrix(runif(4 * n,-pi,pi), nrow = n))
 sobol2007(ishigami.fun, X1, X2)