Boucle programmation R

[inviteb8843a7e]

Bonjour,
je souhaiterais compiler dans une matrice les résultats de cette fonction selon une série de valeurs pris par mes paramètres et (le paramètre dans la fonction est à ignorer)


Étant novice sur R, je ne vois pas ce qui cloche dans mon programme:

# Mes vecteurs contenant les différentes valeurs de ma fonction f():
a <- seq(from = 0.001, to = 0.01, by = 0.0001)
b <- seq(from = 0.1, to = 0.5, by = 0.01)

ascalar <- 0.001
bscalar <- 0.1

# Ma matrice censée contenir mes résultats:
result <- matrix(data = 0, nrow = length(a), ncol = length(b))

Integrand <- function(x) {
      exp(-(ascalar/bscalar) * exp(bscalar*x))
}

ey <- function(y) {

# Balayage sur mes vecteurs a et b
          for (i in 1:length(a)) {
          ascalar <- a[i]
            for (j in 1:length(b)) {
          bscalar <- b[j]

              result[i,j] <- exp((a[i]/b[j])*exp(b[j]*y))*as.numeric(integrate(Integrand,y,Inf)["value"][1])
          }
         }
        }


test <- ey(0)

Au passage, mes résultats sont-ils censés être contenu dans la matrice "result" ou dans l'objet test ??
Je suis perdu
-----

[invitec5eb4b89]

Bonjour,

Pour un novice, je trouve ton code pas mal du tout !
Il y a cependant des variables globales qui sont utilisées en local par tes fonctions, ce qui me dérange un peu.

Pour répondre à ta question, le problème était que tu avais oublié le "return" dans la fonction ey.

Voici comment je propose de modifier ton code :

# Mes vecteurs contenant les différentes valeurs de ma fonction f():
A <- seq(from = 0.001, to = 0.01, by = 0.0001)
B <- seq(from = 0.1, to = 0.5, by = 0.01)

# Déclaration de fonctions
Integrand <- function(x,a,b) exp(-(a/b) * exp(b*x))
temp <- function(y,a,b) exp((a/b)*exp(b*y))*as.numeric(integrate(Integrand,lower=y,upper=Inf,a=a,b=b)["value"][1])

ey <- function(y,veca,vecb) {
  # Matrice censée contenir mes résultats:
  result <- matrix(data = 0, nrow = length(veca), ncol = length(vecb))
  # Balayage sur les vecteurs veca et vecb
  for (i in 1:length(a)) {
   for (j in 1:length(b)) {
    result[i,j] <- temp(y,veca[i],vecb[j])
   }
  }
  # Ajout : pour pouvoir obtenir la matrice à la fin des boucles
  return(result)
}

test <- ey(0,A,B)

J'espère que cela répond à tes questions.
V.

[inviteb8843a7e]

Merci infiniment HigginsVincent,

mieux que ce que " j'ai demandé ", j'ai grâce à ton programme 4 dimensions sur lequelles je peux jouer, à savoir:
- y
- a
- b
- et ma fonction ey elle même ,

cela me permettra sur quelques graphique de voir comment ma fonction évolue en fonction des 3 variables (y,a,b).
Petit détail, il a fallut au passage dans la double boucle écrire length(A) et length(B) au lieu de length(a) et length(b). Sans ça tout colle parfaitement.

Encore merci.

Au passage, si j'ai bien compris donc, ma commande return{result}, une fois la fonction ey exécuté, ne produira jamais une matrice du nom de "result" contenant mais résultat comme je le penser au départ, mais plutôt, l'objet du nom de test <-ey(O,A,B) . C'est bien ça ?
Je veux dire par là que je m'attendais au début à avoir mes résultats dans l'objet "result", alors qu'ils sont dans "test".

[invitec5eb4b89]

Rah zut, effectivement, j'ai mal défini la fin des boucles.
Mais il ne faut pas utiliser length(A), car A est un vecteur qui est défini en dehors de la fonction ey, il faut donc à tout pris éviter de l'utiliser à l'intérieur de la fonction ey. Imaginons par exemple que tu décides de renommer A en la variable toto, ta fonction ey sera devenu inutilisable ! Pour éviter ce genre d'erreur, il vaut mieux déclarer les fonctions avant les variables, ce que j'ai oublié de faire

Bref, voici la nouvelle version :

