simulation de variables aléatoires de Xenakis

La probabilité que $X$, variable aléatoire de Xenakis, soit comprise entre $a$ et $b$ est $\int_{a}^{b} (2-2x) dx$. Ce qui veut dire que pour calculer cette probabilité, on a besoin d’une primitive de $2-2x$ : La probabilité que $X$ soit inférieure à $x$, qui est $F(x)=2x-x^2$. Finalement, la probabilité que $a\leqslant X \leqslant b$ est $F(b)-F(a)=2b-b^2-(2a-a^2)=2(b-a)-(b^2-a^2)=2(b-a)-(b-a)(b+a)=(b-a)(2-a-b)$.

On constate que la primitive est plus facile à calculer que celle pour une loi exponentielle tout en étant plus intéressante que celle pour une loi uniforme. Ce qui motive déjà dès le début l’étude de cette variable aléatoire en Terminale...

Par analogie avec le cours de Première, on définit en Terminale, l’espérance d’une variable aléatoire $X$ de densité $f$ comme étant l’intégrale de $x \times f(x)$ sur le support de $f$. Le problème est que la seule variable aléatoire de support fini au programme est la loi uniforme pour laquelle l’intégrale est un calcul d’aire de trapèze, ce qui rend inintéressant (parce que non nécessaire) le calcul de primitive. Comme les intégrales indéfinies ne sont pas au programme de lycée, on aborde rapidement cette question parce qu’on aborde quelque chose qui n’est pas au programme. De plus, il faut une intégration par parties (ou un logiciel de calcul formel) pour avoir une primitive de $x e^{-\lambda x}$.

Ces problèmes n’existent pas pour une variable $X$ de Xenakis, pour laquelle l’espérance est donc $\int_{0}^{1} x(2-2x) dx = \int_{0}^{1} (2x-2x^2) dx$ : Une primitive du trinôme à intégrer est $G(x)=x^2-\frac{2x^3}{3}$ ; cette primitive s’annulant en 0, l’espérance est donc $G(1)=1^2-\frac{2\times 1^3}{3}=1-\frac{2}{3}=\frac{1}{3}$.

On constate que l’espérance est légèrement supérieure à la médiane $1-\frac{\sqrt{2}}{2}$, ce qui montre que la densité est asymétrique. D’ailleurs sa "skewness" (ou asymétrie (statistiques)) est $\frac{2\sqrt{2}}{5}$ soit environ 0,56 qui est assez élevée.

Il y a deux manières de calculer la variance, selon la définition vue en Première :

1. Ou bien on a vu que $V(X)=E(X-EX)^2$, ce qui ici donne $V(x)=\int_{0}^{1}(x-\frac{1}{3})^2 (2-2x) dx=\int_{0}^{1}(x^2-\frac{2}{3}x+\frac{1}{9})(2-2x) dx$ soit, après développement, $V(X)=\int_{0}^{1}(-2x^3+\frac{10}{3}x^2-\frac{14}{9}x+\frac{2}{9}) dx=[-\frac{x^4}{2}+\frac{10 x^3}{9}-\frac{7 x^2}{9}+\frac{2x}{9}]_{0}^1$ soit, puisque la primitive choisie s’annule en 0, $V(X)=-\frac{1}{2}+\frac{10}{9}-\frac{7}{9}+\frac{2}{9}=\frac{1}{18}$
2. Ou alors on utilise la formule $V(x)=E(X^2)-(EX)^2$ en commençant par calculer l’espérance de $X^2$, soit $\int_{0}^{1}x^2(2-2x) dx=\int_{0}^{1}(2x^2-2x^3) dx$. Comme une primitive de $2x^2-2x^3$ est $\frac{2x^3}{3}-\frac{x^4}{2}$ qui a la gentillesse de s’annuler en 0, on trouve finalement $E(X^2)=\frac{2}{3}-\frac{1}{2}=\frac{1}{6}$. Alors $V(X)=\frac{1}{6}-\frac{1}{9}=\frac{1}{18}$.

On constate que la deuxième version est plus facile à calculer que la première, ce qui montre son intérêt. De toute manière, il est possible, comme dans les sujets de bac, de donner la primitive dans l’énoncé, ce qui amène l’exercice à un calcul de dérivée et un développement.

Il résulte donc de ce qui précède que l’écart-type d’une variable aléatoire de Xenakis est $\sqrt{VX}=\sqrt{\frac{1}{18}}=\frac{\sqrt{2}}{6}$ soit environ 0,24.

