varanal: Essential calculations for Gaussian process variance analysis

In mcrucifix/gp: Gaussian process calibration for emulation

Essential calculations for Gaussian process variance analysis

Description

Essential calculations for Gaussian process variance analysis based on a determinstic integral.

Usage

varanal(E, rhoarray, p, method='deterministic')

Arguments

E	A Gaussian process list as generated by GP_C

rhoarray	An array of n1 x n2 x n3 ... with as many dimensions as the number of design factors used to generate the Gaussian process. The array represent the discretized density function. The densitity must be associated with a list of levels correspnding to the mid-points of the distribution. See examples.
p	A vector of the indices over which the integral of squares are being estimated (see details)
method	The determistic method assuming a sum over all components of rhoarray is currently the only one supported.

Details

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.

Value

Vector with following indices

ev	Estimate of -p x -p integral of GP covariance
sv	Estimate of -p x -p integral of squared emulator mean
svm	Contribution to sv associated with regression (constant, linear, quandratic)
svmg	Contribution to sv associated with covariance between regression and stochastic component
svg	Contribution to sv associated with stochastic component

Author(s)

Michel Crucifix

References

Jeremy E. Oakley and Anthony O'Hagan, Probabilistic sensitivity analysis of complex models: a Bayesian approach, Journal of the Royal Statistical Society: Series B (Statistical Methodology), 66, 751–769 2004

See Also

GP_P

Examples

# 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] } )
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_F = varanal(E, rho, c(1,2,3))
E_13 = varanal(E, rho,  c(1,3))
E_2 = varanal(E, rho, c(2))
E_0 = varanal(E, rho,   NULL)

TotalVariance = (E_F-E_0)['ev'] + (E_F-E_0)['sv'] 

print ('Mean Sens. Index. associated with fact.  2')
print (E_2-E_0)/TotalVariance

print ('Total Sens. Index. associated with fact.  2')
print (E_F-E_13)/TotalVariance

# total sensivity index associated with the second index
print ("observe that in this particular case the total sensitvitiy 
index of 2 is very small due to its non-linear combination with index 1")

mcrucifix/gp documentation built on July 29, 2023, 8:58 p.m.