Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

GSA - Sobol - Question #325

Open
jbensabat opened this issue Jan 1, 2025 · 7 comments
Open

GSA - Sobol - Question #325

jbensabat opened this issue Jan 1, 2025 · 7 comments

Comments

@jbensabat
Copy link

Hello
I am running GSA on a surface flow model (WRF_HYDRO). using the Sobol option, I obtain sensitivity indices that can be negative and greater than 1 (in absolute value).
Since the sensitivity index is the ratio between two variances (a partial and a total) this index should be less than 1 and greater than 0.
How can I check what I am doing wrong ?
best
jac

@jtwhite79
Copy link
Collaborator

Not sure how you got those implausible values. Can you post or send me your control file, rec file, and rmr file?

@jbensabat
Copy link
Author

Good morning and thanks for your reply
I have attached the rec and the control file. I am running pest in serial mode so I have no rmr file
best
jac
HMU17_GSA.Sobol.zip

@jbensabat
Copy link
Author

Hello
I dug a little more into the sobol class in GSA
the sensitivity index s_i is calculated by calling the function
which according to the body of the function can return positive or negative values

double si_saltelli_numer(const std::vector& y_a, const std::vector& y_b, const std::vector& y_c, double missing_val)
{
assert (y_a.size() == y_b.size());
assert (y_a.size() == y_c.size());
int valid_count = 0;
double a, b, c;
double t = 0;
for (int i = 0; i < y_a.size(); i++)
{
a = y_a[i], b = y_b[i], c = y_c[i];
if ((a == missing_val) ||
(b == missing_val) ||
(c == missing_val))
continue;
valid_count++;
t = t + (b * (c - a));

}
t = t / valid_count;

return t;

}
could you please explain ?
thanks
jac

@jtwhite79
Copy link
Collaborator

that formula is from Saltelli's book and its described here: https://en.wikipedia.org/wiki/Variance-based_sensitivity_analysis
under "calculation of indices".

I poked around a little today and didnt see anything obvious (and we have some simple tests that compare results with SALib...). But I need to check what's happening with tied parameters and also whats happening with log transform...

Do you know if you had any runs fail when you ran this originally?

@jbensabat
Copy link
Author

Hi
I did some research outside the scope of pestpp.
negative indices do occur (many do complain). This is attributed numerical errors resulting from either undersampling or the parameter behavior (need to transform the paramater and or apply an affine transformation).
i tried to run with n=4 runs per parameters (4*(np+2)) runs and n=10 runs per parameter (10*(np+2)) run. np is the number of parameters for sensitivity analysis. The obtained sensitivities are different in each run but the number of negative indices increased.
I check - I have no failed runs.
It seems that the problem is far more subtle.
best
jac

@jbensabat
Copy link
Author

Hi
I have further investigated the issue and this is what I have found. Sobol requires hundreds of values per parameter, which is infeasible for the problem I am analyzing (which has ~150 parameters). When working with much less parameters (~10) it is no surprise that we can get negative values of Si (this depends on the parameter values that are picked from the hypercube),
I found in a paper by Saltelli et al a formulation of the indices as suggested by Jansen (1999) which appears to be more robust as it ensures that Si is between 0 and 1
I have attached the paper.
please let me know what you think. I could attempt to add this option to pespp-gsa.
best
jac

PUBLISHED_PAPER.pdf

jansen1999.pdf

@jtwhite79
Copy link
Collaborator

Thanks for digging into this. I'll take a peak at those refs and see what it might take to implement something more clever...

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

No branches or pull requests

2 participants