Determine Sobol indices from PCE solutions

Description

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

Usage

getSobol(d,Index,Coeff,PhiIJ)

Arguments

Value

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')

Results

R version 3.3.1 (2016-06-21) -- "Bug in Your Hair"
Copyright (C) 2016 The R Foundation for Statistical Computing
Platform: x86_64-pc-linux-gnu (64-bit)

R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.

R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.

Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.

> library(GPC)
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: KernSmooth
KernSmooth 2.23 loaded
Copyright M. P. Wand 1997-2009
Loading required package: misc3d
Loading required package: mvtnorm
Loading required package: rgl
Loading required package: lars
Loaded lars 1.2

> png(filename="/home/ddbj/snapshot/RGM3/R_CC/result/GPC/getSobol.Rd_%03d_medium.png", width=480, height=480)
> ### Name: getSobol
> ### Title: Determine Sobol indices from PCE solutions
> ### Aliases: getSobol
> 
> ### ** 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
Quadrature level assumed to be isotropic.
>  
> # 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')
PCE Sobol indices are  10622 0.709061 1.710349e-27 8.9161e-28 9377.792 7.511127e-29 1.331203e-30 
> cat('Sobol absolte errors are ',abs(Sobol$Values-SobolExact),'\n')
Sobol absolte errors are  1186.364 0.415939 1.710349e-27 8.9161e-28 7153.338 7.511127e-29 1.331203e-30 
> 
> 
> 
> 
> 
> dev.off()
null device 
          1 
>