Discussion :

[alan_ker]

Bonjour à tous,
Je commence tout juste à me rendre compte des possibilités qu'offre R pour l'analyse de données expérimentales en biologie mais aussi du casse tête que ça peut parfois représenter. Je lance donc un appel à l'aide sur un problème qui pourra (et j'espère) paraitre bête à un utilisateur chevronné.
Le but est d' établir un graphique de l'équation d'euler 1=somme(exp(-r*x)*lx*mx en fonction de r où
r : taux d'acroisement intrinseque d'une population (valeur inconnue mais comprise entre 0 et 1)
x : âge des bestioles
lx : taux de survie de la classe d'âge "x"
mx : fécondité de la classe d'âge "x"

Après avoir construit ma table de survie, j'obtiens :

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
> table_survie
    x   lx    mx
1   0 0.00  0.00
2   1 0.00  0.00
3   2 0.00  0.00
4   3 0.00  0.00
5   4 0.00  0.00
6   5 0.00  0.00
7   6 0.00  0.00
8   7 0.00  0.00
9   8 0.00  0.00
10  9 0.00  0.00
11 10 1.00 15.89
12 11 0.73  9.90
13 12 0.59  7.84
14 13 0.47  7.84
15 14 0.40  5.88
16 15 0.37  4.62
17 16 0.29  5.33
18 17 0.27  5.02
19 18 0.24  3.15
20 19 0.18  1.71
21 20 0.13  2.28
22 21 0.07  2.70
23 22 0.05  2.50
24 23 0.04  2.33
25 24 0.03  1.00

Si j'échantillonne un "r" au hasard dans une séquence entre 0 et 1 et que je répète l'opération 3 fois ça donne :

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
> sum(exp(-(sample(seq(0,1,0.0001),1)*x))*lx*mx)
[1] 3.600984
> sum(exp(-(sample(seq(0,1,0.0001),1)*x))*lx*mx)
[1] 0.2631770
> sum(exp(-(sample(seq(0,1,0.0001),1)*x))*lx*mx)
[1] 0.002341589

L'idée serait d'arriver à trouver une fonction qui permette de répéter automatiquement cette opération (par exemple 1000 fois) avec à chaque fois un "r" pris au hasard et de tracer ensuite un graphique avec en ordonée y ["sum(exp(-(sample(seq(0,1,0.0001),1)*x))*lx*mx)"] et en abscisse x [ce "r" pris au hasard dans cette séquence].

Au moment ou y se rapproche de 1, je pourrais alors savoir quelle valeur vaut "r" et avoir accès au taux d'accroisement intrinséque de la population.

J'espere avoir été clair, si ce n'est pas le cas je reste à la disposition du premier internaute qui aura une question,
Alan

[vchouraki]

Bonjour,

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

?replicate

Avec votre exemple ("répéter 1000 fois") :

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

replicate(1000,sum(exp(-(sample(seq(0,1,0.0001),1)*x))*lx*mx)

HTH

Vincent Chouraki

[alan_ker]

Bonjour vincent,
Merci pour la réponse, je commençais à deseperer. Effectivement la fonction replicate enlève une sacrée épine du pied. Mais il reste encore à placer ces valeurs sur un graph en fonction du "r" échantillonné dans la séquence entre 0 et 1.
Voili voilou ce que j'ai essayé :

> r<-seq(0,0.999,0.001)
> y<-replicate(1000,sum(exp(-(sample(r,1)*x))*lx*mx))
> plot(y~r,ylim=c(0,2))

Le problème c'est que je n'obtiens pas un graphique avec une courbe qui se rapproche de 1 pour une valeur de "r" se rapprochant du véritable taux d'accroissement.
Je soupçonne ma formule d'affecter un "y" pris au hasard parmi les 1000 répliques de calculs à un "r" lui aussi pris au hasard dans la séquence entre 0 et 1.
Le but du jeu est d'afficher pour un "r" tiré le "y" qui lui correspond.
Je suis encore au labo mais j'essaie de m'y coller en rentrant. Je suis preneur de toute suggestion ...

Quoiqu'il en soit, merci bien
A.

[vchouraki]

Rebonjour,

replicate() est une forme simplifiée d'une des fonctions de la famille "apply" et qui, grosso modo, remplace les boucles "for" dans R. Voir l'aide pour plus de détail.

Concernant le comportement de "plot" dans votre cas, ça ne m'étonnerait pas qu'il fasse le graphique de la dernière réplication.

Une solution peut consister en l'écriture d'une fonction qui fera tous les calculs et le graphique (avec sauvegarde sous forme de fichier pdf par exemple) et ensuite de la passer à la moulinette "lapply" (plus complet que replicate pour générer des noms de fichier par exemple).

Par exemple :

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
test <- function(iter){
 
  a <- sample(1:10000,100)
 
  c <- sample(1:10,1)
 
  b <- a^c
 
  pdf ( file = paste ( "test_" , iter , ".pdf" , sep = ""))
 
  plot(b~a)
 
  dev.off ()
 
}
 
lapply(1:10,test)

La fonction test() crée des données, puis sauvegarde le graphique de ces données au format pdf dans le répertoire courant avec un nom de fichier comprenant le numéro de l'itération.

La fonction lapply() prend chaque élément du premier argument (élément unique s'il s'agit d'un vecteur, vecteur-colonne s'il s'agit d'un data.frame, élément d'une liste quel qu'il soit s'il s'agit d'une liste) et lui applique la fonction test()

En adaptant cet exemple à vos besoins, vous devriez vous en sortir.

HTH

Vincent Chouraki