SciML / GlobalSensitivity.jl

Robust, Fast, and Parallel Global Sensitivity Analysis (GSA) in Julia
https://docs.sciml.ai/GlobalSensitivity/stable/
MIT License
52 stars 20 forks source link

Bootstrap Fails with Homma1996 estimator - Potential fix proposed #128

Closed H-Sax closed 1 year ago

H-Sax commented 1 year ago

When using the Homma1996 estimator to compute Sobol indices and trying to bootstrap the following error occurs.

using QuasiMonteCarlo, GlobalSensitivity
samples = 600000
lb = -ones(4)*π
ub = ones(4)*π
sampler = SobolSample()
A,B = QuasiMonteCarlo.generate_design_matrices(samples,lb,ub,sampler)

function ishi_batch(X)
    A= 7
    B= 0.1
    @. sin(X[1,:]) + A*sin(X[2,:])^2+ B*X[3,:]^4 *sin(X[1,:])
end
res2_HO = gsa(ishi_batch,Sobol(nboot = 1000),A,B,batch=true,Ei_estimator = :Homma1996)
ERROR: BoundsError: attempt to access 2-element Vector{Float64} at index [7]
Stacktrace:
 [1] getindex
   @ ./essentials.jl:13 [inlined]
 [2] (::GlobalSensitivity.var"#37#58"{Int64, Vector{Float64}, Vector{Float64}, Int64, Vector{Vector{Float64}}, Vector{Float64}})(k::Int64)
   @ GlobalSensitivity ./none:0
 [3] iterate
   @ ./generator.jl:47 [inlined]
 [4] collect(itr::Base.Generator{UnitRange{Int64}, GlobalSensitivity.var"#37#58"{Int64, Vector{Float64}, Vector{Float64}, Int64, Vector{Vector{Float64}}, Vector{Float64}}})
   @ Base ./array.jl:782
 [5] gsa_sobol_all_y_analysis(method::Sobol, all_y::Vector{Float64}, d::Int64, n::Int64, Ei_estimator::Symbol, y_size::Nothing, #unused#::Val{false})
   @ GlobalSensitivity ~/.julia/packages/GlobalSensitivity/SJGNh/src/sobol_sensitivity.jl:190
 [6] gsa(f::typeof(ishi_batch), method::Sobol, A::Matrix{Float64}, B::Matrix{Float64}; batch::Bool, Ei_estimator::Symbol, distributed::Val{false}, kwargs::Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
   @ GlobalSensitivity ~/.julia/packages/GlobalSensitivity/SJGNh/src/sobol_sensitivity.jl:139

I have tracked the error down to the definition of the estimator in sobol_sensitivity.jl

            if Ei_estimator === :Homma1996
                push!(Eᵢs,
                    [Varys[i] .- sum(fA .* fAⁱ[k]) ./ (n) + Eys[i] .^ 2 for k in 1:d])

Given I replace Varys[i] and Eys[i] for the estimators as defined in the original paper bootstrapping then computes as expected. Just wanted someones opinon on this before I make a PR. Fix below

            if Ei_estimator === :Homma1996
                push!(Eᵢs, [sum((fA .- (sum(fA) ./ n)).^2) ./ (n-1) .- sum(fA .* fAⁱ[k]) ./ (n) + (sum(fA) ./ n) .^ 2 for k in 1:d])
ChrisRackauckas commented 1 year ago

That looks reasonable to me.