# Déclaration de fonctions
Integrand <- function(x,a,b) exp(-(a/b) * exp(b*x))
temp <- function(y,a,b) exp((a/b)*exp(b*y))*as.numeric(integrate(Integrand,lower=y,upper=Inf,a=a,b=b)["value"][1])

ey <- function(y,veca,vecb) {
  # Matrice censée contenir mes résultats:
  result <- matrix(data = 0, nrow = length(veca), ncol = length(vecb))
  # Balayage sur les vecteurs veca et vecb
  for (i in seq_along(veca)) 
   for (j in seq_along(vecb)) result[i,j] <- temp(y,veca[i],vecb[j])
  # Ajout : pour pouvoir obtenir la matrice à la fin des boucles
  return(result)
}

# Mes vecteurs contenant les différentes valeurs de ma fonction f():
A <- seq(from = 0.001, to = 0.01, by = 0.0001)
B <- seq(from = 0.1, to = 0.5, by = 0.01)

test <- ey(0,A,B)

J'avoue que j'intuite seulement que la variable "result" n'est pas déclarée à la fin de l'exécution de la fonction car l'utilisateur peut tout à fait vouloir exécuter cette fonction plusieurs fois avec des paramètres différents. Retourner "result" à chaque fois aurait comme résultat assez agaçant d'écraser le résultat précédent...

Bon courage !
V.


ps : je suis curieux de voir ton graphe, si jamais tu peux le poster ici, je serais intéressé par voir le code associé

[A voir en vidéo sur Futura]

[inviteb8843a7e]

Merci pour ces précieux conseils, je comprends mieux où tu veux en venir.
Pour l'aspect graphique, voici ce que donnerais la fonction, pour y=0.

Ma fonction est en fait l'espérance de vie à la naissance (y=age=0) étant donnée un
certain profil de mortalité ("à la Gompertz" dans le jargon).

[inviteb8843a7e]

Une dernière question, si je peux me permettre, car en fait je pourrais
éventuellement généraliser un peu plus "ma procédure", en faisant varier y (l'âge) tout comme j'ai fait varier a et b dans ma fonction.
Dans ce cas j'obtiendrais autant de matrice de taille (length(veca),length(vecb)) que j'ai de y (y variant de 0 à 100 ans par exemple).
Je pourrais donc dans le même graph en 3-D "ploter" une fonction pour chaque valeur de y (i.e pour chaque âge, espérance de vide à 0 an, 1 an, 2 ans,...)....cooooooollll

Mais comment m'y prendre ?

[invitec5eb4b89]

Je ne suis pas sûr de moi, mais il me semble que tu peux essayer d'utiliser la fonction "movie3d" de la librairie "rgl". J'ai essayé le code suivant :

# inititialiser
rm(list=ls())
# charger librairies
require(rgl)
# Déclaration de fonctions
Integrand <- function(x,a,b) exp(-(a/b) * exp(b*x))
temp <- function(y,a,b) exp((a/b)*exp(b*y))*as.numeric(integrate(Integrand,lower=y,upper=Inf,a=a,b=b)["value"][1])
ey <- function(y,veca,vecb) {
  result <- matrix(data = 0, nrow = length(veca), ncol = length(vecb))
  for (i in seq_along(veca)) 
   for (j in seq_along(vecb)) result[i,j] <- temp(y,veca[i],vecb[j])
  return(result)
}
 
#### représenter les graphes en 3d pour plusieurs valeurs de y
moove <- function(y) ey(y,seq(from = 0.001, to = 0.01, by = 0.0001),seq(from = 0.1, to = 0.5, by = 0.01))
A <- seq(from = 0.001, to = 0.01, by = 0.0001)
B <- seq(from = 0.1, to = 0.5, by = 0.01)
plot3d( A,B,ey(0,A,B) )
movie3d( moove,duration=4)

Pour être honnête, je ne sais pas vraiment si ca donne un résultat très satisfaisant...

[inviteb8843a7e]

J'aime plutôt bien ce que ça donne , je m'approche de plus en plus de ce que j'ai envie d'obtenir finalement. Sinon, hormis le fait que j'aimerais avoir quelques indications sur les lignes de niveau: savoir par exemple laquelle réfère à ey0, ey1,...autrement je vais voir si il n'y a pas moyen de dessiner les plans au lieu d'avoir des points (comme le graph que j'ai mis en p.j). ça devrait finalement donner un plus joli graphique avec plusieurs surfaces à l'intérieur...
Je vais essayé dans apprendre plus sur cette librairie et savoir ce qu'elle a dans le ventre.

