Tess Final Meeting Capri

1
From Monte Carlo filtering and Regional
Sensitivity Analysis (RSA) to Generalised
Likelihood Uncertainty Estimation (GLUE)
and Tree-Structured Density Estimation
(TSDE) M. Ratto, IPSC-JRC, European Commission
2
Background

• environmental sciences, early 80s
• complex numerical and analytical models, based
  on first principles, conservation laws,
• ill-defined parameters, competing model
  structures
• (different constitutive equations, different
  types of process considered, spatial/temporal
  resolution, )
• need to establish magnitude and sources of
  prediction uncertainty
• Monte Carlo simulation analyses.

to achieve a better understanding of the
simulated system to increase the reliability of
model predictions to define realistic values to
be used in subsequent risk assessment
3
MC filtering and RSA
One of the earliest landmark application of MC
simulation eutrophication in the Peel Inlet, SW
Australia. Regional Sensitivity Analysis
developed and employed. Hornberger G.M., Spear
R.C., An approach to the preliminary analysis of
environmental systems, J. Environm. Management,
12, 7-18, 1981. Hornberger G.M., Spear R.C.,
Eutrophication in Peel Inlet, I. The
problem-defining behavior and a mathematical
model for the phosphorous scenario, Water
Research, 14, 29-42, 1980. Spear R.C. and
Hornberger G.M., Eutrophication in Peel Inlet,
II. Ideintification of Critical Uncertainties via
Generalised Sensitivity Analysis, Water Research,
14, 43-49, 1980.
4
MC filtering and RSA

• Two tasks for RSA
• qualitative definition of the system behaviour
• a set of constraints thresholds, ceilings, time
  bounds based on available information on the
  system
• binary classification of model outputs based on
  the specified behaviour definition.
• qualifies a simulation as behaviour (B) if the
  model output lies within constraints,
  non-behaviour ( ) otherwise

5
MC filtering and RSA
Define a range for k input factors xi 1 ? i ?
k, reflecting uncertainty in parameters and
model constituent hypotheses. Each Monte Carlo
simulation is associated to a vector of values of
the input factors. Classifying simulations as
either B or , a set of binary elements is
defined allowing to distinguish two sub-sets for
each xi
6
MC filtering and RSA
The Kolmogorov-Smirnov two-sample test (two-sided
version) is performed for each factor
independently
Test statistic
where F are marginal cumulative probability
functions, f are probability density functions
7
MC filtering and RSA
At what significance level does the computed
value of dm,n determine the rejection of Ho?
8
MC filtering and RSA
The importance of the uncertainty of each
parameter is (inversely) related to this
significance level. Input factors are grouped
into three sensitivity classes, based on the
significance level for rejecting Ho 1
critical 2 important 3 insignificant.
9
RSA comment
RSA has many global properties, similarly to
variance based methods The whole range of value
of input factors is considered, all factors are
varied at the same time RSA classification is
related to Main Effects of Variance Based
methods (it analyses univariate marginal
distributions). No higher order analysis is
performed with RSA approach, searching for
interaction structure.
10
RSA problems

• Spear et al. (1994) reviewed their experience
  with RSA, highlighting two key drawbacks
• (1) The success rate
• the fraction of B is ?5 over the total
  simulations for large models (kgt20)
• lack in statistical power

11
RSA problems (ctd.)

• (2) Correlation and interaction structures of the
  B subset ? also Becks review, 1987
• Smirnov test is a sufficient test only if Ho is
  rejected (i.e. the factor is IMPORTANT)
• any covariance structure induced by the
  classification is not detected by the univariate
  dm,n statistic.
• e.g. factors combined as products or quotients
  may compensate
• bivariate correlation analysis is not revealing,
  either.

12
RSA problems (ctd.)
Example 1. Y X1X2, X1,X2 ?U-0.5,0.5 Behaviou
ral runs Ygt0.
Scatter plot of the B subset in the X1,X2 plane
Smirnov test fails Correlation analysis would
help (e.g. PCA)
13
RSA problems (ctd.)
Example 1. Cumulative distributions of
X1 (Behavioural runs Ygt0.)
14
RSA problems (ctd.)
Example 2. Y X12 X22, X1,X2
?U-0.5,0.5 Behavioural runs 0.2 lt Y lt 0.25
Scatter plot of the B subset in the X1,X2 plane
Smirnov test fails Correlation fails as
well (sample ??-0.04)
15
Tree-Structured Density Estimation (TSDE)

• TSDE Multivariate analysis of
• posterior B parameter distributions.
• Sequence of recursive binary splits of the B
  region into sub-domains (similarly to peaks and
  tails of a histograms)
• small regions of relatively high-density
• larger sparsely populated regions

Spear R.C., Grieb T.M. and Shang N., Factor
Uncertainty and Interaction in Complex
Environmental Models, Water Resources Research,
30, 3159-3169, 1994.
16
TSDE algorithm
The m vectors of factors (xiB) are samples from
a population of unknown density function. The
aim is to make a local approximation in the
k-dimensional domain ?. The density estimate of
a sub-domain ?i and the loss function L of the
density estimation are
17
TSDE algorithm
Each parameter range in ? is divided into p bins,
so defining (k x p) ways to make a split
For each trial split L(0) - L(1) measures the
increase in accuracy of the estimated
density. The TSDE searches for the optimal split
to maximise (L(0) - L(1) ).
18
TSDE algorithm

• The algorithm stops if
• the maximum increase in accuracy is less than a
  given tolerance
• the sample size in each sub-domain is smaller
  than some critical number.

