Bonjour,
je suis bloqué sur mon script d'analyse de sensibilité par la méthode de Sobol.
Au niveau de : #run du modèle a partir du plan d'experience Sobol-Saltelli
Message d'erreur : Error in run[[2]] : subscript out of bounds
Vous pourriez m'aider svp ?
Script :
#déclaration de la fonction
sz_d<-function(vw,ks,qw){
0.031*vw^(0.22)*(ks/qw)^(-0.17)
}
#création de deux "jeux" de paramètres LHS
n.tranches<-10
n.repet<-100
set.seed(123)
jeu_param1<-NULL
for(i in 1:n.repet)
{ ech.lhs<-randomLHS(n.tranches,3)
ech.lhs[,1]<-runif(n=ech.lhs[,1],min=1,max=100)
ech.lhs[,2]<-runif(n=ech.lhs[,2],min=10^(-7),max=10^(-2))
ech.lhs[,3]<-runif(n=ech.lhs[,3],min=1,max=10)
jeu_param1<-rbind(jeu_param1, ech.lhs)}
jeu_param2<-NULL
for(i in 1:n.repet)
{ ech.lhs<-randomLHS(n.tranches,3)
ech.lhs[,1]<-runif(n=ech.lhs[,1],min=1,max=100)
ech.lhs[,2]<-runif(n=ech.lhs[,2],min=10^(-7),max=10^(-2))
ech.lhs[,3]<-runif(n=ech.lhs[,3],min=1,max=10)
jeu_param2<-rbind(jeu_param2, ech.lhs)}
#création des graphiques
par(mfrow=c(1,1))
par(mar=c(2,2,0,0))
pairs(jeu_param1, main='jeu_param1')
pairs(jeu_param2, main='jeu_param2')
#mise en oeuvre de la methode de Sobol-Saltelli
install.packages("sensitivity_1.4-1.tar.gz", repos=NULL)
library(sensitivity)
#source('sz_d')
jeu_param1<-as.data.frame(jeu_param1)
jeu_param2<-as.data.frame(jeu_param2)
names(jeu_param1)<-c('vw','ks','qw')
names(jeu_param2)<-c('vw','ks','qw')
X <- sobol2002(model = NULL, X1 = jeu_param1, X2 = jeu_param2, nboot=100, conf = 0.95,)
#plan d'expérience
par(mfrow=c(1,1))
par(mar=c(2,2,0,0))
pairs(X$X, main='jeu_param')
#run du modèle a partir du plan d'experience Sobol-Saltelli
Y<-NULL
for (i in 1:nrow(X$X))
{
run<-sz_d(vw=X$X[i,1], ks=X$X[i,2], qw=X$X[i,3])
Y<-c(Y, sum(run[[2]][,2]*18/1000))
}
résultats<-Y
jeu_param<-X$X
#représentation des résultats
par(mfrow=c(1,3))
plot(jeu_param[,1],résultats, xlab='vw', ylab='d (m)')
plot(jeu_param[,2],résultats, xlab='ks', ylab='d (m)')
plot(jeu_param[,3],résultats, xlab='qw', ylab='d (m)')
#Si = V ar(E[Y=Xi])/V ar(Y )
par(mfrow=c(1,3))
plot(jeu_param[,1],résultats, xlab='vw', ylab='d (m)')
abline(v=runif(p=seq(0,1,by=1/n.tranches),min=1,max=100))
clas.vw<-cut(jeu_param[,1],breaks=runif(p=seq(0,1,by=1/n.tranches),
+ min=1,max=100))
mid.vw<-tapply(jeu_param[,1],clas.vw,'mean')
moy.vw<-tapply(résultats,clas.vw,'mean')
points(mid.vw, moy.vw, pch=19, col='red', cex=2)
plot(jeu_param[,2],résultats, xlab='ks', ylab='d (m)')
abline(v=runif(p=seq(0,1,by=1/n.tranches),min=10^(-7),max=10^(-3)))
clas.ks<-cut(jeu_param[,2],breaks=runif(p=seq(0,1,by=1/n.tranches),min=10^(-7),max=10^(-3)))
mid.ks<-tapply(jeu_param[,2],clas.ks,'mean')
moy.ks<-tapply(résultats,clas.ks,'mean')
points(mid.ks, moy.ks, pch=19, col='red', cex=2)
plot(jeu_param[,3],résultats, xlab='qw', ylab='d (m)')
abline(v=runif(p=seq(0,1,by=1/n.tranches),min=1,max=10))
clas.qw<-cut(jeu_param[,3],breaks=runif(p=seq(0,1,by=1/n.tranches),min=1,max=10))
mid.qw<-tapply(jeu_param[,3],clas.qw,'mean')
moy.qw<-tapply(résultats,clas.qw,'mean')
points(mid.qw, moy.qw, pch=19, col='red', cex=2)
var(moy.vw)/var(résultats)
var(moy.ks)/var(résultats)
var(moy.qw)/var(résultats)
#indices de sobol
Y<-as.numeric(Y)
tell(X,Y)
print(X)
Call:sobol2002(model = NULL, X1 = param1, X2 = param2, nboot = 100, conf = 0.95,)
par(mfrow=c(1,1))
plot(X, ylim=c(-1,2))
abline(h=c(0,1), col='red')