[invitec5eb4b89]

Pardon, je me suis complètement planté dans l'utilisation de la fonction plot3d. Du coup j'ai modifié la fonction ey pour qu'elle soit plus adaptée à la représentation en 3 dimensions :

# inititialiser
rm(list=ls())
# charger librairies
require(rgl)
# Déclaration de fonctions
Integrand <- function(x,a,b) exp(-(a/b) * exp(b*x))
temp <- function(y,a,b) exp((a/b)*exp(b*y))*as.numeric(integrate(Integrand,lower=y,upper=Inf,a=a,b=b)["value"][1])

ey <- function(y,veca,vecb) {
  result <- matrix(NA,length(veca)*length(vecb),3)
  colnames(result) <- c("a","b","ey")
  k <- 1
  for (i in seq_along(veca))
   for (j in seq_along(vecb)) { result[k,] <- c(veca[i],vecb[j],temp(y,veca[i],vecb[j])) ; k <- k+1 }
  return(result)
}

# Mes vecteurs contenant les différentes valeurs de ma fonction f():
A <- seq(0.001,0.01, by = 0.001)
B <- seq(0.1,0.5, by = 0.1)

persp3d(ey(0,A,B),col=1)
persp3d(ey(1,A,B),col=2,add=TRUE)
persp3d(ey(2,A,B),col=3,add=TRUE)
persp3d(ey(3,A,B),col=4,add=TRUE)

Malheureusement, c'est pas très beau, il faudrait ajouter de la transparence...

[inviteb8843a7e]

Comme je te l'ai dit je suis novice, je vais donc relire ton programme et comprendre un peu ce qu'il en est.