Le calcul de primitive effectué plus haut montre que, si $X$ est une variable aléatoire de Xenakis de paramètre 1, on a $P(X \leqslant y)=2y-y^2$ pour tout $y \in [0 ;1]$. En posant $X=f(R)$ où $R$ est une variable (pseudo-)aléatoire uniforme comme celle de la calculatrice et $f$ une fonction strictement croissante, on réécrit la définition de la fonction de répartition ci-dessus comme $P(f(R) \leqslant y)=2y-y^2$, soit en posant $y=f(x)$, $P(R \leqslant x)=2 f(x)-f(x)^2$ ; or pour une variable aléatoire uniforme entre 0 et 1, la probabilité qu’elle soit inférieure à $x$ est égale à $x$ ; d’où $x=2y-y^2$. Pour trouver l’expression de $f$, on résout cette équation du second degré, qui s’écrit aussi $-y^2+2y-x=0$ : Son discriminant est $\Delta=2^2-4 (-1)(-x)=4-4x=4(1-x)$ dont la racine carrée est $\sqrt{\Delta}=2\sqrt{1-x}$ ; or une seule des deux solutions de l’équation est inférieure à 1, c’est $1-\sqrt{1-x}$ ; donc $f(x)=1-\sqrt{1-x}$. Pour simuler une variable aléatoire de Xenakis, on peut donc appliquer cette fonction à une variable aléatoire uniforme.

Si on est confronté à une variable de Xenakis de paramètre $a$ inconnu, on peut estimer le paramètre à partir d’un échantillon suffisamment grand de valeurs de $X$. Plusieurs estimateurs peuvent être envisagés :

• le maximum sur un échantillon assez vaste de réalisations de $X$ ; il semble être biaisé (d’ailleurs la probabilité d’avoir une valeur proche du maximum est faible).
• le produit de la médiane par $2+\sqrt{2}$ sur un échantillon assez vaste de réalisations de $X$ ;
• le triple de la moyenne sur échantillon. Il s’agit d’un cas particulier de méthode des moments (statistiques).

Dans son livre, il est clair que cette variable aléatoire est utilisée par Xenakis pour approcher une variable exponentielle en la bornant. Dans le cas du paramètre 1, quel est le paramètre de cette loi exponentielle ?

Pour y répondre, on peut chercher à minimiser l’aire coloriée entre les deux courbes représentant les fonctions de répartition :

Mais depuis Hilbert, on préfère minimiser l’intégrale du carré de la différence. Celle-ci est composée de deux morceaux :

• $\int_{0}^{1}(1-e^{-\lambda x}-2x+x^2)^2 dx=\frac{1}{5}-\frac{3}{2 \lambda}+\frac{4}{\lambda ^2}-\frac{4}{\lambda ^3}-\frac{e^{-2 \lambda}}{2 \lambda}+\frac{4 e^{-\lambda}}{\lambda ^3}$
• $\int_{1}^{+ \infty}(1-e^{-\lambda x}-1)^2 dx=\int_{1}^{+ \infty} e^{-2\lambda x}dx=\frac{e^{-2 \lambda}}{2 \lambda}$

Donc on veut minimiser la fonction $\lambda \mapsto \frac{1}{5}-\frac{3}{2 \lambda}+\frac{4}{\lambda ^2}-\frac{4}{\lambda ^3}-\frac{e^{-2 \lambda}}{2 \lambda}+\frac{4 e^{-\lambda}}{\lambda ^3}+\frac{e^{-2 \lambda}}{2 \lambda}$ ; sa dérivée faisant intervenir des produits d’exponentielles par $\lambda$ et ses puissances, on ne peut que chercher numériquement les valeurs qui l’annulent. D’après Xcas, on trouve environ 2,83173397154 qui est supérieur à la valeur de e (nombre) utilisée par Xenakis.

Pour simuler la longueur d’un segment aléatoire, on manipule directement des listes, c’est ce que permet R (langage de programmation et environnement statistique).

• La première méthode consiste à créer des listes de nombres aléatoires uniformes appelées A et B (car ce sont des points) puis à calculer la distance X comme valeur absolue de leur différence :

N <- 100000
A <- runif(N,0,1)
B <- runif(N,0,1)
X <- abs(A-B)
bornes <- seq(0,1,0.05)
hist(X,breaks=bornes)

L’histogramme obtenu est celui-ci :

• Mais R n’est pas un logiciel de géométrie, alors autant faire une simulation à partir de la fonction de répartition :

N <- 100000
U <- runif(N,0,1)
X <- 1-sqrt(U)
bornes <- seq(0,1,0.05)
hist(X,breaks=bornes)

L’histogramme obtenu est si proche du précédent qu’il ne vaut pas la peine de l’afficher

