Morris Method / Elementary Effects

Morris Method / Elementary Effects#

Morris’ Method of Elementary Effects.

Module includes optimised sampling of trajectories including optional groups as well as calculation of the Morris measures mu, stddev and mu*.

Saltelli, A., Chan, K., & Scott, E. M. (2000). Sensitivity Analysis. Wiley Series in Probability and Statistics, John Wiley & Sons, New York, 1-504, pages 68ff

Provided functions are:

morris_sampling - Sample trajectories in parameter space

elementary_effects - Calculate Elementary Effects from model output on trajectories

Note that the functions morris_sampling and elementary_effects are wrappers for the functions Optimized_Groups and Morris_Measure_Groups of F. Campolongo and J. Cariboni ported to Python by Stijn Van Hoey.

This module was originally written by Stijn Van Hoey as a translation from an original Matlab code of F. Campolongo and J. Cariboni, JRC - IPSC Ispra, Varese, IT. It was adapted by Matthias Cuntz while at Department of Computational Hydrosystems, Helmholtz Centre for Environmental Research - UFZ, Leipzig, Germany, and continued while at Institut National de Recherche pour l’Agriculture, l’Alimentation et l’Environnement (INRAE), Nancy, France.

copyright:

Copyright 2012-2022 Stijn Van Hoey, Matthias Cuntz, see AUTHORS.rst for details.

license:

MIT License, see LICENSE for details.

The following wrappers are provided

morris_sampling(nparam, lb, ub, nt[, ...])

Sample trajectories in parameter space

elementary_effects(nparam, OptMatrix, ...[, ...])

Compute the Morris measures given the Morris sample matrix, the output values and the group matrix.
History

  • Written in Matlab by F. Campolongo, J. Cariboni, JRC - IPSC Ispra, Varese, IT. Last Update 15 November 2005 by J. Cariboni: http://sensitivity-analysis.jrc.ec.europa.eu/software/index.htm now at: https://ec.europa.eu/jrc/en/samo/simlab

  • Translated to Python in May 2012 by Stijn Van Hoey.

  • Adapted to Python 3, etc., Oct 2013, Matthias Cuntz

  • Went from exponential time increase with number of trajectories to linear increase by using in Optimised_Groups one call to cdist from scipy.spatial.distance and removed one loop in a loop over total number of trajectories. Several further little improvements on speed, Dec 2017, Matthias Cuntz

  • Allow single trajectories, Dec 2017, Matthias Cuntz

  • Catch degenerated case where lower bound==upper bound, return 0, Feb 2018, Matthias Cuntz

  • Use integer division operator // instead of / for trajectory length r, Jul 2018, Fabio Genaretti

  • Distance matrix is not done for all trajectories at once because of very large memory requirement, Aug 2018, Matthias Cuntz & Fabio Gennaretti

  • Changed to Sphinx docstring and numpydoc, Dec 2019, Matthias Cuntz

  • Distinguish iterable and array_like parameter types, Jan 2020, Matthias Cuntz

  • Remove np.matrix in Sampling_Function_2, called in Optimized_Groups to remove numpy deprecation warnings, Jan 2020, Matthias Cuntz

  • Plot diagnostic figures in png files if matplotlib installed, Feb 2020, Matthias Cuntz

  • Renamed file to morris_method.py, Feb 2020, Matthias Cuntz

  • Adjusted argument and keyword argument names to be consistent with pyeee, Feb 2020, Matthias Cuntz

  • Make number of final trajectories an argument instead of a keyword argument and sample default of 10*final trajectories, Feb 2020, Matthias Cuntz

  • Sample not only from uniform distribution but allow all distributions of scipy.stats, Mar 2020, Matthias Cuntz

  • More consistent docstrings, Jan 2022, Matthias Cuntz

  • Raise Error if more than one component changed at once, Jul 2023, Matthias Cuntz

elementary_effects(nparam, OptMatrix, OptOutVec, Output, nsteps=4, Group=[], Diagnostic=False)[source]#

