General interactive user guide for GSUA-CSB Toolbox

Global Sensitivity and Uncertainty Analysis - Confidence Subcontour Box (GSUA-CSB) Toolbox is a product developed by Universidad EAFIT for command-line mathematical model validation in both of Simulink or Symbolic Math Toolbox environment . At present time, the toolbox allows to perform the following functions: To apply and visualize several variance-based sensitivity (SA) and uncertainty (UA) analysis, to estimate model parameters (PE) and to estimate confidence subcontour boxes (CSB) for estimated parameters.. This toolbox is based on the previous work of Carlos Mario Vélez: GSUA of dynamical systems using variance-based methods, published in this mathworks file exchange link.

Developers

Daniel Rojas Díaz - drojasd@eafit.edu.co
Carlos Mario Vélez Sánchez - cmvelez@eafit.edu.co

Additional resources

Toolbox paper, Slides on UA/SA, model identifiability paper.

Table of Contents

General interactive user guide for GSUA-CSB Toolbox Developers Additional resources Initial configuration Working with Simulink models Working with Symbolic math models Working with User-defined models Toolbox routines for selected model 1. Simulating the model 1.1 gsua_eval function 2. Sample the space of factors 2.1 gsua_dmatrix function 3. Uncertainty analysis 3.1 gsua_ua function 4. Sensitivity analysis (SA) 4.1 gsua_sa function 4.2 SA results visualization and interpretation 5. Parameter estimation 5.1 gsua_pe function 5.2 Identifiability analysis 6. Uncertainty Interval Calculus 6.1 gsua_oatr function 6.2 gsua_csb function 7. Fit quality, acceptance and confidence intervals 7.1 gsua_covmetric function 7.2 gsua_costcutoff function 7.3 gsua_noisefloor function 7.4 gsua_likelihood function 8. Where the other functions fit 9. Python port Interactive user guide for Toolbox initial configuration (Simulink) 1. Model preparation 2. Environment configuration and additional features 2.1 Additional features: 2.1.1 Range definition 2.1.2 Nominal factor family 2.1.3 Model kind (default: Dynamic) 2.1.4 Fixed step (Default: []) Interactive user guide for Toolbox initial configuration (Symbolic Matlab) 1. Model creation 2. Environment configuration and additional features 2.1 Additional features 2.1.1 Range definition 2.1.2 Nominal factor family 2.1.3 Specific model output (Default: 1:model_order) 2.1.4 Specific ode solver (Default: ode45) 2.1.5 Ode solver configuration (Default: none) 2.1.6 Static model simulation Interactive user guide for Toolbox initial configuration (User-Defined functions) 1. Time-dependent functions 1.1 Model creation 1.1.1 Example 1.2 Model implementation 1.3 Additional features 1.3.1 Range definition 1.3.2 Nominal factor family 1.3.3 Specific model output (Default: 1:model_order) 1.3.4 Specific names for model factors (Default: {}) 1.3.5 Specific names for model outputs (Default: {}) 1.3.6 Run as vectorized function (Default: false) 1.3.7 Any additional user-provided configuration (Default: {}) 2. Time-free functions 2.1 Model creation 1.1.1 Example 2.2 Model implementation 2.3 Additional features

IMPORTANT NOTE: If you lost any of the results in this guide, you can recover it easy with the following command:

load userguidesession.mat

Initial configuration

Working with Simulink models

It is necessary to create and configure a simulink model as indicated in GSUA Simulink userguide. Also, check the guide for additional configurations.

[T,tout]=gsua_dataprep(model_name,Ranges,Parameter_names)

Working with Symbolic math models

It is necessary to create a system of equations as indicated in GSUA Symbolic Math userguide. Also, check the guide for additional configurations.

[T,solver]=gsua_dataprep(odes,vars,domain,model_name)

Working with User-defined models

It is necessary to create a system of equations as indicated in GSUA User-defined functions userguide. Also, check the guide for additional configurations.

[T,sol]=gsua_dataprep(func,ranges,'domain',domain)

Toolbox routines for selected model

For the following examples we will use a symbolic math model (same as constructed in GSUA Symbolic Math userguide). You can execute your own model or the Simulink model example following the respective userguide steps.

1. Simulating the model

%load default example model

default_example

Setting environment to work with symbolic Matlab
Introduce ranges in the following order: 

ans(t) =

T%check the summary table

T = 4×2 table 

	Range	Nominal
	undefined	Nominal
1 S0	100	3000
2 I0	1	500
3 beta	0	1
4 gamma	0	1

1.1 gsua_eval function

help gsua_eval%check features
  Function for few model evaluations
  
  Y=gsua_eval(values,Table)
  Parameters:
  values <-- array of NpxN (number of factors x number of simulations)
  Table  <-- summary table from gsua_dataprep function
  Outputs:
  Y <-- array with model output
  Additional features:
  If you provide an array (xdata) that belongs to the model domain, then
  the model output is interpolated to match xdata. Also you can provide a
  previous model output or experimental data as an array (ydata) that must
  coincide with the xdata array:
  Y=gsua_eval(values,Table,xdata,ydata)

Execute the following line to modify the number of model outputs (valid only for symbolic math and user-defined models)

T.Properties.CustomProperties.output=1:3;%set to visualize the three states of the model

y=gsua_eval(T.Nominal,T);% simple simulation with factor values stored in T.Nominal

%now, we are going to use the additional features of the function

xdata=linspace(0,264,265);%create 265 points between 0 and 264

noise= normrnd(0,0.2,size(xdata));%create noise

%simulate the model with user defined parameters and return values in xdata