• Et pour simuler X comme minimum de deux variables aléatoires uniformes, on peut effectuer une boucle :
  • On crée une liste vide avec X <- c()
  • On boucle sur n allant de 1 à N ; à chaque passage dans la boucle,
    • on crée une liste de deux variables aléatoires uniformes avec runif(2,0,1)
    • on cherche la plus petite des deux valeurs avec la fonction statistique min(runif(2,0,1))
    • on concatène le résultat à X avec c(X,min(runif(2,0,1))) que l’on met dans X

N <- 10000
X <- c()
for (n in 1:N) {
  X <- c(X,min(runif(2,0,1)))
}
bornes <- seq(0,1,0.05)
hist(X,breaks=bornes)

Là encore, l’histogramme est assez proche des deux précédents, la seule différence étant que le script s’exécute plus lentement que les deux autres, ce qui a conduit à choisir N plus petit :

Comme la loi de Xenakis est un cas particulier de loi bêta et que celle-ci est connue de R, on peut simplifier et accélerer ce qui a été vu ci-dessus, en effet il suffit pour simuler une variable de Xenakis, de faire rbeta(x,1,2) où x est une liste de nombres de la taille voulue :

N <- 10000
x <- seq(0,1,len=N+1)
X <- rbeta(x,1,2)
bornes <- seq(0,1,0.05)
hist(X,breaks=bornes,col='lightcyan')
points(x,N*0.05*dbeta(x,1,2),type="l",col="red")

En prime, on a avec dbeta(x,1,2) la fonction de densité, et avec pbeta(x,1,2) la fonction de répartition de la variable de Xenakis. Le graphique obtenu est celui-ci :

L’idée de base était la suivante : Si 10 élèves réussissent une simulation sur 1000 expériences chacun, en moyennant leurs données, on aura l’équivalent d’une simulation sur 10000 données, et en plus, si on transmet les données au serveur sous forme de points, l’alignement des points du nuage moyen sera meilleur que sur chacun des jeux de données. De ces deux points de vue le TP a été décevant :

Au moment de démarrer le serveur, on peut choisir d’activer, ou non, le travail collaboratif :

Dans le premier groupe, ce choix n’ayant pas été fait, chaque figure était envoyée par l’élève l’ayant créée, dans un onglet qui lui était propre, et il n’a pas été possible de saisir les données des différents onglets par script pour calculer la moyenne [1]. Du moins, pas pendant le TP. Le nuage moyen, en rouge, calculé sur 7 des fichiers, n’est pas tellement aligné que ça :

Dans le second groupe, par contre, l’option "travail collaboratif" a été actionnée, et les nuages de points, portant chacun des noms différents, ont bien été envoyés dans la figure "global". Par contre certains élèves avaient laissé le nombre maximal n à 100 au lieu de 1000 et le calcul de la moyenne a été fait en tenant compte de cela :

noms=["Jirayalecochon","Mamdex","Corobizar","Drank","Kakashi","Gros tête","Walter","Atomiiz","Prototype"];
for(p=0; p<=10;p++){
	somme = 0;
	for (n=0; n<=6; n++){
		somme += Y(noms[n]+p);
	}
	for (n=7; n<=8; n++){
		somme += Y(noms[n]+p)/10;
	}
	q = Point("Q"+p,p,somme/9);
	SetPointType(q,"cross");
	SetColor(q,"red");
}

Comme les élèves ont vu qu’en bougeant un point sur leur écran, le point bougeait sur tous les écrans, et que ça les a amusés, ils ont bougé beaucoup de points et le nuage moyen (en rouge ci-dessous) n’est pas très aligné non plus :

Il aurait donc fallu rajouter cette ligne dans le script :

setFixed(p,true);

Il aurait également fallu faire faire des nuages de 100000 points [2] ; mais avec 1000 points, le TP dure 40 minutes...

Une fois faite l’activité sur la loi de Xenakis, il est facile de refaire tout avec la loi exponentielle, même avant le cours sur celle-ci. Toutefois, la sensibilité aux erreurs des élèves rend le TP long (environ une heure) et peu probant. Aussi semble-t-il préférable de faire le TP avec la calculatrice, en entrant

Suite(-10*ln(nbrAléat),I,1,100) →L1

puis en affichant l’histogramme (surtout avec des calculatrices comme la Ti84 qui peuvent aller jusqu’à 1000 données) ; du coup il reste du temps pour la très motivante question : "Quelle fonction passe au plus près de l’histogramme ?". Ce qui a permis de conjecturer une fonction du type "constante fois l’exponentielle de 0,1X" plus ou moins spontanément pour environ le tiers de la classe.