Description Usage Arguments Value Author(s) References See Also Examples
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.
1 2 |
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 |
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 |
QuadType |
Type of quadrature. Options are |
QuadPoly |
The type of one-dimensional quadrature rule to be used. For |
QuadLevel |
Level of quadrature used in the approximation. |
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) |
Jordan Ko
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.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 | # 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)
|
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.