y=gsua_eval([1000 200 0.1 0.06]',T,xdata);

yexp=y+y.*noise;%create a noisy output

%compare results and simulate another factors family

y=gsua_eval([1000 200 0.1 0.06; 2000 200 0.1 0.06]',T,xdata,yexp);

2. Sample the space of factors

2.1 gsua_dmatrix function

help gsua_dmatrix

  Function for design of experiments (factor space sampling)
 
  [M,T2]=gsua_dmatrix(T,N)
  Parameters:
  T <-- summary table from gsua_dataprep
  N <-- number of samples
  Outpus:
  M  <-- design matrix of NxNp for later routines
  T2 <-- summary table with fixed parameters actualized
  Additional features:
  You can choose a method for factor space sampling between uniform
  distribution design and latin hypercube design (default). To switch
  between methods use the paired feature 'Method' and 'Uniform',
  'Sobol', or 'LatinHypercube'.
  Also you can visualize the sampling result using the paired feature
  'Show', 'on'
  M=gsua_dmatrix(T,N,'Method','Sobol','Show','on')

M=gsua_dmatrix(T,500,'Show','on');%design matrix with 500 samples

%we ignore T2 because there is no new fixed parameters

The methods above are marginal: each factor is drawn independently inside its own range. When the factors are correlated -- which is what gsua_ia's correlation heatmap shows, and what a cost_band far above cost_data in gsua_covmetric signals -- independent draws leave the identified manifold, and the band built from them is both too wide and badly centred, with a median that resembles no good fit. 'Method','Joint' instead draws whole factor vectors from an ensemble of accepted estimates, so the correlation structure survives sampling.

Tia = gsua_ia(T,T.Estlsqc); % Tia.Est holds the accepted ensemble

M = gsua_dmatrix(Tia,2000,'Method','Joint');

% 'JointType' -> Bootstrap (default), SmoothBootstrap or Gaussian

% 'Pool' supplies the ensemble explicitly instead of reading Tia.Est

Joint draws also inherit the ensemble's own spread, which matters more than it sounds: gsua_ia replaces T.Range with a confidence interval of the MEDIAN of the pool, and that interval narrows as the pool grows. It describes where the centre lies, not the spread of factor values consistent with the data, so sampling it understates uncertainty. For that reason joint draws are not clipped back to T.Range by default.

3. Uncertainty analysis

3.1 gsua_ua function

help gsua_ua

  Function for uncertainty analysis
 
  Y=gsua_ua(M,T)
  Parameters:
  M <-- design matrix from gsua_dmatrix function
  T <-- summary table from gsua_dataprep function
  Output:
  Y <-- result of Monte-Carlo simulation
  Additional features:
  The function automatically apply a Monte-Carlo filtering over results and
  present it as figures (gsua_MCF function). It is possible to perform UA
  over multiple model outputs. To obtain extrapolated results, use the
  paired feature 'xdata'. To compare results with a specific output use
  the paired feature 'ynom' (ynom length must coindice with xdata length).
  To avoid parallel computing (no speed up), use the paired feature 
  'parallel',false.
  Y=gsua_ua(M,T,'xdata',xdata,'ynom',ynom,'parallel',false)

M=gsua_dmatrix(T,500);%design matrix

xdata=linspace(0,264,265);%265 points between 0 and 264

Y=gsua_ua(M,T,'xdata',xdata);

Progress: 100%
Estimated processing time (h:m:s): 0:0:1
Remaining time (h:m:s): 0:0:0
Elapsed time (h:m:s): 0:0:1
Estimated stop time (h:m:s): 2:31:45
Number of simulations: 500

Since we do not specify a nominal output, the function will take as nominal the output associated with T.Nominal values (red line in first figure above).

The Monte-Carlo filtering interpretation is as follows:

The bigger the difference between low values (blue) and high values (red) lines, the bigger is the factor relevance.
If low values (blue) is the bottom line, then high values for the correspondent factor generates low values for model output and low values for the correspondent factor generates high values for model output..
If high values (red) is the bottom line, then high values for the correspondent factor generates high values for model output and low values for the correspondent factor generates low values for model output.

%if you change ynom, then Monte-Carlo filtering results could vary

ynom=gsua_eval([1000 200 0.1 0.06]',T,xdata);

Y=gsua_ua(M,T,'xdata',xdata,'ynom',ynom);

Progress: 100%
Estimated processing time (h:m:s): 0:0:1
Remaining time (h:m:s): 0:0:0
Elapsed time (h:m:s): 0:0:1
Estimated stop time (h:m:s): 2:31:50
Number of simulations: 500

Now, results interpretation is much more tricky. It is better to perform a SA (next section).

4. Sensitivity analysis (SA)

4.1 gsua_sa function

Sensitivity indices estimators implemented in this toolbox are based on the following works:

[1]: Saltelli, A., Annoni, P., Azzini, I., Campolongo, F., Ratto, M., and Tarantola, S. (2010). Variance based sensitivity analysis of model output. design and estimator for the total sensitivity index. Computer Physics Communications, 181(2):259–270.

[2]: Xiao, S., Lu, Z., and Wang, P. (2018). Multivariate global sensitivity analysis based on distance components decomposition. Risk Analysis, 38(12):2703–2721.

Saltelli method has been configured to work with Minimum Least-Squares outputs, if you do not provide a specific 'ynom' the function assume that 'ynom' is given by the nominal Factors in T.

Xiao method is especially good for time-dependent model responses.

help gsua_sa
  Function for sensitivity analysis
 
  T=gsua_sa(M,T)
  Parameters:
  M <-- design matrix from gsua_dmatrix function
  T <-- summary table from gsua_dataprep function
  Output:
  T <-- summary table with sensitivity indices.
  Additional features:
  It is possible to request 2 extra positional outputs:
  J <-- vector of model outputs in scalar representation (MSE)
  Y <-- model output matrix (same as gsua_ua output)
  It is possible to choose between 6 SA methods using the paired feature
  'SensMethod',method. Where method is one of the following:
  'Saltelli'    <-- requires N*(Np/2+1) simulations
  'Jansen'      <-- requires N*(Np/2+1) simulations
  'Xiao'        <-- newest method (Default), requires N*(Np/2+1)
                    simulations
  'Sobol'       <-- classic method, requires N*(Np+1) simulations
  'brute-force' <-- explores all possible combinations, N + Np2*N^2
                    simulations
  'OAT'         <-- non-global method, requires N*(Np+1) simulations
  Also, it is possible to specify an output for indices estimation using
  the paired feature 'ynom'. Finally, you can avoid parallel speed up 
  for simulation process using the paired feature 'parallel', false.
  [ParT,J,Y] = gsua_sa(M,T,'SensMethod','Xiao','ynom',ynom)
  Note: ynom size must coincide with xdata size, where xdata is the next
  xdata=linspace(Domain(1),Domain(2),Domain(2)-Domain(1)+1) and Domain is
  the array of the model domain

Let's use the same M matrix as in previous section, without a specific ynom. Hence, the toolbox will take as ynom the output associated to T.Nominal values. Sensitivity indices are not affected by ynom selection. Additional outputs J and Y will be used for posterior analysis.

[T,J,Y]=gsua_sa(M,T,'SensMethod','Xiao');
Progress: 17%
Estimated processing time (h:m:s): 0:0:5
Remaining time (h:m:s): 0:0:4
Elapsed time (h:m:s): 0:0:0
Estimated stop time (h:m:s): 2:32:11
Number of simulations: 1500
Progress: 34%
Estimated processing time (h:m:s): 0:0:6
Remaining time (h:m:s): 0:0:4
Elapsed time (h:m:s): 0:0:2
Estimated stop time (h:m:s): 2:32:12
Number of simulations: 1500
Progress: 50%
Estimated processing time (h:m:s): 0:0:8
Remaining time (h:m:s): 0:0:4
Elapsed time (h:m:s): 0:0:4
Estimated stop time (h:m:s): 2:32:14
Number of simulations: 1500
Progress: 67%
Estimated processing time (h:m:s): 0:0:9
Remaining time (h:m:s): 0:0:3
Elapsed time (h:m:s): 0:0:6
Estimated stop time (h:m:s): 2:32:15
Number of simulations: 1500
Progress: 84%
Estimated processing time (h:m:s): 0:0:9
Remaining time (h:m:s): 0:0:1
Elapsed time (h:m:s): 0:0:8
Estimated stop time (h:m:s): 2:32:15
Number of simulations: 1500
Progress: 100%
Estimated processing time (h:m:s): 0:0:9
Remaining time (h:m:s): 0:0:0
Elapsed time (h:m:s): 0:0:9
Estimated stop time (h:m:s): 2:32:15
Number of simulations: 1500

4.2 SA results visualization and interpretation

gsua_plot('UncertaintyAnalysis',T,Y,xdata,ynom); <- Uncertainty analysis graph
gsua_plot('ScatterOutput',T,Y,M,xdata,xref); <- Uncertainty analysis for output in time=xref
gsua_plot('ScatterOutput',T,J,M,xdata); <- Uncertainty analysis for scalar output (MSE of the output with respect to ynom)
gsua_plot('Pie',T,T.STi/Si); <- Pie graph of sensitivity indices
gsua_plot('Bar',T,T.STi/Si); <- Bar graph of sensitivity indices
gsua_plot('FractionalSensitivityArea',T,T.STi_vec/Si_vec,xdata); <- Area graph of vectorial sensitivity indices (a index for each point in xdata)
gsua_plot('TotalSensitivityArea',T,T.STi_vec/Si_vec,xdata); <- Normalized area graph of vectorial sensitivity indices (a index for each point in xdata)
gsua_plot('FractionalSensitivityPlots',T,T.STi_vec/Si_vec,xdata); <- Figure of subgraphs with the sensitivity indices plot for each factor, and for each point in xdata.
gsua_plot('Pie',T,T.STi_vec/Si_vec,xdata,[t1,t2,t3,t4]); <- Figure of pie subgraphs. As many subgraphs as points in [t1,t2,...tn]. Each pie represent the sensitivity indices for a specific time [t1,t2,...tn].

The most relevant index to assess parameter relevance is STi. However, it is necessary to point out a relevant feature about the indices: The sum of STi is always greater or equal to one, while the sum of Si is always lesser or equal to one. If one of the sums is equal to one, the other one must be equal too, and hence, the model is purely additive.

A general advice is to perform a global SA over the model when the interaction among factors it is unknown. If the sum of Si is greater than 0.65, then you can state that the model has weak interactions among factors and said interactions are negligible. Thereby, for subsequent analysis, it is better to perform local SA analysis or analytical SA methods.

IMPORTANT NOTE: Choosing an appropriate N for SA is a challenging task. However, Saltelli and Distance method allow an approximation: A good N has been chosen when there are no negative Si. Thereby, if you get negative Si values, you should consider to try a greater N.

T.Si

ans = 4×1

sum(T.Si)

ans = 0.9270

sum(T.STi)

ans = 1.2884

clf

gsua_plot('Pie',T,T.STi)

colormap jet

clf

gsua_plot('TotalSensitivityArea',T,T.STi_vec,xdata);

colormap jet

clf

gsua_plot('ScatterOutput',T,J,M,xdata);

5. Parameter estimation

5.1 gsua_pe function

help gsua_pe

  Parameter estimation function
  [T,res] = gsua_pe(T,xdata,ydata,N)
  Parameters:
  T     <-- summary table from gsua_dataprep function
  xdata <-- array of points where the model will be evaluated
  ydata <-- array with expected model output
  Outputs:
  T   <-- summary table with parameter estimation results
  res <-- cost functions for each estimation
  Additional paired features:
  'N',N     <-- number of parameter estimations
  'Multistart',k <-- activate multistart feature for lsqcurvefit optimizer,
  perform k parameter estimations in each cycle
  'solver',{'lsqc','lsqn','ga','particle','psearch','surrogate','annealing'}
  <-- allows to choose among several matlab optimizers. Default: lsqc
  'opt',optimoptions(optimizer,...) <-- allows to configure the respective
  matlab optimizer, optimizer must match with the real matlab optimizer
  name.
  'ipoint',point <-- allows to give to the optimizer the initial points for
  estimations. point must have the same number of columns as factors to
  estimate and the same number of rows as N.
  'Show',{'off','on'}. If activated, the function plot optimizer results
  Default: off.
  [T,res] = gsua_pe(T,xdata,ydata,'Show','on','solver','particle','N',3)

solver='lsqc';%this is the default solver

opt=optimoptions('lsqcurvefit','UseParallel',true);%configure optimizer

[T,res]=gsua_pe(T,xdata,yexp,'solver',solver,'Show','on','opt',opt,'N',10);

Generating a valid matrix for estimations

Estimation 1

Local minimum possible.

lsqcurvefit stopped because the final change in the sum of squares relative to 
its initial value is less than the value of the function tolerance.

<stopping criteria details>
Estimation 2

Local minimum possible.

lsqcurvefit stopped because the final change in the sum of squares relative to 
its initial value is less than the value of the function tolerance.

<stopping criteria details>
Estimation 3

Local minimum possible.

lsqcurvefit stopped because the final change in the sum of squares relative to 
its initial value is less than the value of the function tolerance.

<stopping criteria details>
Estimation 4

Local minimum possible.

lsqcurvefit stopped because the final change in the sum of squares relative to 
its initial value is less than the value of the function tolerance.

<stopping criteria details>
Estimation 5

Local minimum possible.

lsqcurvefit stopped because the final change in the sum of squares relative to 
its initial value is less than the value of the function tolerance.

<stopping criteria details>
Estimation 6

Local minimum possible.

lsqcurvefit stopped because the final change in the sum of squares relative to 
its initial value is less than the value of the function tolerance.

<stopping criteria details>
Estimation 7

Local minimum possible.

lsqcurvefit stopped because the final change in the sum of squares relative to 
its initial value is less than the value of the function tolerance.

<stopping criteria details>
Estimation 8

Local minimum possible.

lsqcurvefit stopped because the final change in the sum of squares relative to 
its initial value is less than the value of the function tolerance.

<stopping criteria details>
Estimation 9

Local minimum possible.

lsqcurvefit stopped because the final change in the sum of squares relative to 
its initial value is less than the value of the function tolerance.

<stopping criteria details>
Estimation 10

Local minimum possible.

lsqcurvefit stopped because the final change in the sum of squares relative to 
its initial value is less than the value of the function tolerance.

<stopping criteria details>

Beyond the paired features used above, two are worth knowing. 'margin' and 'alpha' switch gsua_pe to the correlation-penalized multi-objective cost (gsua_costf) instead of plain least squares, which is what you want when one model has to match several signals at once; gsua_pe then records both on the table so that gsua_noisefloor (section 7.3) can recover them. 'ipoint' seeds the multistart from a specific factor family rather than at random, and 'A', 'B', 'Aeq', 'Beq' and 'nonlcon' pass linear and nonlinear constraints through to the solver.

T%look the actual table, there is a lot of information
T = 4×7 table 
RangeNominalSiSTi
undefinedundefinedundefinedundefined
1 S0100300015500.16990.133900.0232
2 I01500250.50000.02200.02140.92160.6237
3 beta010.50000.04560.063500.1204
4 gamma010.50000.68951.069600.2518

	Range	Nominal	Si	STi
1 S0	100	3000	1550	0.1699	0.1339	0	0.0232
2 I0	1	500	250.5000	0.0220	0.0214	0.9216	0.6237
3 beta	0	1	0.5000	0.0456	0.0635	0	0.1204
4 gamma	0	1	0.5000	0.6895	1.0696	0	0.2518

5.2 Identifiability analysis

T_newRange=gsua_ia(T,T.Estlsqc)

T_newRange = 4×7 table 

	Range				Nominal	Si	STi
	undefined	undefined	undefined	undefined	Nominal	Si	STi
1 S0	336.3700	336.4517	336.3700	0.1699	0.1339	0	0.0232
2 I0	171.4599	171.4726	171.4599	0.0220	0.0214	0.9216	0.6237
3 beta	0.1290	0.1291	0.1291	0.0456	0.0635	0	0.1204
4 gamma	0.0326	0.0326	0.0326	0.6895	1.0696	0	0.2518

gsua_ia does considerably more than narrow the ranges. The paired feature 'cost' -- give it the res output of gsua_pe -- screens out runs that converged badly before any statistic is computed (see gsua_costcutoff, section 7.2). 'cluster' runs spectral clustering over the repeated estimations to detect multiple global minima, that is, genuinely distinct basins the optimizer landed in, rather than assuming every run scatters around one point; when a real split is found, the reported ranges, correlations and indices describe the dominant basin only, since pooling separated basins into one interval is not meaningful. gsua_dia reports the same statistics without drawing the plots.

[Tia,clusterInfo] = gsua_ia(T,T.Estlsqc,false,false,true,true,'cost',res);

clusterInfo.NumClusters % greater than 1 means distinct basins were found

clusterInfo.Centers % one candidate global minimum per basin

Tia.Est % the accepted ensemble the statistics came from

6. Uncertainty Interval Calculus

6.1 gsua_oatr function

help gsua_oatr

  Function for once-at-time ranges expantion
 
  T=gsua_oatr(T)
  Parameters:
  T <-- summary table from gsua_dataprep function. The objective output is
  given by T.Nominal.
  Outputs:
  T <-- summary table with new factor intervals in T.Range
  Additional features:
  To define limit for range expantion use the paired feature 'lim'.
  Default:0.3.
  To speed up range calculus, use the paired feature 'parallel',true.
  T=gsua_oatr(T,'lim',0.2,'parallel',true)

T.Nominal=[1000 200 0.1 0.06]';

T_oat=gsua_oatr(T)

Expansion-reduction Method OAT is being launched
------------------------------------
Calculus of gamma Starting
Activating intern_counter for gamma
Activating intern_counter for gamma
Range for gamma Done!
Initial range of gamma(4)--> 0  1
Actual range --> 0       0.135
Calculus of beta Starting
Activating intern_counter for beta
Activating intern_counter for beta
Range for beta Done!
Initial range of beta(3)--> 0  1
Actual range --> 0       0.225
Calculus of I0 Starting
Range for I0 Done!
Initial range of I0(2)--> 1  500
Actual range --> 137.375      265.5931
Calculus of S0 Starting
Activating intern_counter for S0
Activating intern_counter for S0
Range for S0 Done!
Initial range of S0(1)--> 100  3000
Actual range --> 0  2250

T_oat = 4×7 table 

	Range				Nominal	Si	STi
	undefined	undefined	undefined	undefined	Nominal	Si	STi
1 S0	0	2250	1000	0.1699	0.1339	0	0.0232
2 I0	137.3750	265.5931	200	0.0220	0.0214	0.9216	0.6237
3 beta	0	0.2250	0.1000	0.0456	0.0635	0	0.1204
4 gamma	0	0.1350	0.0600	0.6895	1.0696	0	0.2518

6.2 gsua_csb function

Please, do not execute gsua_csb function inside live script environment.

help gsua_csb
  Function for uncertainty-based confidence intervals
 
  T= gsua_csb(T,N)
  Parameters:
  T <-- summary table from gsua_oatr function
  N <-- Number of samples per cycle
  Outputs:
  T <-- summary table with new confidence intervals in T.Range
  Additional features:
  It is possible to request the following additional positional outputs
  New_range <-- array with range modification record.
  J_test <-- MSE output of last iteration
  Y_test <-- output of last iteration
  sup    <-- record of good scalars
  Also, it is possible to apply the next paired features
  'ynom',ynom <-- to especify a new output objective
  'reps',k <--to perform k cycles. Default: 100.
  'recort',r <-- to set another recort criteria (r). Default: 0.5.
  'select',s <-- to set another select criteria (s). Default: 0.5.
  'lim',l <-- to set another distance criteria (l). Default: 0.3.
  'stop',p <-- to set the confidence of the interval in (p). Default: 0.95.
  'parallel',false <-- to avoid the speed up of parallel computing toolbox.
  [T,New_range,J_test,Y_test,sup]= gsua_csb(T,N,'lim',0.5,'select',0.4)

Execute the following in matlab command window

T_uci=gsua_csb(T_oat,100)

7. Fit quality, acceptance and confidence intervals

Sections 5 and 6 give you an estimate and a range around it. The functions below answer the questions that follow: was the fit good enough to believe, which multistart runs deserve to be kept, how wide should the acceptance band actually be, and what is the confidence interval of each factor once the others are allowed to move.

7.1 gsua_covmetric function

help gsua_covmetric

gsua_covmetric scores a Monte-Carlo band against the data with two numbers on the same normalized scale: cost_data, the distance from the band median to the data (accuracy), and cost_band, the distance between the band edges (precision, independent of the data). A cost_band far larger than cost_data is the signature of factor confounding: the band is wide because the sampled factors are inconsistent with each other, not because the data is uninformative. That is exactly the condition the joint sampling in section 2.1 exists to fix.

M = gsua_dmatrix(T,2000);

Y = gsua_ua(M,T,'xdata',xdata);

[cost_data,cost_band] = gsua_covmetric(Y,yexp);

7.2 gsua_costcutoff function

help gsua_costcutoff

A multistart run that converged to a much worse cost than the best one contributes essentially arbitrary factor values. Left in the pool it can make a genuinely well-identified factor look poorly identified, purely because one optimizer run failed. gsua_costcutoff decides which runs to keep, either by relative tolerance (rtol) or by the largest gap in the sorted costs (gap), with a minKeep floor so that it can never empty the pool.

keep = gsua_costcutoff(res,'rtol',0.1);

Est_kept = T.Estlsqc(:,keep);

7.3 gsua_noisefloor function

help gsua_noisefloor

A widespread idiom for deciding which fits to accept is res < 1.5*res(1). That threshold is scale-dependent: it asks how close a run is to the best run, not whether it is statistically distinguishable from the truth. So it tightens as the fit improves, and can end up rejecting curves that describe the data perfectly well.

gsua_noisefloor replaces it with a parametric bootstrap. It simulates the best fit, estimates an observation-noise model from that fit's residuals (Poisson, quasi-Poisson by default, or a global negative-binomial), generates synthetic datasets under that noise, and scores the true model against each one using the same cost function gsua_pe minimized. The upper quantile of the resulting distribution is a data-calibrated acceptance threshold rather than a fit-quality-dependent one.

gsua_pe records the margin and alpha it scored with on the table itself, so gsua_noisefloor recovers them automatically. You do not have to restate them, and therefore cannot accidentally score the floor with a different setting than the fit.

[T,res] = gsua_pe(T,xdata,yexp,'solver','lsqc','N',50,'margin',0.1);

out = gsua_noisefloor(T,xdata,yexp,T.Estlsqc,res);

out.threshold % the calibrated cutoff

out.accepted % logical index, res < out.threshold

out.byModel % all three noise models, for a sensitivity check

% 'cumulative' handles cumulative outputs by differencing to incident series

% 'seed' makes the bootstrap reproducible

If the best fit in the pool itself exceeds the noise floor, the function warns: at that level the model is rejected as a description of the data.

7.4 gsua_likelihood function

help gsua_likelihood
Setting environment to work with simulink
Setting nominal values
All done!
T = 4×2 table 
RangeNominal
undefined
1 beta0.40000.6000
2 delta0.40000.6000
3 ro2.20002.8000
4 P02.70003.2000
time = 391×1
    1.0000
    1.1000
    1.2000
    1.3000
    1.4000
    1.5000
    1.6000
    1.7000
    1.8000
    1.9000

	Range	Nominal
1 beta	0.4000	0.6000
2 delta	0.4000	0.6000
3 ro	2.2000	2.8000
4 P0	2.7000	3.2000

Profile-likelihood confidence intervals. Each factor is scanned across a range while all the others are re-estimated at every step, and the interval is the region where the profiled cost stays below a threshold. This is slower than the distribution-free interval gsua_ia reports, but because it re-fits the remaining factors it accounts for the trade-offs between them, which makes it the better choice precisely when factors are correlated. Pass a vector of factor indices to profile only a subset instead of all of them.

8. Where the other functions fit

The functions above are the main workflow. These fill it in:

gsua_dia - the same identifiability statistics as gsua_ia, without the diagnostic plots.
gsua_oatr2 - second-generation one-at-a-time range refinement, the successor to gsua_oatr in section 6.1.
gsua_medianCI - the distribution-free confidence interval of the median that gsua_ia uses to build T.Range.
gsua_depth - band depth, for ranking curves within an ensemble.
gsua_MCF - Monte-Carlo filtering on its own, outside of gsua_ua.
gsua_userdefined - build the summary table T straight from a user function, bypassing gsua_dataprep.
gsua_deval and gsua_pardeval - single and parallel model evaluation, used internally by the analysis functions.
gsua_costf, gsua_rcostf, gsua_costfMulti, gsua_likecost - the cost functions gsua_pe minimizes.

9. Python port

A Python port of this toolbox is maintained alongside it, in the python folder of the repository (package gsua_csb). It mirrors the functions above under idiomatic names with gsua_-prefixed aliases kept for cross-reference: design_matrix and gsua_dmatrix, parameter_estimation and gsua_pe, identifiability_analysis and gsua_ia, noise_floor and gsua_noisefloor, and so on. Simulink models are MATLAB-only; symbolic ODE and user-defined models are both supported there. See python/USERGUIDE.md for its own guide, and python/README.md for what is and is not ported.

Interactive user guide for Toolbox initial configuration (Simulink)

We will use a preconfigured predator - prey model (order: 2, factors: 5) as example throughout the present guide. To access the model, type:

open predatorprey.slx

1. Model preparation

It is advisable to create a mask for factor (input, parameter or initial condition) value assignation. However, it is not necessary.

It is fundamental to assign the value of each factor as a reference to some variable in the Model workspace, it is necessary because the toolbox will identify and assign values for the factors by name matching. If you want to exclude a factor for the analysis, you must to assign it a constant value instead of a reference name. As you can see in the previous figures, the predator-prey model has 5 factors, and we will use all of them but D0 for analysis. If you have preferred a configuration, apply it to the model (use the Model configuration Parameters option, redbox in next figure), otherwise, the toolbox will apply the Simulink current configuration.

At the moment, it is possible to work only with a single model output, hence, you must to connect an output block to the objective output and a terminator block for any other output, as shown in following figure.

For this example, we will apply analysis over prey population (P).

Resources:

Simulink masks - MathWorks

Espacio de estado y simulación - spanish slides

2. Environment configuration and additional features

[T,tout] = gsua_dataprep(model,Ranges,ParIn)

Parameters:

model: file name of the simulink model
Ranges: Npx2 array or cell array with the range of Np factors (see additional features)
ParIn: cell array with the factor names of the model

Outputs:

T: table array with all required system information
tout: array of output time for the model

model='predatorprey';
ParIn={'beta','delta','ro','P0'};
Ranges=[0.4 0.6; 0.4 0.6; 2.2 2.8; 2.7 3.2];
[T,time] = gsua_dataprep(model,Ranges,ParIn)

2.1 Additional features:

2.1.1 Range definition

Paired feature 'rMethod' (default: range)

You can choose between 4 methods to define the ranges of the model factors, each method requires a different definition of variable Ranges

std <-- {Nominal value, Distance to interval limits}
percent <-- {Nominal value, percent of distance to interval limits}
range <-- {Inferior limit, superior limit}
normal <-- {Nominal value, standard deviation}

Examples:

%use of std

Ranges={0.5 0.1; 0.5 0.1; 2.5 0.3; 3 0.2}

[T,time] = gsua_dataprep(model,Ranges,ParIn,'rMethod','std');

%use of percent

s=10;

Ranges={0.5 s; 0.5 s; 2.5 s; 3 s}%you can define a different value for each s

[T,time] = gsua_dataprep(model,Ranges,ParIn,'rMethod','percent');

%use of normal

Ranges={0.5 0.05; 0.5 0.01; 2.5 0.1; 3 0.1}

[T,time] = gsua_dataprep(model,Ranges,ParIn,'rMethod','normal');

Note: The use of normal distribution affects the generation of samples from the space of factors when using the gsua_dmatrix function.

A 'rMethod' for each single range (do not use normal rMethod)

You can define a different rMethod for each factor range. To do this, just add a third column to Range array and specify the rMethod there.

Example:

Ranges={0.5 0.1 'std'; 0.5 0.1 'percent'; 2.5 0.3 'range'; 3 0.2 'percent'};

[T,time] = gsua_dataprep(model,Ranges,ParIn);

2.1.2 Nominal factor family

You can specify a nominal family of factor values using the paired feature 'nominal'. It could be useful for some toolbox functions.

Example:

nominal=[0.6,0.4,2.5,3.2]%array with nominal values for factors

[T,time] = gsua_dataprep(model,Ranges,ParIn,'nominal',nominal);

2.1.3 Model kind (default: Dynamic)

It is possible to work with static models in simulink (initial time = final time = 0). To active the configuration for static models use the paired feature 'modelkind'

Example:

[T,time] = gsua_dataprep(model,Ranges,ParIn,'modelkind','Static');

2.1.4 Fixed step (Default: [])

It is possible to specify from the toolbox a fixed-step magnitude for the model using the paired feature 'Step'

Example:

step=0.1;

[T,time] = gsua_dataprep(model,Ranges,ParIn,'Step',step);

Return to main page

Go to next toolbox features

Interactive user guide for Toolbox initial configuration (Symbolic Matlab)

We will create and configure a simple SIR model with Symbolic Math Toolbox

1. Model creation

syms S(t) I(t) R(t) P(t) beta gamma %symbolic state variables and parameters

P=S+I+R;

%Defining the system of differential equations

ode1 = diff(S) == -S*I*beta/P;

odes(t) =

ode2 = diff(I) == S*I*beta/P - gamma*I;

ode3 = diff(R) == gamma*I;

%Array with the system

odes=[ode1; ode2; ode3]

%choose an order for state variables

Setting environment to work with symbolic Matlab
Introduce ranges in the following order: 

ans(t) =

vars = [R S I];

If your model is in the form

then use the the following command

[odes, vars]=reduceDifferentialOrder(ode,var)% ode=f(), var=y^n;

2. Environment configuration and additional features

[T,solver] = gsua_dataprep(odes,vars,domain,modelName)

Parameters:

odes: array with the system of differential equations
vars: array with state variables
domain: 1x2 array with the range of the domain
modelName: the function creates a file with the model for later simulations, modelName is the name for said output model file

Outputs:

Table: table array with all required system information
solver: function handle for simulations (not necessary)

Use:

Run the following to obtain the correct order for the factors

domain=[0 264];

Setting environment to work with symbolic Matlab
Introduce ranges in the following order: 

ans(t) =

T = 4×2 table 

	Range	Nominal
	undefined	Nominal
1 S0	100	3000
2 I0	1	500
3 beta	0	1
4 gamma	0	1

solver = function_handle with value:
    @(init,pars)ode45(@(t,Y)odefun(t,Y,pars),domain,[Table.Range(fixvars,1)',init],opt)

modelName='SIR';

[T,solver]=gsua_dataprep(odes,vars,domain,modelName);

Now, you can manually add information to the Table T or use the following additional features:

2.1 Additional features

2.1.1 Range definition

Paired feature 'range' and 'rMethod' (Default rMethod: range)

You can directly provide the range of model factors.

You can choose between 4 methods to define the ranges of the model factors, each method requires a different definition of variable Ranges

std <-- {Nominal value, Distance to interval limits}
percent <-- {Nominal value, percent of distance to interval limits}
range <-- {Inferior limit, superior limit}
normal <-- {Nominal value, standard deviation}

Examples:

%use of std

Ranges=[0 0; 1500 1000; 300 200; 0.5 0.5; 0.5 0.5]

[T,sol] = gsua_dataprep(odes,vars,domain,modelName,'range',Ranges,'rMethod','std');

%use of percent

s=10;

Ranges=[0 s; 100 s; 200 s; 0.5 s; 0.5 s]%you can define a different value for each s

[T,sol] = gsua_dataprep(odes,vars,domain,modelName,'range',Ranges,'rMethod','percent');

%use of normal

Ranges=[0 0; 2000 0.1; 300 0.1; 0.5 0.05; 0.5 0.05]

[T,sol] = gsua_dataprep(odes,vars,domain,modelName,'range',Ranges,'rMethod','normal');

Note: The use of normal distribution affects the generation of samples from the space of factors when using the gsua_dmatrix function.

Ranges=[0 0; 100 3000; 1 500; 0 1; 0 1];
[T,solver] = gsua_dataprep(odes,vars,domain,modelName,'range',Ranges)
Setting environment to work with user-defined function
T = 2×2 table 
RangeNominal
undefined
1 10.10000.3000
2 210.000014.0000
sol = function_handle with value:
    @(pars)user_func(pars',domain)

	Range	Nominal
1 1	0.1000	0.3000
2 2	10.0000	14.0000

Note: To fix a factor, give to the function the respective range as a degenerated one (min=max). As you can see, fixed factors will not appear in table T (see previous table for R0 factor).

2.1.2 Nominal factor family

You can specify a nominal family of factors values using the paired feature 'nominal'. It could be useful for some toolbox functions.

Example:

nominal=[0,2000,50,0.1,0.06]%array with nominal values for factors

[T,sol] = gsua_dataprep(odes,vars,domain,modelName,'nominal',nominal);

2.1.3 Specific model output (Default: 1:model_order)

You can specify the output(s) of interest for the model as an array (it should be only one output for SA, and UCI) using the paired feature 'output'.

Example:

%output=[1,3]

output=3;% selected output

[T,sol] = gsua_dataprep(odes,vars,domain,modelName,'output',output);

2.1.4 Specific ode solver (Default: ode45)

You can choose between a fixed step solver (ode4) and an adaptative ode solver (ode45) to solve your model using the paired feature 'solver'.

Example:

solver='ode4';%change the default solver

domain=[0 264 0.1]%add a third column with fixed step size to domain variable

[T,sol] = gsua_dataprep(odes,vars,domain,modelName,'solver',solver);

2.1.5 Ode solver configuration (Default: none)

You can configure the ode45 solver with all possible options (check Matlab odeset) using the paired feature 'opt'.

Example:

opt = odeset('NonNegative',1:3);%nonnegative feature is the only supported for ode4 configuration

[T,sol] = gsua_dataprep(odes,vars,domain,modelName,'opt',opt);

2.1.6 Static model simulation

You can simulate static models defining the domain variable as a single number.

Example:

domain=1;

[T,sol] = gsua_dataprep(odes,vars,domain,modelName);

Return to main page

Go to next toolbox features

Interactive user guide for Toolbox initial configuration (User-Defined functions)

User defined functions must be stored as .m files.

1. Time-dependent functions

1.1 Model creation

Time-dependent functions must be configured in the form output=user_function(factors,domain). Where:

Parameters:

factors: is a vector of parameters of NpxN, being Np the number of parameters and N the number of simulations to run (if vectorized, otherwise, N=1).
domain: is a vector of 1x2 that represents the time domain for the simulation in the form [initial_time, final_time]

Output:

output is a struct() object with the following fields at least:

output.x: simulation time.
output.y: simulation outputs, as a matrix of nOutputsx|output.x|, where nOutputs is the number of outputs of the model and |output.x| is the length of the simulation time.

In this way, the outputs from matlab ode and dde solvers meet the requirements.

1.1.1 Example

function sol = user_dependent(params,domain,in)

% This problem is an epidemic model due to Cooke et alia, more information

% can be found in 'Interaction of maturation delay and nonlinear birth in

% population and epidemic models' J. Math. Biol., 39 (1999) 332-352.

% (This is reference 3 of the tutorial).

% You can Download the tutorial from:

% https://www.mathworks.com/matlabcentral/fileexchange/3899-tutorial-on-solving-ddes-with-dde23

T = params(1,:);

lambda = params(2,:);

sol = dde23(@prob4f,T,[2; 3.5],domain,[],lambda,T);%the output from dde23 meet the requirements to be the output of the function.

%-----------------------------------------------------------------------

function yp = prob4f(t,y,Z,lambda,T)

%PROB4F The derivative function for Problem 4 of the DDE Tutorial.

a = 1;

b = 80;

d = 1;

d1 = 1;

e = 10;

gamma = 0.5;

I = y(1);

N = y(2);

Nlag = Z(2,1);

dIdt = lambda*(N - I)*(I/N) - ( d + e + gamma)*I;

dNdt = b*exp(-a*Nlag)*Nlag*exp(-d1*T) - d*N - e*I;

yp = [ dIdt; dNdt];

1.2 Model implementation

To prepare the environment to work with time-dependent user-defined functions, you must run the following code, replacing input parameters as indicated

[Table,solver]=gsua_dataprep(func,ranges,'domain',domain)

Parameters:

func: name of the .m file as char
ranges: Npx2 array with the range of the factors
domain: 1x2 array with the range of the domain

Outputs:

Table: table array with all required system information
solver: function handle for simulations (not essential)

Use:

We will expose an example implementing the example function user_defined.m (defined above), you can use the command

open user_defined.m

to access the function

func='user_dependent';

ranges=[0.1 0.3; 10 14];

domain=[0 25];

[T,sol]=gsua_dataprep(func,ranges,'domain',domain)

Test if the model is actually working by running a simulation with the values stored at T.Nominal

gsua_eval(T.Nominal,T);

Change the number of outputs (max 2) with the following command

T.Properties.CustomProperties.output=1:2;

Provide a specific xdata to achieve the model behaviour at those points

xdata=0:0.1:25;
Setting environment to work with user-defined function
T = 8×2 table 
RangeNominal
undefined
1 104.0000
2 204.0000
3 301.0000
4 401.0000
5 512.0000
6 600.3000
7 701.0000
8 800.1000
sol = function_handle with value:
    @(pars)user_func(pars')

	Range	Nominal
1 1	0	4.0000
2 2	0	4.0000
3 3	0	1.0000
4 4	0	1.0000
5 5	1	2.0000
6 6	0	0.3000
7 7	0	1.0000
8 8	0	0.1000

Evaluate the model again

gsua_eval(T.Nominal,T,xdata);

1.3 Additional features

1.3.1 Range definition

Paired feature 'rMethod' (Default rMethod: range)

You can choose between 4 methods to define the ranges of the model factors, each method requires a different definition of variable ranges

std <-- {Nominal value, Distance to interval limits}
percent <-- {Nominal value, percent of distance to interval limits}
range <-- {Inferior limit, superior limit}
normal <-- {Nominal value, standard deviation}

Examples:

%use of std

ranges=[1 0.1; 12 2]

[T,sol] = gsua_dataprep(func,ranges,'domain',domain,'rMethod','std');

%use of percent

s=10;

ranges=[0.2 s; 12 s]%you can define a different value for each s

[T,sol] = gsua_dataprep(func,ranges,'domain',domain,'rMethod','percent');

%use of normal

ranges=[0.3 0.001; 12 0.1]

[T,sol] = gsua_dataprep(func,ranges,'domain',domain,'rMethod','normal');

Note: The use of normal distribution affects the generation of samples from the space of factors when using the gsua_dmatrix function.

Note: To fix a factor, give to the function the respective range as a degenerated one (min=max). As you can see, fixed factors will not appear in table T (see previous table for R0 factor).

1.3.2 Nominal factor family

You can specify a nominal family of factors values using the paired feature 'nominal'. It could be useful for some toolbox functions.

Example:

nominal=[0.2,11]%array with nominal values for factors

[T,sol] = gsua_dataprep(func,ranges,'domain',domain,'nominal',nominal);

1.3.3 Specific model output (Default: 1:model_order)

You can specify the output(s) of interest for the model as an array (it should be only one output for SA, and UCI) using the paired feature 'output'.

Example:

%output=[1,2]

output=2;% selected output

[T,sol] = gsua_dataprep(func,ranges,'domain',domain,'output',output);

1.3.4 Specific names for model factors (Default: {})

You can specify the names of the model factors as an array using the paired feature 'names'. The names of the factors are displayed when plotting some graphs.

Example:

names={'\tau','\lambda'};% selected output

[T,sol] = gsua_dataprep(func,ranges,'domain',domain,'names',names);

1.3.5 Specific names for model outputs (Default: {})

You can specify the names of the model output(s) as an array using the paired feature 'out_names'. The names of the outputs are displayed when plotting some graphs.

Example:

out_names={'I(t)','N(t)'};% selected output

[T,sol] = gsua_dataprep(func,ranges,'domain',domain,'out_names',out_names);

1.3.6 Run as vectorized function (Default: false)

If you have a vectorized function, then you can exploit its capabilities using the paired feature 'vectorized'

Example:

[T,sol] = gsua_dataprep(func,ranges,'domain',domain,'vectorized',true);

1.3.7 Any additional user-provided configuration (Default: {})

If your function requires on three inputs, then the third input must be for additional configuration. You can provide anything you want for this function input using the paired feature 'opt'

Example:

opt=struct();%anything you want

[T,sol] = gsua_dataprep(func,ranges,'domain',domain,'opt',opt);

2. Time-free functions

2.1 Model creation

Time-free functions must be configured in the form output=user_function(factors). Where:

Parameters:

factors: a vector that contains all the parameters needed to run the function

Output:

A numerical array (if vectorized) or a single number.

Tip: If your model is time-free but its output is a vector, treat it as a time-dependent model. Just ignore the domain input inside the function.

1.1.1 Example

function out = user_free(params)

%R_0

% ODES = R_0(BETA_M,BETA_H,THETA_M,THETA_H,ALPHA,MU_M,GAMMA_H,MU_H)

% This function was generated by the Symbolic Math Toolbox version 8.3.

% 04-Aug-2019 17:48:42

beta_m=params(1,:);

beta_h=params(2,:);

theta_m=params(3,:);

theta_h=params(3,:);

alpha=params(3,:);

mu_m=params(3,:);

gamma_h=params(3,:);

mu_h=params(3,:);

t2 = alpha.*mu_m;

out = (beta_h.*beta_m.*theta_h.*theta_m)./(t2.*(gamma_h+mu_h).*(mu_h+theta_h).*(t2+theta_m));

2.2 Model implementation

To prepare the environment to work with time-free user-defined functions, you must run the following code, replacing input parameters as indicated

[T,solver] = gsua_dataprep(Function,Ranges)

Parameters:

Function: Is the name of the user-defined function as char ('my_function')
Ranges: Numerical array of Npx2

Outputs:

T: table array with all required system information
solver: function handle for simulations (not necessary)

func='user_free';
range=[0 4; 0 4; 0 1; 0 1; 1 2; 0 0.3; 0 1; 0 0.1];
[T,sol]=gsua_dataprep(func,range)

2.3 Additional features

For additional features see Time-dependent features

Return to main page

Go to next toolbox features