Je rassemble dans cette entrée quelques faits algorithmiques et informatiques qui sont généralement « bien connus » (et franchement assez basiques) mais souvent utiles, et qu'il est possiblement difficile de trouver rassemblés en un seul endroit. Le problème général est de tirer algorithmiquement des variables aléatoires selon différentes distributions, typiquement à partir d'un générateur aléatoire qui produit soit des bits aléatoires (indépendants et non biaisés) soit des variables aléatoires réelles (indépendantes) uniformément réparties sur [0;1]. Je parle d'algorithmique, mais ce n'est pas uniquement sur un ordinateur : ça peut être utile même dans la vie réelle, par exemple si on a une pièce avec laquelle on peut tirer à pile ou face et qu'on veut s'en servir pour jouer à un jeu qui réclame des dés à 6 faces, ou si on a des dés à 6 faces et qu'on veut jouer à un jeu d'aventure qui réclame des dés à 20 faces.
Comment tirer des nombres aléatoires en conditions
adversariales ? Je commence par ce problème-ci qui n'a pas de
rapport direct avec la suite, mais que je trouve quand même opportun
de regrouper avec : Alice et Bob veulent jouer à pile ou face, ou plus
généralement tirer un dé à n faces, mais ils n'ont pas de
pièce ou de dé en lequel ils fassent tous les deux confiance. Par
exemple, Alice a sa pièce fétiche que Bob soupçonne d'être truquée et
symétriquement (ou peut-être même que chacun est persuadé de pouvoir
tirer des nombres aléatoires dans sa tête mais ne fait évidemment pas
confiance à l'autre). La solution est la suivante : chacun fait un
tirage avec son propre moyen de son côté, sans connaître le
résultat de l'autre, et on combine ensuite les résultats selon
n'importe quelle opération (choisie à l'avance !) qui donne tous
les n résultats possibles pour chaque valeur fixée d'une
quelconque des entrées
(un carré
latin, par exemple une loi de groupe) ; par exemple, s'il s'agit
de tirer à pile ou face, on peut décider (à l'avance !) que le
résultat sera pile
(0) si les deux pièces ont donné le même
résultat et face
(1) si elles ont donné un résultat différent ;
s'il s'agit de dés à n faces donnant un résultat entre
0 et n−1, on fait la somme modulo n
(c'est-à-dire qu'on fait la somme et qu'on soustrait n si
elle vaut au moins n, pour se ramener à un résultat entre
0 et n−1). Bien sûr, il faut un protocole pratique pour
faire en sorte que chacun fasse son tirage sans connaître le résultat
de l'autre (sinon, s'il a moyen de tricher, il pourra adapter le
résultat en conséquence) : physiquement, chacun peut faire son tirage
en secret et écrire le résultat secrètement sur un papier placé dans
une enveloppe scellée, qu'on ouvrira une fois les deux tirages
effectués (en fait, il n'y a que le premier tirage qui a besoin d'être
fait de la sorte) ; cryptographiquement, on procède à
une mise en
gage (typiquement au moyen d'une fonction de hachage, mais je ne
veux pas entrer dans ces questions-là). On peut bien sûr généraliser
à plus que deux joueurs (en faisant la somme modulo n de
nombres tirés par chacun des participants). Le protocole garantit que
le résultat sera un tirage uniforme honnête si l'un au moins
des participants désire qu'il le soit (et a les moyens de réaliser un
tirage honnête) : bien sûr, si aucun des participants ne le souhaite,
c'est leur problème, donc on s'en fout. Même si on abandonne toute
prétention à ce que les participants tirent leur valeur
aléatoirement et qu'on s'imagine qu'ils la choisissent, tant
qu'il s'agit d'un jeu à somme nulle, la stratégie optimale est bien de
tirer au hasard (et encore une fois, s'ils veulent coopérer pour un
autre résultat, tant qu'il n'y a pas d'autre partie impliquée, c'est
leur problème, de même s'ils s'imaginent pouvoir faire mieux que le
hasard en utilisant, par exemple, une prédiction psychologique).
Je ne sais plus où j'avais lu que ce protocole a été découvert (il l'a certainement été de nombreuses fois !) à la renaissance. Dans mon souvenir, le découvreur proposait même que, pour une question de la plus haute importance, on demande au pape de faire un des tirages en plus de tous les autres participants. Je ne sais d'ailleurs pas si ce protocole a un nom standard.
Comment tirer une variable de Bernoulli de
paramètre p à partir de bits aléatoires ?
Autrement dit, ici, on a fixé p, et on veut faire un tirage
aléatoire qui renvoie oui
avec probabilité p
et non
avec probabilité 1−p, et pour ça, on dispose
simplement d'une pièce qui renvoie des bits aléatoires en tirant à
pile (0) ou face (1), et on souhaite effectuer les tirages de façon
économique. Par exemple, combien de tirages de pièce faut-il, en
moyenne, pour générer un événement de probabilité 1/3 ? Il s'avère,
en fait, que quel que soit p on peut s'en tirer avec
deux (2) tirages en moyenne (je veux dire en espérance
). Pour
cela, on peut procéder ainsi : on effectue des tirages répétés et on
interprète les bits aléatoires ainsi produits comme l'écriture binaire
d'un nombre réel x uniformément réparti entre 0 et 1 : on
compare x à 1−p en binaire, c'est-à-dire qu'on
s'arrête dès qu'on dispose d'assez de bits pour pouvoir décider
si x < 1−p
ou x > 1−p (on peut considérer le
cas x = 1−p comme s'il était impossible vu qu'il
est de probabilité 0), en notant qu'on va
avoir x < 1−p lorsque le k-ième
bit tiré est 0 et que le k-ième bit de l'écriture binaire
de p vaut 0, et x > 1−p lorsque
le k-ième bit tiré est 1 et que le k-ième bit
de p vaut 1 ; et si x < 1−p on
renvoie non
, sinon oui
. Concrètement, donc, faire des
tirages aléatoires jusqu'à ce que le k-ième bit tiré soit
égal au k-ième bit de l'écriture de p, et alors
s'arrêter et renvoyer ce bit-là. Il est clair que cet algorithme
fonctionne, mais pour qu'il soit encore plus évident qu'il conduit à
faire deux tirages en moyenne, on peut le reformuler de la façon
encore plus élégante suivante (il suffit d'échanger les résultats
0 et 1 pour x, qui sont complètement symétriques, lorsque
le bit correspondant de p vaut 0) : tirer des bits
aléatoires jusqu'à tomber sur 1, et lorsque c'est le cas, s'arrêter et
renvoyer le k-ième bit de p
(où k est le nombre de bits aléatoires qui ont été tirés).
Je trouve ça incroyablement élégant et astucieux (même si c'est très
facile), et je ne sais pas d'où sort ce truc. (Cela revient encore à
tirer une variable aléatoire k distribuée selon une loi
géométrique d'espérance 1, comme je l'explique plus bas, c'est-à-dire
valant k avec probabilité (½)k+1, et
renvoyer le (k+1)-ième
bit bk+1 de p, ce qui,
quand on écrit p =
∑k=0+∞ bk+1·(½)k+1,
est finalement assez évident.)
Comment tirer un entier aléatoire entre 0 et n−1
à partir de bits aléatoires ? (Ma première réaction en
entendant ce problème a été de dire : considérer x uniforme
dont l'écriture binaire est donnée par la suite des bits tirés,
générer suffisamment de bits pour calculer la valeur de
⌊n·x⌋, où ⌊—⌋ désigne la partie entière, et
renvoyer celle-ci. Ceci fonctionne, mais ce n'est pas le plus
efficace. Un autre algorithme avec rejet consiste à
générer r := ⌈log(n)/log(2)⌉ bits, qui, lus en
binaire, donnent un entier aléatoire c entre 0 et
2r−1, renvoyer c s'il
est <n, et sinon tout recommencer. Mais ce n'est pas
très efficace non plus, quoique dans des cas un peu différents.) Je
décris ce problème plus en détails
dans ce
fil Twitter, mais donnons juste l'algorithme : on utilise deux
variables internes à l'algorithme, notées v
et c, qu'on initialise par v←1 et c←0
(il s'agit d'un réservoir d'entropie
, et la garantie est
que c est aléatoire uniformément réparti entre
0 et v−1) ; puis on effectue une boucle : à chaque étape,
on génère un bit aléatoire b (valant 0 ou 1 avec
probabilité ½ pour chacun, et indépendant de tous les autres, donc) et
on remplace v ← 2v
et c ← 2c+b ; puis on
compare v avec n
et c avec n : si v<n
(ce qui implique forcément c<n) on continue
simplement la boucle (il n'y a pas assez d'entropie) ;
si v≥n et c<n, on
termine l'algorithme en renvoyant la valeur c ; enfin,
si c≥n, on
effectue v ← v−n
et c ← c−n et on continue la boucle.
Le calcul du nombre moyen de tirages effectués est fastidieux
(voir cette référence citée dans le fil Twitter référencée
ci-dessus), mais c'est optimal.
L'algorithme que je viens de décrire s'adapte assez bien pour tirer un entier uniforme entre 0 et n−1 à partir d'une source de entiers uniformes entre 0 et m−1 (le cas que je viens de décrire est le cas m=2), autrement dit : comment fabriquer un dé à n faces à partir d'un dé à m faces ? Je n'ai pas vraiment envie de réfléchir à si c'est optimal (mise à jour : on l'a fait pour moi), mais c'est en tout cas assez élégant : on utilise deux variables internes à l'algorithme, notées v et c, qu'on initialise par v←1 et c←0 ; puis on effectue une boucle : à chaque étape, on génère un tirage aléatoire b entre 0 et m−1 à partir de la source dont on dispose et on remplace v ← m·v et c ← m·c+b ; puis on effectue la division euclidienne de v et de c par n : si les deux quotients calculés sont différents (⌊c/n⌋ < ⌊v/n⌋), on termine l'algorithme en renvoyant le reste c%n := c − n·⌊c/n⌋ de la division de c par n, tandis que si les deux quotients sont égaux, on remplace chacun par son reste, c'est-à-dire v ← v%n et c ← c%n et on continue la boucle.
Introduisons maintenant aussi des tirages continus.
Comment tirer uniformément un point dans le simplexe de dimension d ? Autrement dit, on veut obtenir d+1 réels positifs (z0,…,zd) de somme 1 uniformément répartis sur toutes les possibilités : ce qui revient à dire que disons, les d premiers sont uniformément répartis sur toutes les possibilités, ce qui est cette fois un conditionnement simple : il y a donc un algorithme trivial mais trop inefficace qui consiste à tirer z0,…,zd−1 aléatoirement (indépendamment et uniformément) entre 0 et 1, poser zd := 1 − z0 − ⋯ − zd−1, tester si zd≥0 et recommencer tous les tirages si ce n'est pas le cas ; ceci fonctionne, mais demandera d! essais en moyenne, ce qui n'est pas acceptable.
Voici un algorithme qui marche et qui est efficace : tirer d nombres réels aléatoirement (indépendamment et uniformément) entre 0 et 1, trier ces nombres en ordre croissant, disons qu'on les appelle t1≤⋯≤td, poser t0=0 et td+1=1, et définir zi := ti+1 − ti.
En revanche, que je sache, si on veut tirer un
point
du polytope
de Birkhoff, c'est-à-dire une matrice carrée à coefficients
positifs dont toutes les lignes et toutes les colonnes ont somme 1
(matrice bistochastique
), il n'y a pas d'algorithme intelligent
connu (seulement des méthodes approchées). Ceci illustre le fait que
ce n'est pas parce qu'un ensemble est facile à décrire qu'il est
facile de tirer au hasard uniformément dedans.
Une autre façon de tirer un point dans le simplexe (continu, i.e., réel) de dimension d consiste à tirer d+1 réels positifs selon une loi exponentielle (cf. ci-dessous), toujours la même, mais peu importe son paramètre, et ensuite les diviser par leur somme (i.e., normaliser pour avoir une somme 1).
Ayant parlé du simplexe continu, je ne peux pas ne pas évoquer le simplexe discret et la question de comment tirer uniformément un point dans le simplexe discret de dimension d et de taille N, autrement dit, d+1 entiers positifs (a0,…,ad) de somme N uniformément répartis sur toutes les possibilités. L'adaptation évidente de l'algorithme avec tri que j'ai donné ne marche pas (il faudrait tirer d entiers entre 0 et N non pas indépendamment avant de les trier, mais selon une urne de Pólya ; comme ce n'est pas évident, le plus simple est de faire la chose suivante) : à la place, on va plutôt tirer d entiers distincts entre 1 et N+d (il y a toutes sortes de manières de faire ça, ce n'est pas compliqué et je ne rentre pas dans les détails, le plus évident est simplement de tirer chaque nombre en recommençant tant qu'on tombe sur un nombre déjà tiré), les trier en v1<⋯<vd, poser v0=0 et vd+1=N+d+1, et définir ai := vi+1 − vi − 1. Mais une façon différente et possiblement intéressante dans certains cas d'obtenir ce tirage est de commencer par tirer (z0,…,zd) dans le simplexe (réel) comme je viens de le dire, puis de simuler une élection à la proportionnelle à la plus forte moyenne avec la méthode d'Hondt (concrètement, cela signifie qu'on commence avec (a0,…,ad) ← (0,…,0), puis on trouve le i qui maximise zi/(ai+1), et on l'incrémente, ai ← ai+1, et on répète ça N fois ; la preuve du fait que ceci donne bien une répartition uniforme sur les (a0,…,ad) est donnée ici). Bon, je ne sais pas dans quel cas cet algorithme bizarre pourrait servir (peut-être si on ne connaît pas N à l'avance ?), mais il est tellement remarquable que je ne peux pas ne pas le mentionner.
Comment tirer un réel selon une loi
exponentielle ? Pour en avoir une d'espérance 1, il suffit
de prendre l'opposé du logarithme (naturel) d'un réel aléatoire
uniforme entre 0 et 1. (Ceci pose un problème de précision si le
nombre aléatoire uniforme est trop proche de 0 : je ne veux pas
rentrer dans ces questions-là, mais une façon de corriger est que si
le réel uniforme est, disons, <1/N avec N
fixé à l'avance, par exemple une puissance de 2, on le jette et on
tire un autre réel aléatoire uniforme qu'on divise ensuite
par N : et on fera bien sûr ça de façon récursive,
c'est-à-dire que la fonction calculer un réel uniforme entre 0 et 1
à résolution améliorée proche de 0
commence par tirer un réel
uniforme entre 0 et 1 avec une résolution standard, et s'il est
<1/N, le jette, s'applique elle-même de façon récursive
et renvoie le résultat divisé par N : ceci peut évidemment
être rendu itératif de façon facile.) Pour une loi exponentielle
d'espérance a, il suffit bien sûr simplement de multiplier
par a.
Et pour un entier distribué selon une loi géométrique (exponentielle discrète) ? Si p (la probabilité de succès) n'est pas excessivement petit, le plus simple à programmer est simplement l'algorithme naïf consistant à tirer des variables de Bernoulli de paramètre p (indépendantes) de façon répétée jusqu'à obtenir un succès, et renvoyer le nombre d'échecs ; on peut aussi, pour le même résultat, tirer un réel selon une loi exponentielle d'espérance −1/log(1−p) et prendre sa partie entière : ceci produit une variable de loi géométrique d'espérance (1−p)/p.
Comment tirer un entier distribué selon une loi de Poisson ? Si l'espérance λ désirée n'est pas excessivement grande, le plus simple à programmer est simplement l'algorithme naïf consistant à simuler un processus de Poisson : tirer des variables exponentielles d'espérance 1 (indépendantes) de façon répétée jusqu'à ce que leur somme dépasse λ, et renvoyer le nombre de variables tirées moins 1 ; d'après la manière dont j'ai décrit la façon de tirer des variables aléatoires exponentielles, cela revient aussi à tirer des variables aléatoires uniformes entre 0 et 1 (indépendantes) de façon répétée jusqu'à ce que leur produit passe en-dessous de exp(−λ) (et toujours renvoyer le nombre de variables tirées moins 1).
Je rappelle qu'un processus de Poisson (homogène) s'obtient de la façon suivante : partir de 0 et sommer de façon répétée des variables exponentielles indépendantes (de même espérance, l'inverse de la densité voulue) pour obtenir les valeurs positives du processus, et faire la même chose avec un signe moins pour les négatives, puis oublier la valeur zéro (ça ressemble à une erreur, parce qu'on se dit qu'on crée un trou plus grand en zéro, pourtant c'est bien correct : l'espérance de la plus petite valeur positive moins la plus grande valeur négative est bien le double de l'espérance de deux valeurs consécutives, c'est le fameux paradoxe du temps d'attente du bus).
Pour une loi binomiale, si le nombre n
d'essais n'est pas excessivement grand, le plus simple à programmer
est simplement l'algorithme naïf consistant à tirer n
variables de Bernoulli de paramètre p (indépendantes) et de
compter combien d'entre elles ont donné vrai. Pour une loi
binomiale négative, si le nombre r d'échecs avant
arrêt est entier et n'est pas excessivement grand, de même,
le plus simple à programmer est simplement l'algorithme naïf
consistant à tirer r variables géométriques (indépendantes)
de probabilité de succès (enfin, ici c'est plutôt un échec, mais la
terminologie la plus courante parle de probabilité de succès
)
1−p et de les sommer ; en revanche, je ne sais pas générer
une loi binomiale négative avec r non entier sans faire
appel à des calculs de fonction Γ ou des systèmes de rejets, mais
pour r demi-entier, ce qui est un cas important, je vais en
redire un mot ci-dessous.
Bon, la plus important distribution continue est sans doute
la loi gaussienne. Je rappelle que la somme
des carrés de r variable aléatoire gaussiennes
indépendantes de moyenne 0 et de même espérance définit une loi du χ²
(chi-carré) avec r degrés de liberté
, et que c'est
un cas particulier de la loi Γ (gamma), à savoir celui où le paramètre
de forme
est demi-entier, précisément r/2 ; en
particulier, la somme des carrés de deux (r=2)
variables aléatoires gaussiennes indépendantes centrées de même
espérance est une variable aléatoire exponentielle : ou, si on
préfère, une gaussienne (isotrope, centrée) de
dimension r=2 (i.e., une variable aléatoire dans ℝ² dont
les deux coordonnées sont des gaussiennes réelles indépendantes
centrées de même espérance), si on la regarde en coordonnées polaires,
a une composante radiale qui est la racine carrée d'une variable
aléatoire exponentielle (et la composante angulaire est, bien sûr,
uniforme).
L'observation du paragraphe précédent donne naissance à l'algorithme de Box-Muller pour tirer une variable de loi gaussienne : tirer U et V deux variables uniformes indépendantes sur [0;1], et alors √(−2·log(U))·cos(2πV) et √(−2·log(U))·sin(2πV) fournissent deux variables gaussiennes standard (i.e., de moyenne 0 et variance 1) indépendantes. Ici, −log(U) sert à tirer une variable exponentielle d'espérance 1, qu'on peut éventuellement corriger comme je l'ai expliqué plus haut.
Un des intérêts de la loi gaussienne est, justement, qu'elle permet de fabriquer des lois gaussiennes en toutes dimensions. Ceci répond notamment à la question fondamentale suivante : comment tirer uniformément un point sur la sphère de dimension d ? Le plus simple est de tirer d+1 variables gaussiennes indépendantes centrées de même espérance, ce qui constitue collectivement une variable gaussienne centrée isotrope dans ℝd+1, et de normaliser cette dernière, c'est-à-dire de tout diviser par la norme (la racine carrée de la somme des carrés des variables tirées). (Pour d=1, i.e., pour le cas d'un cercle, cela revient bien sûr simplement à défaire ce qu'on a fait ci-dessus, et on peut court-circuiter ce bazar en calculant simplement (cos(2πV), sin(2πV)) avec V uniforme dans [0;1].)
(En rapport avec ce que je viens de dire, je devrais mentionner le fait notable suivant que si (z0:⋯:zd) est un point aléatoire uniforme de l'espace projectif complexe de dimension (complexe) d pour la mesure associée à la métrique de Fubini-Study, ce qui s'obtient simplement en prenant un point de la sphère de dimension 2d+1, ses 2d+2 coordonnées réelles étant interprétées comme les d+1 composantes réelles et imaginaires des coordonnées complexes, alors une fois normalisé par |z0|² + ⋯ + |zd|² = 1, ce qui est déjà le cas si on a fait comme je viens de le dire, le point (|z0|², …, |zd|²) est uniformément réparti sur le simplexe de dimension d. Ce n'est qu'une reformulation du fait que pour tirer un point de la (2d+1)-sphère selon l'algorithme de Box-Muller, on tire un point du simplexe de dimension d par normalisation de d+1 variables aléatoires exponentielles dinépendantes et on tire d+1 phases uniformes indépendantes. Mais ceci suggère la question suivante : si on tire une matrice unitaire (d+1)×(d+1) uniforme (selon la mesure de Haar de Ud+1) et qu'on fabrique d+1 points du simplexe de dimension d en considérant les carrés des modules des différentes entrées de la matrice, quelle est la loi jointe ainsi obtenue ?)
Mais il faut noter qu'il n'y a pas que la sphère qu'on peut traiter ainsi : pour tirer uniformément une matrice orthogonale de taille m×m (selon la mesure de Haar de SOm, si on veut), le plus simple est de tirer m² variables gaussiennes indépendantes centrées de même espérance, qu'on considérera comme m vecteurs de taille m, et de leur appliquer l'algorithme de Gram-Schmidt pour les transformer en une base orthonormée (qui est uniformément répartie sur toutes les bases orthonormées puisque, justement, toute l'opération est invariante sous l'effet du groupe orthogonal). Il me semble vaguement qu'on peut faire quelque chose du même goût pour tous les groupes de Lie réels compacts, mais je ne sais plus bien d'où je tire cette idée.
(À titre d'exemple, pour tirer un élément g de G₂, on commence par tirer deux octonions imaginaires purs gaussiens indépendants, qu'on orthonormalise, appelons les g(i) et g(j), qui sont donc orthogonaux et chacun uniforme sur la 6-sphère, on pose g(k) = g(i)·g(j) où il s'agit ici de la multiplication octonionique, il est orthogonal à g(i) et g(j), puis on tire un troisième octonion imaginaire pur gaussien, indépendant des tirages précédents,, dont on retire, comme par Gram-Schmidt, les composantes selon g(i), g(j) et g(k), appelons g(ℓ) le résultat : on pose alors g(i·ℓ) = g(i)·g(ℓ) et de même pour g(j·ℓ) et g(k·ℓ) : alors 1, g(i), g(j), g(k), g(ℓ), g(i·ℓ), g(j·ℓ) et g(k·ℓ) non seulement sont une base orthonormée mais l'application g de passage de la base 1, i, j, k, ℓ, i·ℓ, j·ℓ, k·ℓ à elle est uniformément distribuée selon la mesure de Haar de G₂.)
Bien sûr, pour tirer uniformément un point dans la boule unité de dimension d, on pourra tirer un point de la sphère unité de dimension d−1 (on vient d'expliquer comment) et le multiplier par la puissance (1/d)-ième d'une variable uniforme sur [0;1]. La méthode de rejet (tirer un point de [−1;1]d et recommencer jusqu'à ce qu'il soit dans la boule) devient extrêmement inefficace quand d est grand ; notons cependant qu'elle ne l'est pas pour d petit, et pour d=2 il est possiblement plus efficace de procéder par rejet pour tirer un point du disque que comme je viens de dire, voire, plus efficace pour tirer un point du cercle de tirer d'abord un point du disque en procédant par rejet et ensuite de le normaliser, en évitant ainsi d'avoir à évaluer des lignes trigonométriques (ceci peut, justement, servir dans l'implémentation de Box-Muller).
Le fait de savoir tirer une variable gaussienne permet de tirer des
variables de loi Γ de forme demi-entière, comme je l'ai déjà signalé
(puisque c'est une loi du χ²). Il est notable, même si je ne sais pas
si ça peut servir à quelque chose, que quelque chose d'analogue peut
aussi servir dans le cas discret : ① pour tirer une
variable k entre 0 et n (inclus) selon une loi
bêta-binomiale de paramètres α=β=½, par la
définition même de la loi bêta-binomiale, on peut tirer p
comme cos(2πV) avec V uniforme sur [0;1] (ceci
donne p suivant une loi bêta de
paramètres α=β=½) et ensuite k entre
0 et n (inclus) selon une loi binomiale de
paramètre p ; et ② pour tirer une variable distribuée selon
une loi de Pólya (binomiale négative) avec r=½ (ce que j'ai
envie d'appeler une loi demi-géométrique
, parce que la somme de
deux variables indépendantes ainsi distribuées suit une loi
géométrique), on commence par tirer n selon une loi
géométrique (pour le même p), puis on tire k
selon une loi bêta-binomiale avec α=β=½ comme je
viens de le dire.
(Cf. cette
question. Tout ça n'a rien de spécifique à ½ : c'est juste que ½
est le cas analogue discret du carré de la loi gaussienne, et aussi le
seul cas non entier que je sais traiter de façon élégante. Soit dit
en passant, je suis preneur de toute situation « réelle » où la loi de
Pólya de paramètre r=½ ou bien la loi bêta-binomiale
avec α=β=½ joue un rôle particulier.)
Ajout () : On me propose en commentaire le problème pratique suivant (tout à fait dans l'esprit de ce billet de blog) : comment répartir un ensemble de n personnes en deux « équipes », A de k personnes et B de n−k personnes, avec juste une pièce pour tirer à pile ou face ? Je propose la réponse suivante : chacun tire une pièce (disons qu'on écrit ‘0’ pour pile et ‘1’ pour face) :
- s'il y a k personnes ayant tiré 0 et n−k ayant tiré 1, ce sera les équipes A et B recherchés respectivement, et on s'arrête là ;
- s'il y a ℓ > k personnes ayant tiré 0, alors on met d'ores et déjà les n−ℓ < n−k ayant tiré 1 dans l'équipe B, et on répartit les ℓ personnes ayant tiré 0 en une équipe de k, qui sera l'équipe A, et une équipe de ℓ−k, qu'on ajoutera à l'équipe B, en recommençant (récursivement) la procédure qu'on est en train de définir ;
- symétriquement, s'il y a ℓ < k personnes ayant tiré 0, alors on les met d'ores et déjà dans l'équipe A, et on répartit les n−ℓ > n−k ayant tiré 1 en une équipe de k−ℓ, qu'on ajoutera à l'équipe A, et une équipe de n−k, qui sera l'équipe B, en recommençant (récursivement) la procédure qu'on est en train de définir.
Autrement dit, tout le monde tire une première pièce, ce qui donne deux groupes, 0 et 1, qu'on aurait envie d'identifier aux équipes A et B mais ils n'ont pas forcément la bonne taille : du coup, le groupe qui est moins nombreux que l'équipe qu'on voulait constituer va directement dans l'équipe en question, et le groupe trop gros fait un nouveau tirage pour se répartir les places restant à répartir selon le même procédé, jusqu'à ce que tout le monde soit réparti. Il revient au même de dire que chacun tire de façon répétée des bits aléatoires mais qu'on fait juste les tirages suffisants pour pouvoir envoyer les k premiers dans l'équipe A et les n−k derniers dans l'équipe B.
C'est un joli problème, et il me semble que la solution que je propose est tout à fait facile à mener en pratique (elle est plus compliquée à décrire qu'à appliquer), et ça ne m'étonnerait pas qu'elle soit optimale en un certain sens (j'ai la flemme d'y réfléchir ; mais certainement pas par l'espérance du nombre total de lancers puisque pour k=1 c'est beaucoup moins efficace, selon cette métrique, que de générer un nombre uniforme entre 0 et n−1 par la méthode que j'ai décrite plus haut).