/
workflow_fromVBSAtoEET_hymod.m
124 lines (110 loc) · 5.83 KB
/
workflow_fromVBSAtoEET_hymod.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
% This script shows how the samples produced during VBSA can be re-used to
% estimate EET indices.
%
% METHOD
%
% In this script, first VBSA is applied using the 'resampling' strategy
% for approximating VBSA first-order and total-order indices
% (e.g. Saltelli et al. (2010); more references in 'vbsa_indices.m' help).
% Then, the same input/output samples are used to compute the EET indices
% according to a 'radial' design strategy
% (Campolongo et al. (2011); more in 'EET_indices.m' help).
%
% Because EET is a more frugal method, using the same sample of a rather
% limited size might produce very uncertain VBSA indices
% (i.e. with very large confidence intervals) and
% less uncertain EET indices (i.e. with smaller confidence intervals,
% or at least confidence intervals that do not overlap, so that the input
% ranking can be clearly inferred).
%
% If the user already has a dataset of input-output samples generated by
% the SAFE 'vbsa_resampling' function, then skip Step 1 of this
% workflow and directly start from Step 2.
%
% More on the relation between EET and VBSA can be found in
% (Herman et al, 2013), Campolong et al (2011), Sobol and Kucherenko (2009).
%
% MODEL AND STUDY AREA
%
% The model under study is the rainfall-runoff model Hymod
% (see help of function hymod_sim.m for more details)
% applied to the Leaf catchment in Mississipi, USA
% (Sorooshian et al., 1983).
% The inputs subject to SA are the 5 model parameters, and the scalar
% output for SA is a metric of model performance.
%
% REFERENCES
%
% Campolongo F., Saltelli, A. and J. Cariboni (2011), From screening to
% quantitative sensitivity analysis. A unified approach, Computer Physics
% Communications, 182(4), 978-988.
%
% Herman, J. D., Kollat, J. B., Reed, P. M., and Wagener, T., 2013,
% Technical Note: Method of Morris effectively reduces the computational
% demands of global sensitivity analysis for distributed watershed models,
% Hydrol. Earth Syst. Sci., 17, 2893-2903, doi:10.5194/hess-17-2893-2013.
%
% Saltelli et al. (2010), Variance based sensitivity analysis of model
% output. Design and estimator for the total sensitivity index, Computer
% Physics Communications, 181, 259-270.
%
% I.M. Sobol, S. Kucherenko, 2009, Derivative based global sensitivity measures
% and their link with global sensitivity indices, Mathematics and Computers
% in Simulation, 79(10), 3009-3017.
% http://dx.doi.org/10.1016/j.matcom.2009.01.023.
%
% Sorooshian, S., Gupta, V., Fulton, J. (1983). Evaluation of maximum
% likelihood parameter estimation techniques for conceptual rainfall-runoff
% models: Influence of calibration data variability and length on model
% credibility. Water Resour. Res., 19, 251-259.
% This script prepared by Francesca Pianosi, University of Bristol
% and Corrado Chisari, Universita' di Salerno
% mail to: francesca.pianosi@bristol.ac.uk
%%%% STEP 1: Generate input/output dataset for VBSA %%%%
% Set-up GSA:
load -ascii LeafCatch.txt
rain = LeafCatch(1:365,1) ;
evap = LeafCatch(1:365,2) ;
flow = LeafCatch(1:365,3) ;
M = 5 ; % number of uncertain parameters [ Sm beta alfa Rs Rf ]
% Parameter distribution:
DistrFun = 'unif' ;
% Parameter ranges:
xmin = [ 0 0 0 0 0.1 ] ; % minimum values
xmax = [ 400 2 1 0.1 1 ] ; % maximum values
% Output function:
myfun = 'hymod_nse' ;
% Sample input space according to 'resampling' strategy for approximating
% VBSA first-order and total-order indices:
SampStrategy = 'lhs' ;
N = 100 ; % Base sample size.
% COMMENT: keep this reasonably low (i.e. below 100) so that you should get
% VBSA indices with very large confidence intervals
% (while in Step 2 you should get EET indices with less uncertainty)
for i=1:M; DistrPar{i} = [xmin(i) xmax(i) ] ; end
X = AAT_sampling(SampStrategy,M,DistrFun,DistrPar,2*N);
[ XA, XB, XC ] = vbsa_resampling(X) ; % XA, XB = (N,M); XC = (N*M,M)
% Run the model:
YA = model_evaluation(myfun,XA,rain,evap,flow) ; % size (N,1)
YB = model_evaluation(myfun,XB,rain,evap,flow) ; % size (N,1)
YC = model_evaluation(myfun,XC,rain,evap,flow) ; % size (N*M,1)
% Compute VBSA first-order and total-order indices
% (and associated confidence intervals):
Nboot = 30 ;
[ Si, STi, Si_sd, STi_sd, Si_lb, STi_lb, Si_ub, STi_ub ] = vbsa_indices(YA,YB,YC,Nboot);
figure % plot indices and their confidence bounds:
X_Labels = {'Sm','beta','alfa','Rs','Rf'} ;
boxplot2([Si; STi],X_Labels,[ Si_lb; STi_lb ],[ Si_ub; STi_ub ])
legend('main effects','total effects')
%%%% STEP 2: Apply EET to the dataset generated during VBSA %%%%
% rearrange input and output samples to the format required by EET
% functions:
[ X_EET_radial, Y_EET_radial] = fromVBSAtoEET(XB,YB,XC,YC) ;
% Compute EET indices (and associated confidence intervals):
[mi,sigma,EE,mi_sd,sigma_sd,mi_lb,sigma_lb,mi_ub,sigma_ub] = ...
EET_indices(r,xmin,xmax,X_EET_radial,Y_EET_radial,'radial',Nboot);
% COMMENT: note that EET_indices can be used only with the option
% of 'radial' design in this case, because of the way that XC is obtained
% from XA and XB.
% Plot bootstrapping results in the plane (mean(EE),std(EE)):
EET_plot(mi,sigma,X_Labels,mi_lb,mi_ub,sigma_lb,sigma_ub)