marginals | R Documentation |
Estimate simple and double integrals for estimators of marginal mens and variances calculations for Gaussian process main effect and associated variance analysis
marginals(E, rhoarray, p, method='deterministic')
E |
A Gaussian process list as generated by |
rhoarray |
An array of |
p |
A vector of the indices over which the integral of squares are being estimated (see details) |
method |
The |
varanal
will estimate the emulator variances (more precisely, integral
of squared quantities !) associated with the variance of indices
complementary to p
(-p
in the notation of Oakley and Ohagan, 2004, p. 761),
and then take the expectation of this quantity over p
.
the less indices in "p", the more computationally intensive
given that all (-p) x (-p) combinations have to be summed up.
variances measures will then be obtained as follows,
total sensitivty index associated with, e.g., vars 1 and 2 (assuming 5 factors in total: varanal(E, rho, c(1,2,3,4,5)) - varanal(E, rho, c(3,4,5))
mean sensitivty index associated with the same variables varanal(E, rho, c(1,2) ) - varanal (E, rho, NULL )
Incidentially, varanal(E, rho, c(1,2,3,4,5)) - varanal(E, rho, NULL) will return the total variance.
A list with following elements:
simple_integrals |
Matrix with as many columns as levels supplied along |
simple_integrals |
Matrix with as many columns as levels supplied along |
Michel Crucifix
varanal
# generate data
X = as.matrix( expand.grid(seq(3), seq(3), seq(3)))
Y = apply (X, 1, function(x) { x[1] - 2*x[1]*x[2] + x[3] - 2} )
Y = Y + rnorm(length(Y), sd=0.1)
E = GP_C(X, Y, lambda=list(theta=c(1,1,1), nugget=0.1))
# define rho function : uniform distribution
n = 8
levels = list ( seq(0,3, l=n), seq(0,3, l=n), seq(0,3, l=n) )
rho = array(1, c(n,n,n))
rho = rho / sum(rho) # not necessary.
attr(rho, "levels") <- levels # this is necessary !
# mean sensitivity index associated with the index 2
E_2 = marginals(E, rho, c(2))
plot(E_2$simple_integrals['y',], type='l', main='Main effect of factor 2')
lines(E_2$simple_integrals['y',] +sqrt ( E_2$double_integrals['Sp',] ),
lty=2, col='blue')
# variances associated with Gaussian process variance
lines(E_2$simple_integrals['y',] -sqrt ( E_2$double_integrals['Sp',] ),
lty=2, col='blue')
# variances associated with process mean :
lines(E_2$simple_integrals['y',] +sqrt (
E_2$simple_integrals['yyt',] - E_2$simple_integrals['y',]^2 ),
lty=2, col='red')
lines(E_2$simple_integrals['y',] -sqrt (
E_2$simple_integrals['yyt',] - E_2$simple_integrals['y',]^2 ),
lty=2, col='red')
legend('bottomleft', c('sd associated with GP variance', 'sd associated with process mean'),
lty=2, col=c('red','blue'))
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.