diff options
Diffstat (limited to 'Master/texmf-dist/doc/latex/pbsheet/xpl/pgm/probadis.m')
-rw-r--r-- | Master/texmf-dist/doc/latex/pbsheet/xpl/pgm/probadis.m | 65 |
1 files changed, 65 insertions, 0 deletions
diff --git a/Master/texmf-dist/doc/latex/pbsheet/xpl/pgm/probadis.m b/Master/texmf-dist/doc/latex/pbsheet/xpl/pgm/probadis.m new file mode 100644 index 00000000000..7f472ad4045 --- /dev/null +++ b/Master/texmf-dist/doc/latex/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; |