Analyse de sensibilité Sobol

Postez ici vos questions, réponses, commentaires ou suggestions - Les sujets seront ultérieurement répartis dans les archives par les modérateurs

Modérateur : Groupe des modérateurs

Rémi Scotto di Liguori
Messages : 1
Enregistré le : 27 Juin 2019, 07:49

Analyse de sensibilité Sobol

Messagepar Rémi Scotto di Liguori » 27 Juin 2019, 09:04

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')

Mickael Canouil
Messages : 1315
Enregistré le : 04 Avr 2011, 08:53
Contact :

Re: Analyse de sensibilité Sobol

Messagepar Mickael Canouil » 27 Juin 2019, 09:53

Bonjour,

pour partager du code, c'est mieux d'utiliser les balises "Code" (outon dans l’éditeur).

Votre erreur est run[[2]], il n'y a pas de second élément (si tant est que run soit une liste).

Code : Tout sélectionner

x <- list(1, 2)
x[[3]]
#> Error in x[[3]]: subscript out of bounds


Allez donc à la ligne contenant "run[[2]]" et exécuter votre code pas-à-pas pour identifier là où vous vous êtes trompé.

PS: c'est osé d'utiliser des accents pour des noms de variables "résultats".

PS2: pour une meilleure lisibilité de votre code (je ne l'ai pas lu), je vous recommande par exemple https://style.tidyverse.org/

Cordialement,
Mickaël
mickael.canouil.fr | rlille.fr


Retourner vers « Questions en cours »

Qui est en ligne

Utilisateurs parcourant ce forum : Google [Bot] et 1 invité