GPCE.quad: Solve polynomial chaos expansion with numerical quadrature

In GPC: Generalized Polynomial Chaos

Description

The GPCE.quad function implements a polynomial chaos expansion of a given model or an external model. The strategy for the expansion of the model into a polynomial chaos basis is the Gauss quadrature method where the Gauss quadrature points are used to estimate the integrales corresponding to the coefficients of the expansion. A statistical and a global sensitivity analysis of the model is then carried out.

Usage

GPCE.quad(InputDim,PCSpace,InputDistrib,ParamDistrib,
Model,ModelParam,Output,DesignInput,p,ExpPoly,QuadType,QuadPoly,QuadLevel)

Arguments

InputDim	Number of input random variables
PCSpace	To decide if it should be removed after discussion with Miguel..
InputDistrib	Distribution of the input. Options are Gaussian, Uniform, Beta, Gamma
ParamDistrib	Shape parameters of the random variables.
Model	A function defining the model to analyze or NULL if the model is external
ModelParam	Optional parameters for function evaluation.
Output	The results of the model evaluation at the quadrature points.
DesignInput	To decide if it should be removed after discussion with Miguel.
p	The order of the polynomial chaos expansion.

order)

ExpPoly	The polynomial used in the expansion. Options are Hermite, Legendre, Jacobi, Laguerre
QuadType	Type of quadrature. Options are SPARSE or FULL
QuadPoly	The type of one-dimensional quadrature rule to be used. For SPARSE, one can use ClenshawCurtis or Fejer. For FULL, one could choose Hermite,Legendre,Jacobi or Laguerre
QuadLevel	Level of quadrature used in the approximation.

Value

Designs	A list containing the quadrature size, n, a vector of the n quadrature weights, and a n * d matrix of the quadrature points
Output	Vector of the model output
PCEcoeff	Matrix of the kept sparse polynomial basis. TruncSet_ij is the jth polynomial degree associated to the ith variable
Moments	A list containing the fourth first moments of the output: mean, variance, standard deviation, skweness and kurtosis
Sensitivity	A list containing the sobol sensitivity indices (Values) and the normalized Sobol nominal and total sensitivity indices (S and ST)

Author(s)

Jordan Ko

References

J. Ko, D. Lucor and P. Sagaut, 2008, On Sensitivity of Two-Dimensional Spatially Developing Mixing Layer With Respect to Uncertain Inflow Conditions, Physics of Fluids, 20(7), 07710201-07710220.

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

J. Ko, D. Lucor and P. Sagaut, 2011, Effects of base flow uncertainty on Couette flow stability, Computers and Fluids, 43(1), 82-89.

See Also

tell.GPCE.quad

Examples

# CASE 1: Validating Legendre PCE with Ishigami function definition
rm(list = setdiff(ls(), lsf.str()))
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)
}

# Determine the 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)

d = 3
l = 4
ResultObject=GPCE.quad(Model=Model,InputDim=d,PCSpace="Uniform",InputDistrib=rep('Uniform',d),
                       DesignInput=NULL,p=c(l-1),ExpPoly=rep("LEGENDRE",d),QuadType=c("FULL"),
                       QuadPoly=rep("LEGENDRE",d),QuadLevel=c(l),ParamDistrib=NULL,Output=NULL)
names(ResultObject)
cat('Mean and variance normalized absolute errors are', 
  abs(ResultObject$Moments$PCEMean-MeanExact)/MeanExact,
  'and', abs(ResultObject$Moments$PCEVar-VarExact)/VarExact,'\n')
cat('PCE Sobol indices are ',ResultObject$Sensitivity$Values,'\n')
cat('Sobol absolte errors are ',abs(ResultObject$Sensitivity$Values-SobolExact),'\n')

# CASE 2: Validating Hermite PCE with orthogonal Hermite functions
### Model is a R function as a sum of multivariate Hermite polynomials 