Sinon, au passage, pour le graph je le trouve très beau (même si je suis un peu surpris de la forme qu'il prend=> je vérifierai ça en relisant ton programme). En fait si tu ne vois pas très bien ce qui s'y passe, c'est parce que tu as utilisé des valeurs trop rapprochées l'une de l'autre pour "y" ("y" qui est, je te le rappelle, ma variable âge dans ma fonction de départ). Autrement dit, tu as "ploté" l'espérance de vie à 0 an, 1 an, 2 ans, 3 ans, c'est donc normalement de ne pas y voir grand chose puisque l'espérance de vie ne change pas beaucoup d'une année à une autre, du moins pour ce modèle. C'est pour ça qu'il est difficile de faire la distinction entre les différentes lignes de niveaux; des valeurs plus espacées l'une de l'autre aide à y voir plus clair en fait (par, exple y=0, y=10, y=20, y=30).

Merci beaucoup pour "les coups de pouce".
Je vais relire ton prog.

[inviteb8843a7e]

Encore moi HigginsVincent ,

j'ai finalement considéré le programme suivant car il me semble y avoir une erreur dans le dernier que tu m'a concocté (la fonction "ploté" semble ne pas être correcte):

# inititialiser
rm(list=ls())
# charger librairies
require(rgl)

# Déclaration de fonctions
Integrand <- function(x,a,b) exp(-(a/b) * exp(b*x))
temp <- function(y,a,b) exp((a/b)*exp(b*y))*as.numeric(integrate(Integrand,lower=y,upper=Inf,a=a,b=b)["value"][1])

ey <- function(y,veca,vecb) {
  # Matrice censée contenir mes résultats:
  result <- matrix(data = 0, nrow = length(veca), ncol = length(vecb))
  # Balayage sur les vecteurs veca et vecb
  for (i in seq_along(veca))
   for (j in seq_along(vecb)) result[i,j] <- temp(y,veca[i],vecb[j])
  # Ajout : pour pouvoir obtenir la matrice à la fin des boucles
  return(result)
}

# Mes vecteurs contenant les différentes valeurs de ma fonction f():
A <- seq(from = 0.001, to = 0.01, by = 0.0001)
B <- seq(from = 0.1, to = 0.5, by = 0.01)

#test <- ey(0,A,B)
persp3d(ey(0,A,B),col=5)
persp3d(ey(10,A,B),col=6,add=TRUE)
persp3d(ey(20,A,B),col=7,add=TRUE)
#persp3d(ey(30,A,B),col=8,add=TRUE)
#persp3d(ey(40,A,B),col=8,add=TRUE)

Et en p.j, ce jolie graphique que je voulais obtenir (j'ai fait une capture d'écran, car la fenêtre RGL device n'a pas de menu pour pourvoir enregistrer le graph).
Mais 2 dernières ligne de commande sont en commentaire, parce que lorsque que je veux ajouter ey(30,A,B) et/ou ey(40,A,B), la fenêtre graphique bug car, pour certaines valeurs de a et b ma fonction ey n'est pas "définie" (l'intégral à l'intérieur de cette fonction est divergente pour ces valeurs je crois que c'est la raison). En jetant un oeil à ey(30,A,B) par exemple, certaines cellules de la matrice contiennent des valeurs numériques et d'autres contiennent soit "Inf" soit "NaN": c'est ce qui crée le bug.
Connaitrais-tu un moyen de "ploter" uniquement les points pour lesquelles ma fonction ey (par exemple ey(30,A,B)) est définie ?

[invitec5eb4b89]

Tu avais tout à fait raison de te méfier de la fonction représentée, je viens de m'apercevoir avec horreur que la fonction plt3d n'admet pas le même genre d'arguments que la fonction persp3d

Très beau graphe !

Pour la présence de valeurs singulières dans la matrice calculée, voici ce que je te propose :

# inititialiser
rm(list=ls())
# charger librairies
require(rgl)

# Déclaration de fonctions
Integrand <- function(x,a,b) exp(-(a/b) * exp(b*x))
temp <- function(y,a,b) exp((a/b)*exp(b*y))*as.numeric(integrate(Integrand,lower=y,upper=Inf,a=a,b=b)["value"][1])

ey <- function(y,veca,vecb) {
  # Matrice censée contenir mes résultats:
  result <- matrix(data = 0, nrow = length(veca), ncol = length(vecb))
  # Balayage sur les vecteurs veca et vecb
  for (i in seq_along(veca))
   for (j in seq_along(vecb)) result[i,j] <- temp(y,veca[i],vecb[j])
  # Ajout : pour pouvoir obtenir la matrice à la fin des boucles
  return(result)
}



# Mes vecteurs contenant les différentes valeurs de ma fonction f():
A <- seq(from = 0.001, to = 0.01, by = 0.0001)
B <- seq(from = 0.1, to = 0.5, by = 0.01)

persp3d(A,B,ey(0,A,B),col=5)
persp3d(A,B,ey(10,A,B),col=6,add=TRUE)
persp3d(A,B,ey(20,A,B),col=7,add=TRUE)


ey2 <- function(y,veca,vecb) {
  result <- matrix(NA,length(veca)*length(vecb),3)
  colnames(result) <- c("a","b","ey")
  k <- 1
  for (i in seq_along(veca))
   for (j in seq_along(vecb)) { 
    val <- temp(y,veca[i],vecb[j])
    if ( ! is.nan(val) & ! is.infinite(val)) result[k,] <- c(veca[i],vecb[j],val) ; k <- k+1 
    }
  return(result)
}

plot3d(ey2(30,A,B),col=8,add=TRUE)
plot3d(ey2(40,A,B),col=9,add=TRUE)

Ce n'est pas vraiment optimal, mais au moins tu peux avoir une idée de l'allure de la surface ! Je n'ai pas réussi à trouver un moyen de limiter le tracé de la surface aux zones dans lesquelles la fonction est définie, d'où la bidouille avec la fonction ey2 et le tracé complémentaire avec plot3d...

[inviteb8843a7e]

Merci beaucoup HigginsVincent ,

à vrai dire tu m'as montré le moyen de faire. Concrètement je ne savais pas comment remplacer des valeurs à l'intérieur d'une matrice en utilisant ces fameuses conditions . Il suffisait en fait de remplacer mes potentielles valeurs "-Inf", "Inf" et "NaN" qui me gênaient par la valeur "NA", ensuite la procédure graphique considère ainsi ces nouvelles valeurs comme "non-données" et ne va rien "ploter". Tu peux jeter un œil à mes 2 graphs: l'un propose de "ploter" les différences surfaces en 1 graph, l'autre est dynamique (il est permet de voir la courbe évoluer avec les changements de valeurs de y=>j'ai simplement fait ça en réalisant une boucle). Je pense que tout est fait maintenant

# inititialiser
rm(list=ls())
# charger librairies
require(rgl)

