Obtenir un vecteur de valeurs inverses de distributions cumulées de probabilité

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

Boury frederic
Messages : 12
Enregistré le : 15 Jan 2018, 08:09

Obtenir un vecteur de valeurs inverses de distributions cumulées de probabilité

Messagepar Boury frederic » 03 Fév 2018, 12:17

Bonjour :
Mon problème est le suivant...
Dans mon programme, j'ai n variables aléatoires qui ont toutes une distribution de probabilité distincte (représentée par une distribution cumulée sous la forme d'un vecteur de p positions F(1) =0 et F(p) =1).
J'effectue ensuite un tirage aléatoire uniforme pour chaque variable (soit X <- runif(n)) et je souhaite obtenir en retour le vecteur Y(n) qui donne les valeurs telles que X = F(Y)
Je sais programmer ceci en mode indiciel (cf l'exemple ci dessous) mais je pense qu'il est possible de le faire de façon beaucoup plus concise et performante en mode vectoriel...
J'ai bien évidemment pensé à utiliser la fonction apply (du genre apply(X, fonction) ) mais je ne vois pas dans la documentation comment passer la matrice F en plus de X comme paramètres (faire cbind (X, F) ne marche pas).
Merci e vos conseils...

========== EXEMPLE avec n = 20 ===========
# gènère 20 distributions cumulées distinctes F
F <- matrix(data =1, nrow=20, ncol = 10)
for (i in 1:20) {
for (j in 1:9) { F[i,10-j] = F[i, 10 - (j-1)]-runif(1)/5}}
F <- cbind(rep(0,20), matrix(data =pmax(0,F), nrow=20, ncol = 10))

# la fonction ord renvoie le N° d'ordre d'une valeur comprise entre 0 et 1 dans la distribution D
ord <- function (x, D) { max(which(x > D)) }

# obtention de Y = F-1(X) en notation indicée
X <- runif(20) ; Y <- rep(0, 20)
for (i in 1:20) {Y[i]= ord(X[i],F[i,])}

Pierre-Yves Berrard
Messages : 1029
Enregistré le : 12 Jan 2016, 23:30

Re: Obtenir un vecteur de valeurs inverses de distributions cumulées de probabilité

Messagepar Pierre-Yves Berrard » 03 Fév 2018, 13:17

À regarder simplement le résultat (et pas l'intérieur de la boucle), j'ai l'impression que la fonction replicate pourrait vous être utile :

Code : Tout sélectionner

mat <- replicate(20, c(0, sort(runif(9)), 1))
t(mat)
PY

Boury frederic
Messages : 12
Enregistré le : 15 Jan 2018, 08:09

Re: Obtenir un vecteur de valeurs inverses de distributions cumulées de probabilité

Messagepar Boury frederic » 03 Fév 2018, 14:03

Bonjour Monsieur Berrard :

La solution que vous indiquez est en effet beaucoup plus efficace que la mienne pour générer les 20 distributions cumulées F...

Mais ma question était différente.

La matrice F(n, p) des n distributions cumulées quelconques (réparties sur p colonnes) est réputée être donnée en entrée. De même, le vecteur aléatoire X(n). Mon souci est d'obtenir Y(n) = F-1(X)

Comme je l'indiquais, je sais y arriver en mode indiciel mais je pense que ça n'est pas optimal. Ci après la manière que j'ai trouvée à ce sujet :

ord <- function (x, D) { max(which(x > D)) }
X <- runif(20) ; Y <- rep(0, 20)
for (i in 1:20) {Y[i]= ord(X[i],F[i,])}

Je cherche comment obtenir directement le résultat avec des instructions vectorielles... peut être avec la fonction apply, mais je n'arrive pas à la faire fonctionner sur ce cas. Je pensais à quelque chose du genre Y <- apply(cbind(X,F),1, ord) mais ça ne marche pas !

Pierre-Yves Berrard
Messages : 1029
Enregistré le : 12 Jan 2016, 23:30

Re: Obtenir un vecteur de valeurs inverses de distributions cumulées de probabilité

Messagepar Pierre-Yves Berrard » 03 Fév 2018, 14:28

Pardon j'ai mal lu.
Donc on peut utiliser mapply (une variante de apply) qui permet de parcourir en parallèle plusieurs arguments. Ici, X et les colonnes de F transposé en data.frame.

Code : Tout sélectionner

mapply(
  ord,
  X,
  as.data.frame(t(F))
)
PY

Boury frederic
Messages : 12
Enregistré le : 15 Jan 2018, 08:09

Re: Obtenir un vecteur de valeurs inverses de distributions cumulées de probabilité

Messagepar Boury frederic » 03 Fév 2018, 17:14

Bonjour M. Berrard :

Merci encore pour votre suggestion : cela marche parfaitement. Néanmoins, je cherchais à vectoriser l'expression pour avoir de meilleures performances mais ça n'est pas le cas... Sur 300.000 runs et un vecteur de 53 valeurs, le temps de traitement est de 260 secondes avec l'expression 'indicielle' contre 480 secondes avec mapply.

C'est bizarre, car normalement les expressions vectorisées sont plus performantes. Avez vous une explication ?

Pour information, la matrice v est de 60 lignes et 23 colonnes

expression indicielle :
uni <- runif(nbcp)
for (i in 1:nbcp) {u2[i]= rnk_UV(uni[i],v[i,])}

expression vectorielle :
uni <- runif(nbcp)
u2 <- mapply(rnk_UV, uni, as.data.frame(t(v)))

# fonction rnk_UV donne le rang de l'aléa U sur un vecteur de probabilités cumulées V
rnk_UV <- function (u, V) {max( which( V < u))}

Pierre-Yves Berrard
Messages : 1029
Enregistré le : 12 Jan 2016, 23:30

Re: Obtenir un vecteur de valeurs inverses de distributions cumulées de probabilité

Messagepar Pierre-Yves Berrard » 03 Fév 2018, 17:28

Bonjour,

Les fonctions *apply sont souvent plus élégantes mais pas systématiquement plus rapides.

Dans ma solution, j'ai recours à as.data.frame() et t() qui prennent du temps. Il y a sûrement moyen d'optimiser cette partie.
Est-ce plus rapide avec

Code : Tout sélectionner

purrr::map2(X, as.data.frame(t(F)), ord)

Ensuite, votre boucle est assez efficace car vous avez pris la peine d'initialiser un vecteur au préalable.
PY

Boury frederic
Messages : 12
Enregistré le : 15 Jan 2018, 08:09

Re: Obtenir un vecteur de valeurs inverses de distributions cumulées de probabilité

Messagepar Boury frederic » 03 Fév 2018, 19:37

Bien reçu...
J'ai essayé de lancer votre instruction, mais j'ai un message d'erreur... Cf ci dessous.
je pense que je n'ai pas compris ce que vous voulez faire...
Pouvez vous m'clairer davantage ?
;-)

Y <- map2(X, as.data.frame(t(F)), ord)
Error in map2(X, as.data.frame(t(F)), ord) : impossible de trouver la fonction "map2"

> Y <- purrr::map2(X, as.data.frame(t(F)), ord)
Error in loadNamespace(name) : aucun package nommé ‘purrr’ n'est trouvé

> ? purrr
No documentation for ‘purrr’ in specified packages and libraries:
you could try ‘??purrr’

> ? map2
No documentation for ‘map2’ in specified packages and libraries:
you could try ‘??map2’

Pierre-Yves Berrard
Messages : 1029
Enregistré le : 12 Jan 2016, 23:30

Re: Obtenir un vecteur de valeurs inverses de distributions cumulées de probabilité

Messagepar Pierre-Yves Berrard » 09 Fév 2018, 10:44

Bonjour,

Les messages d'erreurs signalent qu'il faut installer le package purrr.
PY

Logez Maxime
Messages : 3138
Enregistré le : 26 Sep 2006, 11:35

Re: Obtenir un vecteur de valeurs inverses de distributions cumulées de probabilité

Messagepar Logez Maxime » 09 Fév 2018, 13:43

Bonjour,

Si j'ai bien compris la question :

Code : Tout sélectionner

rowSums(F < X)

library(microbenchmark)
microbenchmark(map2(X, as.data.frame(t(F)), ord),
  rowSums(F < X),
  mapply(ord, X, as.data.frame(t(F))), times = 1e3)
Unit: microseconds
                                expr     min      lq       mean  median      uq      max neval
   map2(X, as.data.frame(t(F)), ord) 157.448 164.244 178.994001 166.887 173.684 2441.378  1000
                      rowSums(F < X)   3.776   4.532   5.979415   6.042   6.796   30.962  1000
 mapply(ord, X, as.data.frame(t(F))) 149.897 155.937 173.075962 158.580 165.755 1067.018  1000
Cordialement,
Maxime

Pierre-Yves Berrard
Messages : 1029
Enregistré le : 12 Jan 2016, 23:30

Re: Obtenir un vecteur de valeurs inverses de distributions cumulées de probabilité

Messagepar Pierre-Yves Berrard » 09 Fév 2018, 16:16

Logez Maxime a écrit :

Code : Tout sélectionner

rowSums(F < X)

Imbattable !
Quand X est répliqué pour s'ajuster à la longueur de F, c'est fait par colonne puis par ligne, c'est bien pour cela que ça marche ?
PY

Logez Maxime
Messages : 3138
Enregistré le : 26 Sep 2006, 11:35

Re: Obtenir un vecteur de valeurs inverses de distributions cumulées de probabilité

Messagepar Logez Maxime » 10 Fév 2018, 07:48

Bonjour,

ça marche parce qu'il y a autant de valeurs de X que de lignes de F. Dans une matrice les données sont stockées en colonnes. Si on a un vecteur de la même longueur que le nombre de lignes comme ici, alors il compare la première valeur du vecteur avec la première valeur de la première colonne, puis la deuxième valeur du vecteur avec la deuxième valeur de la première colonne, etc. Une fois le vecteur finit il le recycle et ça va tomber pile poile avec la première valeur de la deuxième colonne.

Cordialement,
Maxime

Boury frederic
Messages : 12
Enregistré le : 15 Jan 2018, 08:09

Re: Obtenir un vecteur de valeurs inverses de distributions cumulées de probabilité

Messagepar Boury frederic » 11 Juin 2018, 08:06

Bonjour Monsieur Logez :

Je suis désolé de revenir aussi tardivement vers vous, mais j'avais manqué votre contribution à l'époque (et ce n'est que de façon récente, à travers une nouvelle question sur le forum) que j'ai enfin vu votre réponse.

Tout ça pour dire que j'ai regardé attentivement vos explications ci dessus et que je n'ai retrouvé à ce sujet que les caractéristiques de la fonction que j'avais plus ou moins bien implémentée soit : function rank_uV (u,V) { pmax(which(V<u))}

Mon souci est donc le suivant...

Rien dans la documentation de la fonction rowSums n'indique qu'en introduisant un test dedans soit a <- rowSums( F<x) on va obtenir le quantile de x sur la distribution F.

D'ailleurs, en toute rigueur, on s'attendrait plutôt à obtenir la distribution cumulée de F (la somme de F) jusqu'en x, pas la valeur du quantile inférieur de x sur F...

D'où ma question : où peut figurer la documentation avancée de rowSums qui explique cette fonctionnalité ?

Désolé d'insister, mais le fait de savoir si il existe un ou des sites avec des niveaux de documentation plus avancées sur les fonctions de R est un sujet sensible... En effet, on achoppe souvent sur le caractère sybilin et peu compréhensibles des documentations en ligne de R.


Bioen à vous
Frédéric Boury


Retourner vers « Questions en cours »

Qui est en ligne

Utilisateurs parcourant ce forum : Aucun utilisateur enregistré et 1 invité