summaryrefslogtreecommitdiff
path: root/macros/latex/contrib/pbsheet/xpl/pgm/probadis.m
diff options
context:
space:
mode:
Diffstat (limited to 'macros/latex/contrib/pbsheet/xpl/pgm/probadis.m')
-rw-r--r--macros/latex/contrib/pbsheet/xpl/pgm/probadis.m65
1 files changed, 65 insertions, 0 deletions
diff --git a/macros/latex/contrib/pbsheet/xpl/pgm/probadis.m b/macros/latex/contrib/pbsheet/xpl/pgm/probadis.m
new file mode 100644
index 0000000000..7f472ad404
--- /dev/null
+++ b/macros/latex/contrib/pbsheet/xpl/pgm/probadis.m
@@ -0,0 +1,65 @@
+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;