getSobol: Determine Sobol indices from PCE solutions

In GPC: Generalized Polynomial Chaos

Description

Computes the Sobol sensitivity indices from the coefficient of the polynomial chaos expansion.

Usage

getSobol(d,Index,Coeff,PhiIJ)

Arguments

d	Number of random inputs
Index	A (m x d) matrix that donates the polynomial order for each random variables in the canonical construction of PCE
Coeff	A m-tuple vector containing the PCE coefficients, ordered according to the canonical sequence of the multivariate polynomial expansion
PhiIJ	A m-tuple vector containing the variance of each.

Value

Names	The
Values

Author(s)

Jordan Ko

References

I. M. Sobol', 1993, Sensitivity estimate for non-linear mathematical models. Mathematical modelling and computational experiments 1, 407-414.

J. Ko, 2009, Applications of the generalized polynomial chaos to the numerical simulation of stochastic shear flows, Doctoral thesis, University of Paris VI.

See Also

getM,CreateQuadrature

Examples

### Ishigami function definition
Model <- function(x){
  a <- 3
  b <- 7
  res <- sin(pi*x[1,]) + a*sin(pi*x[2,])^2 + b*(pi*x[3,])^4*sin(pi*x[1,])
  return(res)
}

# Find exact solution for the Ishigami test function 
a <- 3
b <- 7
MeanExact <- a/2
VarExact <- a^2/8 + b*pi^4/5 + b^2*pi^8/18 + 1/2
SobolExact <- c(b*pi^4/5 + b^2*pi^8/50 + 1/2, a^2/8, 0, 0, 8*b^2*pi^8/225, 0, 0)

# random variable definition
d <- 3              # number of random variables
L <- 4              # quadrature level in each dimension. 
                    # could be anisotropic eg c(3,4,5) for full quadrature 
P <- L-1            # maximum polynomial expansion (cardinal order)
M <- getM(d,P)      # number of PCE terms
ParamDistrib <- NULL
# ParamDistrib<- list(beta=rep(0.0,d),alpha=rep(0.0,d))

# PCE definition
QuadType <- "FULL"              # type of quadrature
QuadPoly <- rep("LEGENDRE",d)   # polynomial to use
ExpPoly <- rep("LEGENDRE",d)    # polynomial to use

# QuadType <- "SPARSE"                                  # type of quadrature
# QuadPoly <- 'ClenshawCurtis'                          # polynomial to use
# ExpPoly <- rep("LEGENDRE",d)    # polynomial to use

Quadrature = CreateQuadrature(d,L,QuadPoly,ExpPoly,QuadType,ParamDistrib) # quadrature
 
# function sampling
y <- Model(Quadrature$QuadNodes) # Ishigami function d=3 

# generate PCE coefficients
PCE = generatePCEcoeff(M,Quadrature$QuadSize,y,Quadrature$PolyNodes,Quadrature$QuadWeights) # PCE

# Sobol' sensitivity analysis
Index <- indexCardinal(d,P)
Sobol = getSobol(d,Index,PCE$PCEcoeff,PCE$PhiIJ)

cat('PCE Sobol indices are ',Sobol$Values,'\n')
cat('Sobol absolte errors are ',abs(Sobol$Values-SobolExact),'\n')

Example output

Loading required package: randtoolbox
Loading required package: rngWELL
This is randtoolbox. For overview, type 'help("randtoolbox")'.
Loading required package: orthopolynom
Loading required package: polynom
Loading required package: ks
Loading required package: lars
Loaded lars 1.2

Quadrature level assumed to be isotropic.
PCE Sobol indices are  10622 0.709061 1.710349e-27 8.9161e-28 9377.792 7.511127e-29 1.331203e-30 
Sobol absolte errors are  1186.364 0.415939 1.710349e-27 8.9161e-28 7153.338 7.511127e-29 1.331203e-30 

GPC documentation built on May 30, 2017, 12:50 a.m.