{"id":416,"date":"2014-01-02T16:40:21","date_gmt":"2014-01-02T21:40:21","guid":{"rendered":"http:\/\/homepages.uc.edu\/~yaozo\/wordpress\/?p=416"},"modified":"2014-01-02T16:40:21","modified_gmt":"2014-01-02T21:40:21","slug":"monte-carlo-bayesian-sensitivity-analysis","status":"publish","type":"post","link":"https:\/\/zhuoyao.net\/index.php\/2014\/01\/02\/monte-carlo-bayesian-sensitivity-analysis\/","title":{"rendered":"Monte Carlo Bayesian Sensitivity Analysis"},"content":{"rendered":"\n\n\n
sens<\/td>\nR Documentation<\/td>\n<\/tr>\n<\/tbody>\n<\/table>\n

Monte Carlo Bayesian Sensitivity Analysis<\/h2>\n

Description<\/h3>\n

Fully Bayesian Monte Carlo sensitivity analysis scheme, based upon any of the regression models contained in the\u00a0tgp\u00a0package. 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.<\/p>\n

Usage<\/h3>\n
sens(X, Z, nn.lhs, model = btgp, ngrid = 100, span = 0.3,\n     BTE = c(3000,8000,10), rect = NULL, shape = NULL, mode = NULL,\n     ...)<\/pre>\n
Arguments<\/h3>\n\n\n\n\n\n\n\n\n\n\n\n\n\n
X<\/code><\/td>\ndata.frame<\/code>,\u00a0matrix<\/code>, or vector of inputs\u00a0X<\/code><\/td>\n<\/tr>\n
Z<\/code><\/td>\nVector of output responses\u00a0Z<\/code>\u00a0of length equal to the leading dimension (rows) of\u00a0X<\/code>, i.e.,\u00a0length(Z) == nrow(X)<\/code><\/td>\n<\/tr>\n
nn.lhs<\/code><\/td>\nSize of each Latin hypercube drawn for use in the Monte Carlo integration scheme. Total number of locations for prediction is\u00a0nn.lhs*(ncol(X)+2)<\/code><\/td>\n<\/tr>\n
model<\/code><\/td>\nEither the regression model used for prediction, or\u00a0NULL<\/code>. If\u00a0model=NULL<\/code>, then the function just returns the\u00a0sens.p<\/code>\u00a0vector of parameters to be passed with a regression model call. This can be used to perform sensitivity analysis through thepredict.tgp<\/code>\u00a0framework<\/td>\n<\/tr>\n
ngrid<\/code><\/td>\nThe number of grid points in each input dimension upon which main effects will be estimated.<\/td>\n<\/tr>\n
span<\/code><\/td>\nSmoothing parameter for main effects integration: the fraction of\u00a0nn.lhs<\/code>\u00a0points that will be included in a moving average window that is used to estimate main effects at the\u00a0ngrid<\/code>\u00a0locations in each input dimension.<\/td>\n<\/tr>\n
BTE<\/code><\/td>\n3-vector of Monte-Carlo parameters (B)urn in, (T)otal, and (E)very. Predictive samples are saved every E MCMC rounds starting at round B, stopping at T<\/td>\n<\/tr>\n
rect<\/code><\/td>\nRectangle 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\u00a0matrix<\/code>\u00a0or\u00a0data.frame<\/code>\u00a0withncol(rect) = 2<\/code>, and number of rows equal to the dimension of the domain. For 1-d data, a vector of length 2 is allowed. Defaults to the input data range (X<\/code>).<\/td>\n<\/tr>\n
shape<\/code><\/td>\nOptional 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\u00a0rect<\/code>. Only concave Beta distributions with\u00a0shape<\/code>\u00a0> 1 are supported. If unspecified, the uncertainty distribution is uniform over\u00a0rect<\/code>. The specification\u00a0shape[i]=0<\/code>\u00a0instructs\u00a0sens<\/code>\u00a0to treat the i’th dimension as a binary variable. In this case,\u00a0mode[i]<\/code>\u00a0is the probability parameter for a bernoulli uncertainty distribution, and we must also have\u00a0rect[i,]=c(0,1)<\/code>.<\/td>\n<\/tr>\n
mode<\/code><\/td>\nOptional 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\u00a0rect<\/code>. If\u00a0shape<\/code>\u00a0is specified, but this is not, then the scaled Beta distributions will be symmetric.<\/td>\n<\/tr>\n
...<\/code><\/td>\nExtra arguments to the\u00a0tgp\u00a0model<\/code>.<\/td>\n<\/tr>\n<\/tbody>\n<\/table>\n
Details<\/h3>\n
Saltelli (2002) describes a Latin Hypercube sampling based method for estimation of the ‘Sobal’ sensitivity indices:<\/p>\n
1st Order for input\u00a0i<\/i>,<\/p>\n
S(i) = var(E[f|x[i]])\/var(f),<\/i><\/p>\n
where\u00a0x[i]<\/i>\u00a0is the\u00a0i<\/i>-th input.<\/p>\n
Total Effect for input\u00a0i<\/i>,<\/p>\n
T(i) = E[var(f|x[-i])]\/var(f),<\/i><\/p>\n
where\u00a0x[-i]<\/i>\u00a0is all inputs except for the\u00a0i<\/i>-th.<\/p>\n
All moments are with respect to the appropriate marginals of the uncertainty distribution\u00a0U<\/i>\u00a0\u2013 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\u00a0nn.lhs<\/code>\u00a0is the sample size for the Monte Carlo estimate, this scheme requires\u00a0nn.lhs*(ncol(X)+2)<\/code>\u00a0function evaluations.<\/p>\n
The\u00a0sens<\/code>\u00a0function implements the method for unknown functions\u00a0f<\/i>, through prediction via one of the\u00a0tgp\u00a0regression models conditional on an observed set of\u00a0X<\/code>\u00a0locations. At each MCMC iteration of the\u00a0tgp\u00a0model fitting, the\u00a0nn.lhs*(ncol(X)+2)<\/code>\u00a0locations 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 non-parametric regression (a moving average).<\/p>\n
Please see\u00a0vignette(\"tgp2\")<\/code>\u00a0for a detailed illustration<\/p>\n
Value<\/h3>\n
The output is a\u00a0\"tgp\"<\/code>-class object. The details for\u00a0btgp<\/code>\u00a0contain a complete description of this output. The list entry that is relevance to sensitivity analysis is\u00a0sens<\/code>, which itself has entries:<\/p>\n\n\n\n\n\n\n\n\n\n
par<\/code><\/td>\nThis contains a\u00a0list<\/code>\u00a0of the input parameters used in the sensitivity analysis, as outlined above.<\/td>\n<\/tr>\n
Xgrid<\/code><\/td>\nA\u00a0matrix<\/code>\u00a0containing a grid in each input dimension (by column) over which the main effects are estimated.<\/td>\n<\/tr>\n
ZZ.mean<\/code><\/td>\nA\u00a0matrix<\/code>, where each column contains the mean main effects over the corresponding column of\u00a0sens.Xgrid<\/code>.<\/td>\n<\/tr>\n
ZZ.q1<\/code><\/td>\nA\u00a0matrix<\/code>, where each column contains the 5th percentile main effects over the corresponding column of\u00a0sens.Xgrid<\/code>.<\/td>\n<\/tr>\n
ZZ.q2<\/code><\/td>\nA\u00a0matrix<\/code>, where each column contains the 5th percentile main effects over the corresponding column of\u00a0sens.Xgrid<\/code>.<\/td>\n<\/tr>\n
S<\/code><\/td>\nA\u00a0matrix<\/code>, where each column contains the posterior sample of 1st order sensitivity indices for the corresponding input dimension.<\/td>\n<\/tr>\n
T<\/code><\/td>\nA\u00a0matrix<\/code>, where each column contains the posterior sample of total sensitivity indices for the corresponding input dimension.<\/td>\n<\/tr>\n<\/tbody>\n<\/table>\n
Note<\/h3>\n
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<\/code>) must increase with the dimensionality of the problem. The total sensitivity indices\u00a0T<\/i>\u00a0are forced non-negative, and if negative values occur it is necessary to increase\u00a0nn.lhs<\/code>. The\u00a0plot.tgp<\/code>\u00a0function replaces negative values with zero for illustration.<\/p>\n
Author(s)<\/h3>\n
Robert B. Gramacy,\u00a0rbgramacy@chicagobooth.edu<\/a>, and Matt Taddy,\u00a0taddy@chicagobooth.edu<\/a><\/p>\n
References<\/h3>\n
R.D. Morris, A. Kottas, M. Taddy, R. Furfaro, and B. Ganapol. (2009)\u00a0A statistical framework for the sensitivity analysis of radiative transfer models.<\/em>\u00a0IEEE Transactions on Geoscience and Remote Sensing, to appear.<\/p>\n
Saltelli, A. (2002)\u00a0Making best use of model evaluations to compute sensitivity indices.<\/em>\u00a0Computer Physics Communications, 145, 280-297.<\/p>\n
See Also<\/h3>\n
btgp<\/code>,\u00a0plot.tgp<\/code>,\u00a0predict.tgp<\/code>,\u00a0lhs<\/code><\/p>\n
Examples<\/h3>\n
# Take a look at the air quality in New York: Sensitivity of\n# ozone levels with respect to solar radiation, wind, and\n# temperature. See help(airquality) for details.  \nX <- airquality[,2:4]\nZ <- airquality$Ozone\n\n# Uncertainty distribution is the default: uniform over range(X)\n# There is missing data, which is removed automatically by tgp\n# range(X).\ns <- suppressWarnings(sens(X=X, Z=Z, nn.lhs=300, model=btgp,\n          ngrid=100, span=0.3, BTE=c(5000,10000,10)))\n\n# plot the results\nplot(s, layout=\"sens\", ylab=\"Ozone\", main=\"main effects\") \n\n# plot only the sensitivity indices\nplot(s, layout=\"sens\", ylab=\"Ozone\", maineff=FALSE) \n\n# plot only the main effects, side by side\nplot(s, layout=\"sens\", ylab=\"Ozone\", main=\"\", maineff=t(1:3))\n\n# build a 'sens.p' parameter vector for a data-dependent\n# informative uncertainty distribution.  For each variable,\n# the input distribution will be a scaled Beta with shape=2,\n# and mode equal to the data mean\nrect <- t(apply(X, 2, range, na.rm=TRUE))\nmode <- apply(X , 2, mean, na.rm=TRUE)\nshape <- rep(2,3)\n\n# plot a sample from the marginal uncertainty distribution.\nUdraw <- lhs(300, rect=rect, mode=mode, shape=shape)\npar(mfrow=c(1,3))\nfor(i in 1:3) hist(Udraw[,i], breaks=15,xlab=names(X)[i]) \n\n# build sens.p with the 'sens' function.\nsens.p <- suppressWarnings(sens(X=X,Z=Z,nn.lhs=300, model=NULL,\n               ngrid=100, rect=rect, shape=shape, mode=mode))\n\n# Use predict.tgp to quickly analyze with respect to this new\n# uncertainty distribution without re-running the MCMC, then\n# plot the results.\ns.new <- predict(s, BTE=c(1,1000,1), sens.p=sens.p, verb=1) \nplot(s.new, layout=\"sens\", ylab=\"Ozone\", main=\"main effects\")<\/pre>\n
Results<\/h3>\n
<\/pre>\n
 <\/p>\n\n\n\n
\"\"<\/a>\"\"<\/a>\"\"<\/a>\"\"<\/a>\"\"<\/a><\/td>\n<\/tr>\n<\/tbody>\n<\/table>\n","protected":false},"excerpt":{"rendered":"
sens R Documentation Monte Carlo Bayesian Sensitivity Analysis Description Fully Bayesian Monte Carlo sensitivity analysis scheme, based upon any of the regression models contained in… <\/p>\n","protected":false},"author":1,"featured_media":0,"comment_status":"open","ping_status":"open","sticky":false,"template":"","format":"standard","meta":{"footnotes":""},"categories":[20],"tags":[],"class_list":["post-416","post","type-post","status-publish","format-standard","hentry","category-r"],"_links":{"self":[{"href":"https:\/\/zhuoyao.net\/index.php\/wp-json\/wp\/v2\/posts\/416","targetHints":{"allow":["GET"]}}],"collection":[{"href":"https:\/\/zhuoyao.net\/index.php\/wp-json\/wp\/v2\/posts"}],"about":[{"href":"https:\/\/zhuoyao.net\/index.php\/wp-json\/wp\/v2\/types\/post"}],"author":[{"embeddable":true,"href":"https:\/\/zhuoyao.net\/index.php\/wp-json\/wp\/v2\/users\/1"}],"replies":[{"embeddable":true,"href":"https:\/\/zhuoyao.net\/index.php\/wp-json\/wp\/v2\/comments?post=416"}],"version-history":[{"count":0,"href":"https:\/\/zhuoyao.net\/index.php\/wp-json\/wp\/v2\/posts\/416\/revisions"}],"wp:attachment":[{"href":"https:\/\/zhuoyao.net\/index.php\/wp-json\/wp\/v2\/media?parent=416"}],"wp:term":[{"taxonomy":"category","embeddable":true,"href":"https:\/\/zhuoyao.net\/index.php\/wp-json\/wp\/v2\/categories?post=416"},{"taxonomy":"post_tag","embeddable":true,"href":"https:\/\/zhuoyao.net\/index.php\/wp-json\/wp\/v2\/tags?post=416"}],"curies":[{"name":"wp","href":"https:\/\/api.w.org\/{rel}","templated":true}]}}