summaryrefslogtreecommitdiff
path: root/macros/latex/contrib/pbsheet/xpl/pgm
diff options
context:
space:
mode:
Diffstat (limited to 'macros/latex/contrib/pbsheet/xpl/pgm')
-rw-r--r--macros/latex/contrib/pbsheet/xpl/pgm/probadis.m65
-rw-r--r--macros/latex/contrib/pbsheet/xpl/pgm/rdiscr.m72
-rw-r--r--macros/latex/contrib/pbsheet/xpl/pgm/rint.m25
-rw-r--r--macros/latex/contrib/pbsheet/xpl/pgm/simbin.m55
4 files changed, 217 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;
diff --git a/macros/latex/contrib/pbsheet/xpl/pgm/rdiscr.m b/macros/latex/contrib/pbsheet/xpl/pgm/rdiscr.m
new file mode 100644
index 0000000000..b3346178e4
--- /dev/null
+++ b/macros/latex/contrib/pbsheet/xpl/pgm/rdiscr.m
@@ -0,0 +1,72 @@
+function X = rdiscr(num,x,p)
+%X = RDISCR(num,x,p)
+% Input num positive integer or a vector [lig,col] of integers
+% 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 X num-vector or a num-matrix of random numbers
+% 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 num réalisations ou une matrice [lig,col] de réalisations
+% pseudo-i.i.d 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 3
+% - 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
+% - le paramètre num est un entier positif non nul une un couple de ce type
+%
+% See also PROBADIS.
+%
+% ### 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
+%
+
+% Voir les commentaires dans le code de la fonction probadis.
+% Faire 'type probadis' pour cela.
+% Le code qui suit pourrait beaucoup gagner en rapidité sur une machine //
+% Il est possible de l'améliorer en stockant un arbre construit avec les
+% tests utilisés par les valeurs déjà générées. Il est surtout aussi possible
+% d'imiter le code matriciel de la fonction de répartition inverse binomiale
+% de la fonction qbinom de Stixbox, appelée par rbinom. Exercice !
+%
+% D'autres algorithmes sont possibles. Par exemple, on pourrait procéder
+% avec des réalisations de Bernoulli i.i.d. (obtenues facilement avec rand)
+% pour choisir l'intervalle :
+% La probabilité d'être <= x(1) est p(1)
+% Sinon, la probabilité d'être <= x(2) est p(2)/sum(p(2:n))
+% etc. Nul besoin de commencer par x(1), et l'on peut adapter l'arbre utilisé
+% aux poids p(i) de façon à faire un nombre de tests optimal. Exercice !
+
+if length(num) == 1
+ num = [num 1];
+else
+ num = reshape(num,1,2);
+end
+
+n = length(x); % le nombre d'atomes.
+U = rand(num); % réalisations de loi uniforme sur [0,1].
+X = repmat(x(n),num(1),num(2)); % par défaut, la valeur est la plus grande
+for l = 1:num(1) % lignes
+for c = 1:num(2) % colonnes
+ a = 0; b = p(1); % [a,b] = sous-inter. de proba p(i) pour l'uniforme.
+ for i = 1:n-1 % parcours de tous les sous-intervalles juxtaposés.
+ if ((U(l,c) >= a) & (U(l,c) < b))
+ X(l,c) = x(i);
+ break; % on a trouvé le bon intervalle, on sort.
+ end
+ a = b; b = b + p(i+1); % on passe à l'intervalle suivant.
+ end
+end
+end
+return;
diff --git a/macros/latex/contrib/pbsheet/xpl/pgm/rint.m b/macros/latex/contrib/pbsheet/xpl/pgm/rint.m
new file mode 100644
index 0000000000..9e8dd3c850
--- /dev/null
+++ b/macros/latex/contrib/pbsheet/xpl/pgm/rint.m
@@ -0,0 +1,25 @@
+function X = rint(num,k);
+% X = RINT(num,k)
+% Input num positive integer or a vector [lig,col] of integers
+% k positive integer
+% Output X num-vector or a num-matrix of random numbers
+% chosen from uniform distribution on {1,...,k}
+%
+% Renvoie num réalisations ou une matrice [lig,col] de réalisations
+% pseudo-indépendantes de la loi uniforme discrète sur les k premiers
+% entiers non nuls {1,...,k}.
+%
+% La méthode consiste à considérer la partie entière de 1+kU où U suit une loi
+% uniforme sur [0,1]. Cette dernière s'obtient via la fonction rand, et les
+% valeurs renvoyées par des appels successifs sont donc pseudo-indépendantes.
+%
+% ### 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
+%
+
+if length(num) == 1
+ num = [num 1];
+end
+X = ceil(k * rand(num));
+return;
diff --git a/macros/latex/contrib/pbsheet/xpl/pgm/simbin.m b/macros/latex/contrib/pbsheet/xpl/pgm/simbin.m
new file mode 100644
index 0000000000..f9c618a8a6
--- /dev/null
+++ b/macros/latex/contrib/pbsheet/xpl/pgm/simbin.m
@@ -0,0 +1,55 @@
+%
+% ### 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
+%
+% Ce bout de code permet de simuler une loi binomiale B(n,p) et de comparer
+% graphiquement les moyennes empiriques avec la moyenne théorique.
+% Tous les vecteurs sont des vecteurs ligne.
+%
+% Nota bene : la bibliothèque Stixbox fournit la fonction rbinom qui permet
+% de s'affranchir des calculs détaillés ici, cf. fin du présent fichier.
+%
+clear r n p q P Q C B X
+clf
+%
+r = input('Nombre maximal de réalisations ? = ');
+n = input('Taille n de la loi binomiale B(n,p), qui a n+1 atomes ? = ');
+p = input('Valeur du paramètre p de la loi binomiale B(n,p) ? = ');
+% p = proba de gagner à pile ou face = proba de 1 dans la Bernoulli sur {0,1}
+% q = proba de perdre à pile ou face = proba de 0 dans la Bernoulli sur {0,1}
+% B(n,p) est la loi de la somme de n v.a. i.i.d. de Bernoulli de ce type,
+% i.e. la nième puissance de convolution de cette loi de Bernoulli. Elle
+% représente la loi du nombre de succés à pile ou face en n lancés.
+%
+disp(sprintf('Calcul de la loi binomiale B(%d,%f)',n,p))
+q = 1 - p;
+P = [1 , cumprod(p * ones(1,n))]; % puissances croissantes de p
+Q = [fliplr(cumprod(q * ones(1,n))), 1]; % puissances décroissantes de q
+C = [1 , cumprod([n:-1:1]) ./ cumprod([1:n])]; % coef. du binôme de i = 1 à n
+B = C .* P .* Q; % vecteur des poids de B(n,p)
+%
+disp(sprintf('Génération aléatoire de %d réalisations de B(%d,%f)',r,n,p))
+X = rdiscr([1,r],[0:n],B); % échantillonnage
+% alternative :
+%X=[]
+%for i=1:r
+% X = [X probadis([0:n],B)];
+%end
+%
+disp(sprintf('Tracé des graphique.'))
+plot(cumsum(X) ./ [1:length(X)],'b') % tracé des moyennes empiriques
+title('Loi des Grands Nombres')
+xlabel('Nombre de realisations')
+ylabel('Moyennes empiriques')
+hold on
+plot(n * p * ones(1,r),'r--') % tracé de la moyenne théorique
+legend('Empirique','Theorique')
+hold off
+%
+% Avec Stixbox, inutile de calculer B et d'appeler rdiscr([1,r],[0:n],B)
+% puisque qu'un simple rbinom([1,r],n,p) suffit.
+% Même si l'on décide d'utiliser quand même prodadis, les coefficients
+% binomiaux nécessaires au calcul de B peuvent se calculer beaucoup plus vite
+% en utilisant la fonction Stixbox bincoef qui fait appel à la fonction gamma.
+%