function res = probadis(x,p) %res = PROBADIS(x,p) % Input x vector of real numbers (support points in IR) % p vector of probability weights associated to x % i.e. non-negative real numbers such that sum(p) == 1 % Output res random number chosen from the finite discrete % distribution on x(1),...,x(n) with probability weights % p(1),...,p(n) where n == length(x) == length(p) % % Renvoie UNE réalisation de la loi discrète à support fini sur IR dont les % points de support sont les composantes du vecteur x, et les poids sont les % composantes du vecteur p. Donc x(i) a une probabilité p(i) d'être renvoyé. % % La méthode consiste à écrire l'intervalle [0,1] comme une réunion disjointe % d'intervalles de longueurs p(1),..., p(n) puis à regarder au quel appartient % la réalisation d'une loi uniforme obtenue par un appel à rand. Les valeurs % renvoyées par des appels successifs sont donc pseudo-indépendantes. % Pour des raisons d'efficacité, les conditions nécessaires suivantes ne sont % pas contrôlées par cette fonction : % - le nombre de paramètres passés est exactement 2 % - les deux paramètres x et p sont bien des vecteurs et sont de même longueur % - les composantes de p sont positives ou nulles et leur somme vaut 1. % % See also RDISCR. % % ### Copyright (C) D. Chafaï, 2003-12-06. % ### http://www.lsp.ups-tlse.fr/Chafai/agregation.html % ### Licence GNU General Public License http://www.gnu.org/copyleft/gpl.html % % On pourrait implémenter cette fonction de la façon suivante : % % INDICES = find(cumsum(p) >= rand); % res = x(INDICES(1)); % return; % % Cette méthode est correcte mais inefficace car elle ne tient pas compte de % la monotonie de cumsum(p), ce qui entraine des tests inutiles une fois que % la valeur critique a été franchie. En jargon, c'est un 'firstmatch' qu'il % nous faut, pas un 'matchall'. Le nombre de if impliqués dans find est % toujours égal à la taille de ce qu'on lui passe en paramètre. % % Morale de l'histoire : la brièveté d'un code n'assure pas sa performance ! % Et l'abscence de if dans un code ne signifie pas qu'il ne fait pas appel % indirectement à des if, et encore moins que cela est fait de manière % optimale ! La fonction find n'est qu'une boucle for contenant un if. % % Ci-dessous, nous utilisons une version plus rapide dans le plus pur style % for-if avec un nombre de if optimal. On pourrait adapter l'ordre des % intervalles testés (et donc l'arbre associéé) aux poids p(i) de façon à % tester d'abord les intervalles les plus probables. Est-ce vraiment mieux ? % Exercice ! n = length(x); % le nombre d'atomes. r = rand; % une réalisation de loi uniforme sur [0,1]. a = 0; b = p(1); % [a,b] = sous-intervalle de proba p(i) pour l'uniforme. for i = 1:n-1 % parcours de tous les sous-intervalles juxtaposés. if ((r >= a) & (r < b)) res = x(i); return; % on a trouvé le bon intervalle, on sort. end a = b; b = b + p(i+1); % on passe à l'intervalle suivant. end res = x(n); % le bon intervalle est le dernier. return;