diff options
author | Karl Berry <karl@freefriends.org> | 2006-01-09 00:49:07 +0000 |
---|---|---|
committer | Karl Berry <karl@freefriends.org> | 2006-01-09 00:49:07 +0000 |
commit | 007f67a693e4d031fd3d792df8e4d5f43e2cb2e7 (patch) | |
tree | 90d17e00e572ecb1e24764b6f29c80e098b08d29 /Master/texmf-dist/doc/latex/pbsheet/xpl/pgm | |
parent | 950209b26f70aa87ed07c54f82a95b6f03b7c3a0 (diff) |
doc/latex
git-svn-id: svn://tug.org/texlive/trunk@84 c570f23f-e606-0410-a88d-b1316a301751
Diffstat (limited to 'Master/texmf-dist/doc/latex/pbsheet/xpl/pgm')
4 files changed, 217 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; diff --git a/Master/texmf-dist/doc/latex/pbsheet/xpl/pgm/rdiscr.m b/Master/texmf-dist/doc/latex/pbsheet/xpl/pgm/rdiscr.m new file mode 100644 index 00000000000..b3346178e47 --- /dev/null +++ b/Master/texmf-dist/doc/latex/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/Master/texmf-dist/doc/latex/pbsheet/xpl/pgm/rint.m b/Master/texmf-dist/doc/latex/pbsheet/xpl/pgm/rint.m new file mode 100644 index 00000000000..9e8dd3c8502 --- /dev/null +++ b/Master/texmf-dist/doc/latex/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/Master/texmf-dist/doc/latex/pbsheet/xpl/pgm/simbin.m b/Master/texmf-dist/doc/latex/pbsheet/xpl/pgm/simbin.m new file mode 100644 index 00000000000..f9c618a8a69 --- /dev/null +++ b/Master/texmf-dist/doc/latex/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. +% |