Model <- function(x,param){
    d <- param$d
    p <- param$p
    PCETrue <- param$PCETrue
    
    n <- dim(x)[2]
    index <- indexCardinal(d,p)
    PHerm <- hermite.he.polynomials(p, normalized=FALSE)
    y <- rep(0,n)
    
    for (nn in seq(1,n)){
        for (mm in seq(1,getM(d,p))){
            tmp <- 1;
            for (dd in seq(1,d))
            {
            tmp = tmp * unlist(polynomial.values(PHerm[index[dd,mm]+1],x[dd,nn]))
            }     
            y[nn] = y[nn] + PCETrue[mm]*tmp
        }
    }
    return(y)
}

## Problem definition
d = 2;          # random dimension 
l = 3;          # quadrature level
p = l - 1;      # polynomial order of expansion
m = getM(d,p);  # size of polynomial expansion

## Model definition
ModelParam <- NULL
ModelParam$d <- d
ModelParam$p <- p
ModelParam$PCETrue <- sample(seq(1,m),m,replace = FALSE)

## CASE 1: The model is directly evaluated from the GPCE.quad function
ResultObject=GPCE.quad(InputDim=d,PCSpace="Normal",InputDistrib=rep('Gaussian',d),
                       DesignInput=NULL,p=c(p),ExpPoly=rep("HERMITE",d),QuadType=c("FULL"),
                       QuadPoly=rep("HERMITE",d),QuadLevel=c(l),ParamDistrib=NULL,Output=NULL,
                       Model=Model,ModelParam=ModelParam)
cat("The exact PCE coefficients are: \n")
cat(ModelParam$PCETrue,"\n")
cat("The estimated PCE coefficients are: \n")
cat(ResultObject$PCEcoeff,"\n")

## CASE 2: Model is evaluated separately from the GPCE.quad function
# First, the quadrature points are determined from the GPCE.quad function
NoModelResult=GPCE.quad(InputDim=d,PCSpace="Normal",InputDistrib=rep('Gaussian',d),
                       DesignInput=NULL,p=c(p),ExpPoly=rep("HERMITE",d),QuadType=c("FULL"),
                       QuadPoly=rep("HERMITE",d),QuadLevel=c(l),ParamDistrib=NULL,Output=NULL)
cat("The quadrature points can be determined from the Design variable of the output below: \n")
cat(names(NoModelResult),"\n")

# Second, the model is evalauted at the  quadrature points and stored in Output
Output <- Model(NoModelResult$Design$QuadNodes,ModelParam)