Compute the Morris measures given the Morris sample matrix, the output values and the group matrix.

Parameters:

  • nparam (int) – Number of parameters / factors list [OptMatrix, OptOutVec] with OptMatrix((nparam+1)*nt,nparam) and OptOutVec(nparam*nt)

  • OptMatrix (ndarray) – ((nparam+1)*nt,nparam) Matrix of the Morris sampled trajectories from morris_sampling

  • OptOutVec (ndarray) – (nparam*nt,) Matrix with the parameter / factor changings from morris_sampling

  • Output (ndarray) – ((nparam+1)*nt,) Matrix of the output values of each point of each trajectory

  • nsteps (int, optional) – Number of levels, i.e. intervals in trajectories (default: 4)

  • Group (ndarray, optional) – (nparam,NumGroups) Matrix describing the groups. (default: []) Each column represents a group. The elements of each column are zero if the parameter / factor is not in the group, otherwise it is 1.

  • Diagnostic (bool, optional) – Print out diagnostics if True (default: False)

Returns:

SA, OutMatrix – SA(nparam*Output.shape[1],N) individual sensitivity measures, OutMatrix(nparam*Output.shape[1], 3) = [Mu*, Mu, StDev] Morris Measures

It gives the three measures of each parameter / factor for each output.

Return type:

list of ndarrays

References

Saltelli, A., Chan, K., & Scott, E. M. (2000). Sensitivity Analysis. Wiley Series in Probability and Statistics, John Wiley & Sons, New York, 1-504. - on page 68ff

Examples

>>> import numpy as np
>>> seed = 1023
>>> np.random.seed(seed=seed)
>>> npara = 10
>>> x0    = np.ones(npara)*0.5
>>> lb    = np.zeros(npara)
>>> ub    = np.ones(npara)
>>> mask  = np.ones(npara, dtype=bool)
>>> mask[5::2] = False
>>> nmask = np.sum(mask)
>>> nt     = npara
>>> ntotal = max(nt**2, 10*nt)
>>> nsteps = 6
>>> tmatrix, tvec = morris_sampling(nmask, lb[mask], ub[mask], nt,
...                                 ntotal=ntotal, nsteps=nsteps,
...                                 Diagnostic=False)
>>> # Set input vector to trajectories and masked elements to x0
>>> x = np.tile(x0, tvec.size).reshape(tvec.size, npara)  # default to x0
>>> x[:,mask] = tmatrix  # replaced unmasked with trajectorie values
>>> func = np.sum
>>> fx = np.array(list(map(func,x)))
>>> out = np.zeros((npara,3))
>>> sa, res = elementary_effects(nmask, tmatrix, tvec, fx, nsteps=nsteps,
...                              Diagnostic=False)
>>> out[mask,:] = res
>>> print(out[:,0])
[1. 1. 1. 1. 1. 0. 1. 0. 1. 0.]
morris_sampling(nparam, lb, ub, nt, nsteps=6, ntotal=None, dist=None, distparam=None, GroupMat=array([], dtype=float64), Diagnostic=0)[source]#

Sample trajectories in parameter space

Optimisation in the choice of trajectories for Morris experiment, that means elementary effects.