# Déclaration de fonctions
Integrand <- function(x,a,b) exp(-(a/b) * exp(b*x))
temp <- function(y,a,b) exp((a/b)*exp(b*y))*as.numeric(integrate(Integrand,lower=y,upper=Inf,a=a,b=b)["value"][1])

ey <- function(y,veca,vecb) {
  # Matrice censée contenir mes résultats:
  result <- matrix(data = 0, nrow = length(veca), ncol = length(vecb))
  # Balayage sur les vecteurs veca et vecb
  for (i in seq_along(veca))
   for (j in seq_along(vecb)) result[i,j] <- temp(y,veca[i],vecb[j])
   result[result==NaN | result==Inf | result==Inf] <- NA
  # Ajout : pour pouvoir obtenir la matrice à la fin des boucles
  return(result)
}


# Mes vecteurs contenant les différentes valeurs de ma fonction f():
A <- seq(from = 0.001, to = 0.01, by = 0.0001)
B <- seq(from = 0.1, to = 0.5, by = 0.01)

# Multi-plot
persp3d(ey(0,A,B),col=5,xlab="a",ylab="b", zlab="e(0)")
persp3d(ey(10,A,B),col=6,add=TRUE)
persp3d(ey(20,A,B),col=7,add=TRUE)
persp3d(ey(30,A,B),col=8,add=TRUE)
persp3d(ey(40,A,B),col=9,add=TRUE)
persp3d(ey(50,A,B),col=10,add=TRUE)


# Dynamic Multi-plot
for(i in 0:100) {
ey(i,A,B)
persp3d(ey(i,A,B),col=5,xlab="a",ylab="b", zlab=paste("e(",i,")"),zlim=c(0,50))
}

Reste maintenant à faire le même travail pour 2 autres fonctions que j'ai en réserve. ça ne devrait pas me poser plus de problème que ceux rencontrer jusqu'ici, il suffit en théorie de changer l'expression de mes fonctions du programme et le tour sera joué. Du moins j'espère...

Encore merci pour ta patience.

[invitec5eb4b89]

Aaah, d'accord, la fonction persp3d accepte les NA !
Ben merci à toi pour ton problème, j'ai appris des choses sur les fonctions graphiques 3D.

[inviteb8843a7e]

Aaah, d'accord, la fonction persp3d accepte les NA !
Ben merci à toi pour ton problème, j'ai appris des choses sur les fonctions graphiques 3D.

En fait, je ne pense que ça soit une question de savoir si la fonction persp3d accepte ou non les "NA", mais plus généralement, R comme d'autres logiciels statistiques, ne va pas représenter un point pour lequel il manque une coordonnée, ainsi par exemple:
le point de coordonnées (x,y,z)=(1.3,0.52,45) sera représenté mais le point (x,y,z)=(NA,0.14,36) ne le sera pas car R comprend qu'il n'a pas de valeur (mais il y a une donnée qui est "NA") sur la première coordonnée ( c'est comme si conceptuellement t'avais un blanc à la place disons (,0.14,36), mais attention un blanc pour R ne veut rien dire, j'ai bien dit conceptuellement ). Mais comme avant on avait des valeurs tel que NaN (Not a Number) ou Inf ou encore -Inf qui veulent dire respectivement que l'élément de la coordonnée en question n'est soit pas un nombre, soit l'infinie ou moins l'infini, ces dernières ne sont pas représentable graphiquement => d'où le bug graphique. Et ce bug ne dépend pas, si je ne me trompe pas (mais j'en suis quasi-sûr), de la procédure graphique. C'est plus lié au fonctionnement interne de R et de la signification des valeurs tel que NaN, Inf ou -Inf.
Pour finir, histoire de m'assurer d'avoir été compris , si par exemple j'avais dans ma matrice un point tel que (x,y,z)= (Inf ,0.14,36), le "Inf" ne veut pas dire que la coordonnées Inf n'existe pas ou et "non-donnée", ça n'est pas une erreur, ce point existe mais il est naturellement non représentable (l'infini n'est pas représentable) mais cela reste quant même le résultat de ma fonction.

Voilà, voilà pour l'interprétation du bug graphique. Et comme mon souci c'était d'empêcher ce bug en représentant les points pour laquelle la fonction était "graphiquement présentable", j'ai triché un peu en remplaçant mais vrai valeur Inf, -Inf ou NaN, par de fausse données "NA" qui veut dire pour R "pas de valeur, donc pas de point pour cette coordonnée".