# Third, the model output is passed back to GPCE.quad, along with DesignInput and Output
cat("After Design$QuadNodes are evaluated and stored in Output, 
the results is passed back to GPCE.quad:\n")
NoModelResult=GPCE.quad(InputDim=d,PCSpace = "Normal",InputDistrib=rep('Gaussian',d),
                        DesignInput=NoModelResult$Design$QuadNodes,p=c(p),
                        ExpPoly=rep("HERMITE",d),QuadType=c("FULL"),
                        QuadPoly=rep("HERMITE",d),QuadLevel=c(l),
                        ParamDistrib=NULL,Output=Output)
cat("The exact PCE coefficients are:\n")
cat(ModelParam$PCETrue,"\n")
cat("The estimated PCE coefficients are:\n")
cat(NoModelResult$PCEcoeff,"\n")

### Run the algorithm with the lar regression method
# ResultObject=GPCE.lar(Model=Model, PCSpace="Gaussian", InputDim=d, InputDistrib=rep("Gaussian",d))


# CASE 3: Validating Jacobi PCE with orthogonal Jacobi functions 
rm(list = setdiff(ls(), lsf.str()))

# Multivariate Jacobi polynomial test
Model <- function(x){
    a = 2;
    b = 3;
    res <- 2 + 
        3*(2*(a+1)+(a+b+2)*(x[1,]-1))/2 + 
        5*(2*(a+1)+(a+b+2)*(x[2,]-1))/2 + 
        7*(2*(a+1)+(a+b+2)*(x[3,]-1))/2 + 
        11.*(4*(a+1)*(a+2)+4*(a+b+3)*(a+2)*(x[1,]-1)+(a+b+3)*(a+b+4)*(x[1,]-1)^2)/8 + 
        13*(2*(a+1)+(a+b+2)*(x[1,]-1))*(2*(a+1)+(a+b+2)*(x[2,]-1))/4 + 
        17*(2*(a+1)+(a+b+2)*(x[1,]-1))*(2*(a+1)+(a+b+2)*(x[3,]-1))/4 +         
        19.*(4*(a+1)*(a+2)+4*(a+b+3)*(a+2)*(x[2,]-1)+(a+b+3)*(a+b+4)*(x[2,]-1)^2)/8 + 
        23*(2*(a+1)+(a+b+2)*(x[2,]-1))*(2*(a+1)+(a+b+2)*(x[3,]-1))/4 + 
        29.*(4*(a+1)*(a+2)+4*(a+b+3)*(a+2)*(x[3,]-1)+(a+b+3)*(a+b+4)*(x[3,]-1)^2)/8
    return(res)
}

# x is a random variabls following the beta distribution with the pdf
# f(x,a,b) = (1-x)^a(1+x)^b/2^(a+b+1)/Beta(a+1,b+1), where Beta is the beta function 

d = 3
l = 3
ParamDistrib <- NULL
ParamDistrib$alpha <- rep(2,d)  # 1st shape parameters for the random variables
ParamDistrib$beta <- rep(3,d)   # 2nd shape parameters for the random variables

ResultObject=GPCE.quad(Model=Model,InputDim=d,PCSpace="Beta",InputDistrib=rep('Beta',d),
                       DesignInput=NULL,p=c(l-1),ExpPoly=rep("JACOBI",d),QuadType=c("FULL"),
                       QuadPoly=rep("JACOBI",d),QuadLevel=c(l),ParamDistrib,Output=NULL)

cat('Estimated PCE coefficients are', ResultObject$PCEcoeff) 

# CASE 4: Validating Laguerre PCE with orthogonal Laguerre functions 

rm(list = setdiff(ls(), lsf.str()))

# Multivariate Laguerre polynomial test
Model <- function(x){
    a <- 2
    res <- 2 + 
      3*(-x[1,]+a+1) +
      5*(-x[2,]+a+1) + 
      7*(-x[3,]+a+1) +
      11*(x[1,]^2/2-(a+2)*x[1,]+(a+2)*(a+1)/2) + 
      13*(-x[1,]+a+1)*(-x[2,]+a+1) + 
      17*(-x[1,]+a+1)*(-x[3,]+a+1) + 
      19*(x[2,]^2/2-(a+2)*x[2,]+(a+2)*(a+1)/2) + 
      21*(-x[2,]+a+1)*(-x[3,]+a+1) +
      29*(x[3,]^2/2-(a+2)*x[3,]+(a+2)*(a+1)/2) 
    return(res)
}

# x is a random variabls following the gamma distribution with scale parameter theta = 1 and shape 
# parameter k = a + 1 giving the pdf f(x,a) = x^(k-1)*exp(-x)/Gamma(k) = x^a*exp(-x)/Gamma(a+1), 
# where Gamma is the gamma function 

a = 2
d = 3
l = 3

ParamDistrib = NULL
ParamDistrib$alpha = rep(a,d)

ResultObject=GPCE.quad(Model=Model,InputDim=d,PCSpace="Gamma",InputDistrib=rep('Gamma',d),
                       DesignInput=NULL,p=c(l-1),ExpPoly=rep("LAGUERRE",d),QuadType=c("FULL"),
                       QuadPoly=rep("LAGUERRE",d),QuadLevel=c(l),ParamDistrib,Output=NULL)

cat('Estimated PCE coefficients are', ResultObject$PCEcoeff) 

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