19
TSDE algorithm
Example from Spear et al., Wat. Res. Research,
1994.
20
TSDE algorithm
Example from Spear et al., 1994.
21
TSDE algorithm
Example from Spear et al., 1994.
22
TSDE algorithm
From the binary tree structure we can make useful
inferences

• of high-density terminal nodes indicate the
  of multiple optima their combined volume of ?,
  indicates the probability of realising B.

• The of high-density terminals that each factor
  defines indicates its importance in matching B
  the higher in the tree, the more it is important

• tracing down to each high-D terminal node will
  reveal the key factors, illustrating the
  interaction structure among the B-giving factors.

23
Generalised Likelihood Uncertainty Estimation
(GLUE)
GLUE is an extension of RSA
Beven K.J., Binley A., The Future of Distributed
Models Model Calibration and Uncertainty
Prediction, Hydrological Processes, 6, 279-298,
1992. Romanowicz R., Beven K., Tawn J.,
Evaluation of predictive uncertainty in nonlinear
hydrological models using a Bayesian approach. In
Statistics for the Environment 2, Water Related
Issues. Barnett V., Turkman F. (eds) pp.
297-315, Wiley New York, 1994.
24
GLUE
- developed from the acceptance of the possible
equifinality of models
- works with multiple sets of factors (via MC
sampling)
- model runs weighted and ranked on a likelihood
scale via conditioning on observations and the
weights are used to formulate a cumulative
distribution of predictions.
Comparing with RSA, no filtering is made
(on-off), BUT a smooth function is used to
evaluate behaviour of simulations.
25
What is GLUE
data
MC model runs
Y(t)
Likelihood measure
Uncertainty estimation
26
What is GLUE
Most typical likelihood measures where
Y(t) is the set of observations ?i is the i-th
set of parameters/structure
is the mean square
error
for the i-th run ...
27
A role for GSA?
As in RSA, complex interaction structures are
typical in GLUE complex structured factor
subsets having similar likelihood of being
simulators of reality are usually identified.
28
A role for GSA?

• When trying to map back the output uncertainty to
  the input factors, as for RSA

• conventional multivariate statistics (PCA ) is
  not very revealing (Spear et al., 1994)

• complex parametric interactions do not become
  evident from looking at univariate marginal
  probability densities (Spear, 1997)

? perform GSA on the likelihood measure GSA-GLUE
Ratto et al., Computer Physics Communications
(2001)
29
Methods GSA-GLUE
Likelihood weights methodological problems for SA

• non monotonic input-output mapping

• high level of interaction between input factors.

? variance-based Global Sensitivity Analysis
(GSA) quantitative methods, without restrictions
about monotony or additivity of the model.
30
Why GSA?

• The determination of both main effect (as other
  methods) and total effect allows a
  classification, e.g.

• factors with high main effect
• affect model output singularly, independently of
  interaction

• factors with small main effect but high total
  effect
• influence on the model output mainly through
  interaction

• factors with small main and total effect
• have a negligible effect on the model output and
  can be fixed or neglected (? more powerful than
  RSA!).

31
How do we do GSA-GLUE?
The sample generated for the GLUE analysis is
designed also for GSA. (Sobol or FAST sample)
In this way, with the same set of model runs
- predictive uncertainty can be estimated
- sensitivity indices computed and
- analysis of the interaction structure
performed.
32
Comments on the likelihood in GLUE
The definition of likelihood function is
somehow arbitrary (to be rigorous, we should
better refer to it as a loss function).
This arbitrariness is sometimes seen as a
limitation for the robustness of the overall
analysis. Specifically, to extend GLUE in the
econometric context, we deem necessary to
introduce a rigorous definition and unambiguous
definition.
33
Comments on the likelihood in GLUE
In our macro-economic applications, we apply
the true likelihood function, as defined by
the FIML estimator. In particular, we found that
a very good compromise is to apply the
concentrated likelihood, which allows a rigorous
and unambiguous definition of the weights for
GLUE analysis.
34
TEST CASE
Chemical system an isothermal first order
irreversible reaction in a batch system A?B
35
TEST CASE
Model solution
Model parameters
36
TEST CASE
Input factors Outputs 1) yB(t) (time
dependent output) 2) a likelihood weight, as
defined in GLUE
37
TEST CASE
Likelihood weights ( better to say loss
function)
where
is the mean square error for the i-th run
38
TEST CASE
Time dependent output Experimental
time series and 5 and 95 confidence bounds for
the output yB.
39
TEST CASE
Sensitivity of yB(t)
main effects total effects
? main effects sum up to 1 (small interaction) ?
kinetics mainly dominates the process
40
Memo

• Total sensitivity measure

example for three factors YY(A,B ,C)
41
TEST CASE
Sensitivity of likelihood 1/s2
?no predominance in main effect ?kinetic
parameters have big interaction
42
TEST CASE
Scatter plots ? Main effect sensitivity indices
? RSA
43
Conclusions

• initial condition is unimportant (smallest ST)
• chemical rate factors mainly drive the model fit
  to the experimental data (highest Si and SiT)
• BUT, chemical rate factors cannot be precisely
  estimated, because their effect is dominated by
  interaction terms
• high correlation between estimates of k? and E
  is usually detected

44
What does GSA tell us?
? an analysis limited to main effects would tell
us nothing more than dotty plots/marginals, etc.
? an high difference between main and total
effects measures the level of interaction in the
model
45
Interaction

• Possible approaches
• Correlation analysis of the posterior pdf (B
  subset)
• Variance based GSA
• Spears Tree-Structured Densitiy Estimation
  technique
• Bootstrapping
• Graphical modelling / Bayesian Networks (GKSS)