Monte Carlo Bayesian Sensitivity Analysis
Description
Fully Bayesian Monte Carlo sensitivity analysis scheme, based upon any of the regression models contained in the tgp package. Random Latin hypercube samples are drawn at each MCMC iteration in order to estimate main effects as well as 1st order and total sensitivity indices.
Usage
1 2 3 
Arguments
X 

Z 
Vector of output responses 
nn.lhs 
Size of each Latin hypercube drawn for
use in the Monte Carlo integration scheme. Total number of locations
for prediction is 
model 
Either the regression model used for prediction, or

ngrid 
The number of grid points in each input dimension upon which main effects will be estimated. 
span 
Smoothing parameter for main effects integration:
the fraction of 
BTE 
3vector of MonteCarlo parameters (B)urn in, (T)otal, and (E)very. Predictive samples are saved every E MCMC rounds starting at round B, stopping at T 
rect 
Rectangle describing the domain of the uncertainty
distribution with respect to which the sensitivity is to be
determined. This defines the domain from which the LH sample
is to be taken. The rectangle should be a 
shape 
Optional vector of shape parameters for the Beta
distribution. Vector of length equal to the dimension of the domain,
with elements > 1. If specified, the uncertainty distribution
(i.e. the LH sample) is proportional to a joint pdf formed by
independent Beta distributions in each dimension of the domain,
scaled and shifted to have support defined by 
mode 
Optional vector of mode values for the Beta
uncertainty distribution. Vector of length equal to the dimension
of the domain, with elements within the support defined by

... 
Extra arguments to the tgp 
Details
Saltelli (2002) describes a Latin Hypercube sampling based method for estimation of the 'Sobal' sensitivity indices:
1st Order for input i,
S(i) = var(E[fx[i]])/var(f),
where x[i] is the ith input.
Total Effect for input i,
T(i) = E[var(fx[i])]/var(f),
where x[i] is all inputs except for the ith.
All moments are with respect to the appropriate marginals of the
uncertainty distribution U – that is, the probability
distribution on the inputs with respect to which sensitivity is being
investigated.
Under this approach, the integrals involved are approximated through
averages over properly chosen samples based on two LH samples
proportional to U. If nn.lhs
is the sample size for the
Monte Carlo estimate, this scheme requires nn.lhs*(ncol(X)+2)
function evaluations.
The sens
function implements the method for unknown functions
f, through prediction via one of the tgp regression
models conditional on an observed set of X
locations.
At each MCMC iteration of the tgp model fitting,
the nn.lhs*(ncol(X)+2)
locations are drawn randomly from the
LHS scheme and realizations of the sensitivity indices are
calculated. Thus we obtain a posterior sample of the indices,
incorporating variability from both the Monte Carlo estimation and
uncertainty about the function output. Since a subset of the
predictive locations are actually an LHS proportional to the
uncertainty distribution, we can also estimate the main effects
through simple nonparametric regression (a moving average).
Please see vignette("tgp2")
for a detailed illustration
Value
The output is a "tgp"
class object. The details for btgp
contain a complete description of this output. The list entry that
is relevance to sensitivity analysis is sens
, which itself has entries:
par 
This contains a 
Xgrid 
A 
ZZ.mean 
A 
ZZ.q1 
A 
ZZ.q2 
A 
S 
A 
T 
A 
Note
The quality of sensitivity analysis is dependent on the size of
the LH samples used for integral approximation; as with any Monte
Carlo integration scheme, the sample size (nn.lhs
) must
increase with the dimensionality of the problem. The total
sensitivity indices T are forced nonnegative,
and if negative values occur it is necessary to increase
nn.lhs
. The plot.tgp
function replaces negative values with zero
for illustration.
Author(s)
Robert B. Gramacy, rbgramacy@chicagobooth.edu, and Matt Taddy, taddy@chicagobooth.edu
References
R.D. Morris, A. Kottas, M. Taddy, R. Furfaro, and B. Ganapol. (2009) A statistical framework for the sensitivity analysis of radiative transfer models. IEEE Transactions on Geoscience and Remote Sensing, to appear.
Saltelli, A. (2002) Making best use of model evaluations to compute sensitivity indices. Computer Physics Communications, 145, 280297.
See Also
btgp
, plot.tgp
,
predict.tgp
, lhs
Examples
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  # Take a look at the air quality in New York: Sensitivity of
# ozone levels with respect to solar radiation, wind, and
# temperature. See help(airquality) for details.
X < airquality[,2:4]
Z < airquality$Ozone
# Uncertainty distribution is the default: uniform over range(X)
# There is missing data, which is removed automatically by tgp
# range(X).
s < suppressWarnings(sens(X=X, Z=Z, nn.lhs=300, model=btgp,
ngrid=100, span=0.3, BTE=c(5000,10000,10)))
# plot the results
plot(s, layout="sens", ylab="Ozone", main="main effects")
# plot only the sensitivity indices
plot(s, layout="sens", ylab="Ozone", maineff=FALSE)
# plot only the main effects, side by side
plot(s, layout="sens", ylab="Ozone", main="", maineff=t(1:3))
# build a 'sens.p' parameter vector for a datadependent
# informative uncertainty distribution. For each variable,
# the input distribution will be a scaled Beta with shape=2,
# and mode equal to the data mean
rect < t(apply(X, 2, range, na.rm=TRUE))
mode < apply(X , 2, mean, na.rm=TRUE)
shape < rep(2,3)
# plot a sample from the marginal uncertainty distribution.
Udraw < lhs(300, rect=rect, mode=mode, shape=shape)
par(mfrow=c(1,3))
for(i in 1:3) hist(Udraw[,i], breaks=15,xlab=names(X)[i])
# build sens.p with the 'sens' function.
sens.p < suppressWarnings(sens(X=X,Z=Z,nn.lhs=300, model=NULL,
ngrid=100, rect=rect, shape=shape, mode=mode))
# Use predict.tgp to quickly analyze with respect to this new
# uncertainty distribution without rerunning the MCMC, then
# plot the results.
s.new < predict(s, BTE=c(1,1000,1), sens.p=sens.p, verb=1)
plot(s.new, layout="sens", ylab="Ozone", main="main effects")
