rm(list=ls(all.names=TRUE)) library(lhs) library(DiceKriging) library(rgl) library(sensitivity) library(boot) library(numbers) #### Learning basis and test basis #### N_BA <- 150; ## Number of points in the learning basis d <- 3; ## Inputs space dimension S_BA <- randomLHS(N_BA,d) ## Construction of the learning basis ## S_BA : inputs for the learning basis ## Y_BA : corresponding outputs of the emulator S_BT <- randomLHS(100,3) ## Construction of the test basis ishigami <- function(X,coeff){ ## Analytical Function of Ishigami for scaled inputs => inputs in [0:1] ## Ishigami function : Y = sinX1+a*(sinX2).^2 + b*X3.^4*sinX1 for Xi ~ U[-pi;+pi] ## with a = coeff(1); b=coeff(2) ## Generally : coeff = [7 0.1]; X <- 2*pi*X-pi sin(X[,1]) + coeff[1]*sin(X[,2])^2 + coeff[2]*X[,3]^4*sin(X[,1]) } coeff_fx_test = c(7, 0.1) Y_BA <- ishigami(S_BA,coeff_fx_test) Y_BT <- ishigami(S_BT,coeff_fx_test) #### METAMODELING #### metamodel <- km(formula=~1,design=S_BA,response=Y_BA,covtype="matern5_2") plot(metamodel) # Affichage des resultats ## Computation of Q2 on the test basis Q2 <- function(y,yhat) 1-mean((y-yhat)^2)/var(y) prediction <- predict(metamodel, S_BT, 'UK', checkNames=FALSE)$mean Q2_BT<-Q2(Y_BT,prediction) print(paste('Q2 sur Base de test :',Q2_BT)) ## Computation of Q2 by cross validation prediction_LOO <- leaveOneOut.km(metamodel,'UK')$mean Q2_LOO<- Q2(Y_BA,prediction_LOO) print(paste('Q2 sur Base de test :',Q2_LOO)) #### Estimation of first-order and total Sobol' indices #### # Building of the gaussian process metamodel used to compute Sobol' indices f <- function(x) predict(metamodel,x,type='UK',checkNames=FALSE)$mean ## Computation of Sobol' indices by SobolEff Janon et al. (2014) n <- 1000 sample1 <- data.frame(X1=runif(n,0,1),X2=runif(n,0,1),X3=runif(n,0,1)) sample2 <- data.frame(X1=runif(n,0,1),X2=runif(n,0,1),X3=runif(n,0,1)) indices.sobol <- sobolEff(f,sample1,sample2,order=1) indices.sobol.totaux <- sobolEff(f,sample1,sample2,order=0) indices.sobol indices.sobol.totaux ## Computation of first- and closed second-order Sobol' indices with replicated designs indices.sobol.rep.o1 <- sobolroalhs(model = f, factors = 3, N = 1000, order = 1, nboot=100) print(indices.sobol.rep.o1) plot(indices.sobol.rep.o1) indices.sobol.rep.o2 <- sobolroalhs(model = f, factors = 3, N = 1000, order = 2, nboot=100) print(indices.sobol.rep.o2) plot(indices.sobol.rep.o2) ## Computation of Sobol' indices by Fast (Saltelli et al. [1999]),spectral approach, which requises regularity ## in the sense fast decrease of Fourier coefficients n_RBDFast <- 1000 M_RBDFast <- 4 indices.fast <- fast99(model=f,factors=d,n=n_RBDFast,M=M_RBDFast,q='qunif',q.arg=list(min=0,max=1)) #### Theoretical values #### ishigami.sobol.indices <- function(coeff){ D <- (coeff[1]^2)/8 + coeff[2]*(pi^4)/5 + coeff[2]^2*(pi^8)/18 + 0.5 D1 = coeff[2]*(pi^4)/5 + coeff[2]^2*(pi^8)/50 + 0.5 D2 = (coeff[1]^2)/8 D13 = coeff[2]^2*(pi^8)*(1/18-1/50) S1 <- St <- array(0,3) S1[1] = D1/D S1[2] = D2/D S1[3] = 0 St[1] = (D1+D13)/D St[2] = D2/D St[3] = D13/D list(S11=S1[1],S12=S1[2],S13=S1[3],St11=St[1],St12=St[2],St13=St[3]) } # Calcul des indices th?oriques indices.theo <- ishigami.sobol.indices(coeff_fx_test) # Comparaison des indices estim?s et th?oriques print('Theoretical indices :') indices.theo print('First-order and total indices estimated by SobolEff :') indices.sobol indices.sobol.totaux print('First-order and total indices estimated by fast1999 :') indices.fast print('First-order indices estimated by sobolroalhs :') indices.sobol.rep.o1 print('Closed second-order indices estimated by sobolroalhs :') indices.sobol.rep.o2 par(mfrow=c(1,2)) plot(indices.sobol) abline(h=indices.theo$S11) abline(h=indices.theo$S12) abline(h=indices.theo$S13) plot(indices.sobol.totaux) abline(h=indices.theo$St11) abline(h=indices.theo$St12) abline(h=indices.theo$St13) ## Play with N_BA (predictivity of the metamodel) and ## with n (the samples size for estimating DSobol' indices) # valeurs théoriques S3=0,S2 environ égal à 0.4424, S1 environ égal à 0.3139