Parameters:

  • nparam (int) – Number of parameters / factors

  • lb (array_like) – (nparam,) Lower bound of the uniform distribution for each parameter / factor or lower fraction of percent point function ppf if distribution given. Be aware that the percent point function ppf of most continuous distributions is infinite at 0.

  • ub (array_like) – (nparam,) Upper bound of the uniform distribution for each parameter / factor or upper fraction of percent point function ppf if distribution given. Be aware that the percent point function ppf of most continuous distributions is infinite at 1.

  • nt (int) – Final number of optimal trajectories

  • nsteps (int, optional) – Number of levels, i.e. intervals in trajectories (default: 6)

  • ntotal (int, optional) – Total number of sampled trajectories. If None: ntotal=10*nt (default: None)

  • dist (list, optional) – List of None or scipy.stats distribution objects for each factor having the method ppf, Percent Point Function (Inverse of CDF) (default: None) If None, the uniform distribution will be sampled from lower bound lb to upper bound ub. If dist is scipy.stats.uniform, the ppf will be sampled from the lower fraction given in lb and the upper fraction in ub. The sampling interval is then given by the parameters loc=lower and scale=interval=upper-lower in distparam. This means dist=None, lb=a, ub=b corresponds to lb=0, ub=1, dist=scipy.stats.uniform, distparam=[a,b-a]

  • distparam (list, optional) – List with tuples with parameters as required for dist (default: (0,1)). All distributions of scipy.stats have location and scale parameters, at least. loc and scale are implemented as keyword arguments in scipy.stats. Other parameters such as the shape parameter of the gamma distribution must hence be given first, e.g. (shape,loc,scale) for the gamma distribution. distparam is ignored if dist is None. The percent point function ppf is called like this: dist(*distparam).ppf(x)

  • GroupMat (ndarray, optional) – (nparam,ngroup) Matrix describing the groups. (default: np.array([])) Each column represents a group. The elements of each column are zero if the parameter / factor is not in the group, otherwise it is 1.

  • Diagnostic (bool, optional) – Plot the histograms and compute the efficiency of the sampling if True (default: False)

Returns:

traj – list [OptMatrix, OptOutVec] with OptMatrix((nparam+1)*nt,nparam) and OptOutVec(nparam*nt)

Return type:

list

References

Saltelli, A., Chan, K., & Scott, E. M. (2000). Sensitivity Analysis. Wiley Series in Probability and Statistics, John Wiley & Sons, New York, 1-504. - on page 68ff

Examples

>>> import numpy as np
>>> seed = 1023
>>> np.random.seed(seed=seed)
>>> npara = 10
>>> x0    = np.ones(npara)*0.5
>>> lb    = np.zeros(npara)
>>> ub    = np.ones(npara)
>>> mask  = np.ones(npara, dtype=bool)
>>> mask[5::2] = False
>>> nmask = np.sum(mask)
>>> nt     = npara
>>> ntotal = max(nt**2, 10*nt)
>>> nsteps = 6
>>> tmatrix, tvec = morris_sampling(nmask, lb[mask], ub[mask], nt,
...                                 nsteps=nsteps, ntotal=ntotal,
...                                 Diagnostic=False)
>>> # Set input vector to trajectories and masked elements to x0
>>> x = np.tile(x0, tvec.size).reshape(tvec.size, npara)  # default to x0
>>> x[:,mask] = tmatrix  # replaced unmasked with trajectorie values
>>> print(x[0,:])
[0.6 0.4 0.8 0.6 0.6 0.5 0.4 0.5 0.  0.5]

>>> import scipy.stats as stats
>>> seed = 1023
>>> np.random.seed(seed=seed)
>>> npara = 10
>>> x0    = np.ones(npara)*0.5
>>> lb    = np.zeros(npara)
>>> ub    = np.ones(npara)
>>> dist  = [ stats.uniform for i in range(npara) ]
>>> distparam = [ (lb[i], ub[i]-lb[i]) for i in range(npara) ]
>>> mask  = np.ones(npara, dtype=bool)
>>> mask[5::2] = False
>>> nmask = np.sum(mask)
>>> nt     = npara
>>> ntotal = max(nt**2, 10*nt)
>>> nsteps = 6
>>> tmatrix, tvec = morris_sampling(nmask, lb[mask], ub[mask], nt,
...                                 nsteps=nsteps, ntotal=ntotal,
...                                 dist=dist, distparam=distparam,
...                                 Diagnostic=False)
>>> # Set input vector to trajectories and masked elements to x0
>>> x = np.tile(x0, tvec.size).reshape(tvec.size, npara)  # default to x0
>>> x[:,mask] = tmatrix  # replaced unmasked with trajectory values
>>> print(x[0,:])
[0.6 0.4 0.8 0.6 0.6 0.5 0.4 0.5 0.  0.5]
Contents