# NOT RUN {
# Model: 3D function
distribution = list()
for (i in 1:3) distribution[[i]]=list("norm",c(0,1))
# Monte Carlo sampling
N = 10000
X = matrix(0,ncol=3,nrow=N)
for(i in 1:3) X[,i] = rnorm(N,0,1)
Y = 2 * X[,1] + X[,2] + X[,3]/2
q95 = quantile(Y,0.95)
nboot=200
# sensitivity indices with perturbation of the mean
v_delta = seq(-1,1,1/10)
toto = PLIquantile(0.95,X,Y,deltasvector=v_delta,
InputDistributions=distribution,type="MOY",samedelta=TRUE,
percentage=FALSE,nboot=nboot)
# Plotting the PLI
par(mar=c(4,5,1,1))
plot(v_delta,toto$PLI[,2],ylim=c(-1.5,1.5),xlab=expression(delta),
ylab=expression(hat(PLI[i*delta])),pch=19,cex=1.5)
points(v_delta,toto$PLI[,1],col="darkgreen",pch=15,cex=1.5)
points(v_delta,toto$PLI[,3],col="red",pch=17,cex=1.5)
lines(v_delta,toto$PLICIinf[,2],col="black")
lines(v_delta,toto$PLICIsup[,2],col="black")
lines(v_delta,toto$PLICIinf[,1],col="darkgreen")
lines(v_delta,toto$PLICIsup[,1],col="darkgreen")
lines(v_delta,toto$PLICIinf[,3],col="red")
lines(v_delta,toto$PLICIsup[,3],col="red")
abline(h=0,lty=2)
legend(0.8,1.5,legend=c("X1","X2","X3"),
col=c("darkgreen","black","red"),pch=c(15,19,17),cex=1.5)
# Plotting the perturbed quantiles
par(mar=c(4,5,1,1))
plot(v_delta,toto$quantile[,2],ylim=c(1.5,6.5),xlab=expression(delta),
ylab=expression(hat(q[i*delta])),pch=19,cex=1.5)
points(v_delta,toto$quantile[,1],col="darkgreen",pch=15,cex=1.5)
points(v_delta,toto$quantile[,3],col="red",pch=17,cex=1.5)
lines(v_delta,toto$quantileCIinf[,2],col="black")
lines(v_delta,toto$quantileCIsup[,2],col="black")
lines(v_delta,toto$quantileCIinf[,1],col="darkgreen")
lines(v_delta,toto$quantileCIsup[,1],col="darkgreen")
lines(v_delta,toto$quantileCIinf[,3],col="red")
lines(v_delta,toto$quantileCIsup[,3],col="red")
abline(h=q95,lty=2)
legend(0.5,2.4,legend=c("X1","X2","X3"),
col=c("darkgreen","black","red"),pch=c(15,19,17),cex=1.5)
# Plotting the PLI in percentage with refined confidence intervals
toto = PLIquantile(0.95,X,Y,deltasvector=v_delta,
InputDistributions=distribution,type="MOY",samedelta=TRUE,
percentage=TRUE,nboot=nboot,bootsample=FALSE)
par(mar=c(4,5,1,1))
plot(v_delta,toto$PLI[,2],ylim=c(-0.6,0.6),xlab=expression(delta),
ylab=expression(hat(PLI[i*delta])),pch=19,cex=1.5)
points(v_delta,toto$PLI[,1],col="darkgreen",pch=15,cex=1.5)
points(v_delta,toto$PLI[,3],col="red",pch=17,cex=1.5)
lines(v_delta,toto$PLICIinf[,2],col="black")
lines(v_delta,toto$PLICIsup[,2],col="black")
lines(v_delta,toto$PLICIinf[,1],col="darkgreen")
lines(v_delta,toto$PLICIsup[,1],col="darkgreen")
lines(v_delta,toto$PLICIinf[,3],col="red")
lines(v_delta,toto$PLICIsup[,3],col="red")
abline(h=0,lty=2)
legend(0.,1.,legend=c("X1","X2","X3"),
col=c("darkgreen","black","red"),pch=c(15,19,17),cex=1.5)
# }
Run the code above in your browser using DataLab